# Building a function to calculate the series solution to this Solved Boundary value problem

I have the following code which solves for $$(1)$$ (i.e. solves for $$C_1$$ and $$C_2$$)

$$T(x,y,z)=\sum_{n,m=0}^{\infty}(C_1 e^{\gamma z}+C_2 e^{-\gamma z})\sin\bigg(\frac{\alpha_n x}{L}+\beta_n\bigg)\sin\bigg(\frac{\delta_m y}{l}+\theta_m\bigg)+T_a \tag 1$$

T[x_, y_, z_] = (C1* E^(γ z) + C2 E^(- γ z))*Sin[(α x/L) + β]*Sin[(δ y/l) + θ] + Ta
tc[x_, y_] = E^(-NTUC* y/l)*{tci + (NTUC/l)*Integrate[E^(NTUC*s/l)*T[x, s, 0], {s, 0, y}]};
tc[x_, y_] = tc[x, y][[1]];
bc1 = (D[T[x, y, z], z] /. z -> 0) == pc (T[x, y, 0] - tc[x, y]);
ortheq1 = Integrate[bc1[[1]]*Sin[(α x/L) + β]*Sin[(δ y/l) + θ], {x, 0, L}, {y, 0, l}] == Integrate[bc1[[2]]*Sin[(α x/L) + β]*Sin[(δ y/l) + θ], {x, 0, L}, {y, 0, l}];
ortheq1 = ortheq1 // Simplify;
th[x_, y_] = E^(-NTUH*x/L)*{thi + (NTUH/L)*Integrate[E^(NTUH*s/L)*T[s, y, w], {s, 0, x}]};
th[x_, y_] = th[x, y][[1]];
bc2 = (D[T[x, y, z], z] /. z -> w) == ph (th[x, y] - T[x, y, w]);
ortheq2 = Integrate[bc2[[1]]*Sin[(α x/L) + β]*Sin[(δ y/l) + θ], {x, 0, L}, {y, 0, l}] == Integrate[bc2[[2]]*Sin[(α x/L) + β]*Sin[(δ y/l) + θ], {x, 0, L}, {y, 0, l}];
ortheq2 = ortheq2 // Simplify;
soln = Solve[{ortheq1, ortheq2}, {Subscript[C, 1], Subscript[C, 2]}];
CC1 = C1 /. soln[[1, 1]];
CC2 = C2 /. soln[[1, 2]];
expression1 := CC1;
c1[α_, β_, δ_, θ_, γ_] := Evaluate[expression1];
expression2 := CC2;
c2[α_, β_, δ_, θ_, γ_] := Evaluate[expression2];


The following relations hold, $$\beta_n=\tan^{-1}(1.66\times10^4 \alpha_n)$$ and $$\delta_m=\tan^{-1}(8.33\times10^3 \theta_m)$$

The n=0 values is $$\alpha_0=0.01095$$ and m=0 value is $$\delta_0=0.01549$$.

Subsequently from n=1 and m=1 it is known that $$\alpha_n=n\pi$$ and $$\delta_m=m\pi$$.

I want to build a function such that this summation can be automatically performed for the desired values of $$n$$ and $$m$$.

$$T_a$$ is added only once in the final $$T(x,y,z)$$. The rest of the constants along with the other functions I wish to calculate are given below:

L = 0.9; l = 1.8; w = 0.0003; NTUH = 17.394; NNTUC = 22.151; ph = 8.6; pc = 13.93;
γ = Sqrt[(α/L)^2 + (δ/l)^2];
thi=460;tci=300;Ta=380;
tc1[x_, y_] = E^(-NTUC* y/l)*{tci + (NTUC*/l)*Integrate[E^(NTUC* s/l)*(TWnet /. {y -> s, z -> 0}), {s, 0, y}]};
th1[x_, y_] = E^(-NTUH* x/L)*{thi + (NTUH/L)*Integrate[E^(NTUH* s/L)*(TWnet /. {x -> s, z -> w}), {s, 0, x}]};
Plot[tc1[x, l], {x, 0, L}]
Plot[th1[L, y], {y, 0, l}]
THotAvg = Integrate[th1[x, y]/l, {y, 0, l}];
TColdAvg = Integrate[tc1[x, y]/L, {x, 0, L}];
THotAvg /. x -> L
TColdAvg /. y -> l
Plot[THotAvg, {x, 0, L}]
Plot[TColdAvg, {y, 0, l}]


The term TWnet in the above code section is the final $$T(x,y,z)$$ function I desire. So if someone can make the final distribution as a function then terms like TWnet /. {y -> s, z -> 0} would be something like TWnet[x,s,0]

I hope I was able to clearly explain the requirements here.

NOTE: The first code section takes some time to execute

CONTEXTUAL INFORMATION

I am trying to solve $$\nabla^2 T(x,y,z)=0$$ defined on $$x\in[0,L], y\in[0,l]$$ and $$z\in[0,w]$$ subjected to the following boundary conditions:

$$k(\frac{\partial T(0,y,z)}{\partial x})=h_a(T(0,y,z)-T_a) \tag A$$

$$-k(\frac{\partial T(L,y,z)}{\partial x})=h_a(T(L,y,z)-T_a) \tag B$$

$$k(\frac{\partial T(x,0,z)}{\partial y})=h_a(T(x,0,z)-T_a)\tag C$$

$$-k(\frac{\partial T(x,l,z)}{\partial y})=h_a(T(x,l,z)-T_a) \tag D$$

$$\frac{\partial T(x,y,0)}{\partial z} = p_c\bigg(T(x,y,0)-e^{-\beta_c y/l}\left[t_{ci} + \frac{\beta_c}{l}\int_0^y e^{\beta_c s/l}T(x,s,0)ds\right]\bigg) \tag E$$

$$\frac{\partial T(x,y,w)}{\partial z} = p_h\bigg(e^{-\beta_h x/L}\left[t_{hi} + \frac{\beta_h}{L}\int_0^x e^{\beta_h s/L}T(x,s,w)ds\right]-T(x,y,w)\bigg) \tag F$$

Now under the conditions $$A,B,C,D$$, the solution form of the three-dimensional Laplacian is given by $$(1)$$

$$\gamma=\sqrt{(\alpha/L)^2 + (\delta/L)^2}$$ (Have not mentioned this explicitly in the original question, so I wrote it here).

In the first section of the code I apply the $$z$$ boundary conditions and use orthogonality to determine the constants $$C_1, C_2$$. I must mention here that I have already proven the orthogonality of $$\sin\bigg(\frac{\alpha_n x}{L}+\beta_n\bigg)$$ under the boundary conditions $$A-D$$ The values of $$\alpha$$ and $$\beta$$ are to be calculated using the following transcendental equation:

$$2\cot{\alpha}=\frac{k\alpha}{h_a L}-\frac{h_aL}{k\alpha}\tag G$$ $$\beta=\tan^{-1}(\frac{k \alpha}{h_a L})\tag H$$

Similar set of equation exists for $$\delta$$ and $$\theta$$

I only want solution in the limit of very small $$h_a \rightarrow 0$$ for which except the first $$\alpha$$ value all other values are $$n\pi$$. I have derived an expression to calculate the first value which is:

$$\alpha=\frac 1{\sqrt a} \left( 1+\frac{1}{3 a}-\frac{8}{45 a^2}+\frac{53}{630 a^3}+O\left(\frac{1}{a^4}\right)\right)$$

where $$a=k/(2h_a L)$$. But in any case, I have posted the numerical values in the original question.

Once I get the $$T(x,y,z)$$ my objective is to calculate $$t_h$$ and $$t_c$$ which are given by:

$$t_h=e^{-\beta_h x/L}\bigg(t_{hi} + \frac{\beta_h}{L}\int_0^x e^{\beta_h s/L}T(x,s,w)ds\bigg) \tag I$$

$$t_c=e^{-\beta_c y/l}\bigg(t_{ci} + \frac{\beta_c}{l}\int_0^y e^{\beta_c s/l}T(x,s,0)ds\bigg) \tag J$$

Origins of the b.c.$$E,F$$

Actual bc(s): $$\frac{\partial T(x,y,0)}{\partial z}=p_c (T(x,y,0)-t_c) \tag K$$ $$\frac{\partial T(x,y,w)}{\partial z}=p_h (t_h-T(x,y,w))\tag L$$

where $$t_h,t_c$$ are defined in the equation:

$$\frac{\partial t_c}{\partial y}+\frac{\beta_c}{l}(t_c-T(x,y,0))=0 \tag M$$ $$\frac{\partial t_h}{\partial x}+\frac{\beta_h}{L}(t_h-T(x,y,0))=0 \tag N$$

It is known that $$t_h(x=0)=t_{hi}$$ and $$t_c(y=0)=t_{ci}$$. I had solved $$M,N$$ using the method of integrating factors and used the given conditions to reach $$I,J$$ which were then substituted into the original b.c.(s) $$K,L$$ to reach $$E,F$$.

My attempt I have written the following script to carry out the summation:

γ[α_, δ_] = Sqrt[(α/L)^2 + (δ/l)^2];
L = 0.9; l = 1.8; w = 0.0003; NTUH = 17.394; NTUC = 22.151; ph = 8.6; pc = 13.93;
α0 = 0.01095439637; δ0 = 0.0154917784; β0 = 1.56532; θ0 = 1.56305;
thi = 460; tci = 300; Ta = 380;
V0 = ((c1[α0, β0, δ0, θ0, γ[α0, δ0]] *E^(γ[α0, δ0] *z) + c2[α0, β0, δ0, θ0, γ[α0, δ0]]* E^(-γ[α0, δ0] *z))*Sin[δ0*y/l + θ0] + Sum[(c1[α0, β0, m*\[Pi], 1.5708,γ[α0, m*\[Pi]]] *E^(γ[α0, m*\[Pi]] *z) + c2[α0, β0, m*\[Pi], 1.5708, γ[α0, m*\[Pi]]]*E^(-γ[α0, m*\[Pi]]* z))*Sin[m*\[Pi]*y/l + 1.5708], {m, 1, 5}])*Sin[α0*x/L + β0];
Vn = Sum[((c1[n*\[Pi], 1.5708, δ0, θ0, γ[n*\[Pi], δ0]] *E^(γ[n*\[Pi], δ0] *z) + c2[n*\[Pi], 1.5708, δ0, θ0, γ[n*\[Pi], δ0]]* E^(-γ[n*\[Pi], δ0]* z))*Sin[δ0*y/l + θ0] + Sum[(c1[n*\[Pi], 1.5708, m*\[Pi], 1.5708, γ[n*\[Pi], m*\[Pi]]] *E^(γ[n*\[Pi], m*\[Pi]] *z) + c2[n*\[Pi], 1.5708, m*\[Pi], 1.5708, γ[n*\[Pi], m*\[Pi]]]* E^(-γ[n*\[Pi], m*\[Pi]]* z))*Sin[m*\[Pi]*y/l + 1.5708], {m, 1, 5}])*Sin[n*\[Pi]*x/L +1.5708], {n, 1, 5}];
Vnet = V0 + Vn + Ta;
tcf[x_, y_] = E^(-NTUC* y/l)*{tci + (NTUC/l)*Integrate[E^(NTUC* s/l)*(Vnet /. {y -> s, z -> 0}), {s, 0, y}]};
thf[x_, y_] = E^(-NTUH* x/L)*{thi + (NTUH/L)*Integrate[E^(NTUH* s/L)*(Vnet /. {x -> s, z -> w}), {s, 0, x}]};
tcfavg = Integrate[tcf[x, y], {x, 0, L}]/L;
thfavg = Integrate[thf[x, y], {y, 0, l}]/l;
tcfavg /. y -> l // Chop
thfavg /. x -> L // Chop


The tcfavg and thfavg plots i get are also weird And the outlet temperatures are

tcfavg /. y -> l // Chop
401.984
thfavg /. x -> L // Chop
344.348

• Your T function is on the left hand side dependent on {x,y,z} but on the right hand side not a y. – Steffen Jaeschke Jun 8 '20 at 17:22
• @user2432923 Sorry, that was a typo. Corrected now. – Indrasis Mitra Jun 8 '20 at 17:55
• @AlexTrounev I can understand the frustration on seeing such code. But actually inside mathematica I used $c_1$, but whenever I copy and paste it here on MMA SE, it changes to these expanded notations. Even Greek letters would become \beta. Apologies for this inconvenience. But I have observed that when this code is copied from here to MMA, it is interpreted correctly – Indrasis Mitra Jun 8 '20 at 19:04
• @IndrasisMitra It looks like you try to solve 3D Laplace equation with Dirichlet boundary condition? – Alex Trounev Jun 8 '20 at 19:24
• @AlexTrounev It is the 3D Laplace equation but with similar Robin boundaries on $z$ coordinate and a set of similar Robin conditions on the $x,y$ direction. If it helps should I post the full set of equations that describe the problem ? – Indrasis Mitra Jun 8 '20 at 19:27

We need some numerical model for comparison, so this is one of them based on FEM. First we make sufficient mesh for this problem:

Needs["NDSolveFEM"];Needs["MeshTools"];
L = .90; l = 1.80; w = 0.0003; bh = 17.394;
bc = 22.151; ph = 8.6;
pc = 13.93; pa = 10; n = 10;
thi = 460; tci = 300; Ta = 380; region = Rectangle[{0, 0}, {L, l}];
mesh2D = ToElementMesh[region, MaxCellMeasure -> 5 10^-3 ,
"MeshOrder" -> 1];
mesh3D = ExtrudeMesh[mesh2D, w, 5];
mesh = HexToTetrahedronMesh[mesh3D];

mesh["Wireframe"]


Now we solve the problem by iteration. I have optimized this code, thus it takes about 5 sec:

TC[x_, y_] := tci; TH[x_, y_] := thi;
Do[U[i] =
NDSolveValue[{-Laplacian[u[x, y, z], {x, y, z}] ==
NeumannValue[-pa (u[x, y, z] -
Ta) , (x == 0 || x == L || y == 0 || y == l) & 0 <= z <=
w] + NeumannValue[-pc (u[x, y, z] - TC[x, y]), z == 0] +
NeumannValue[-ph (u[x, y, z] - TH[x, y]), z == w]},
u, {x, y, z} ∈ mesh];
tc[i] = ParametricNDSolveValue[{t'[y] +
bc/l (t[y] - U[i][x, y, 0]) == 0, t[0] == tci},
t, {y, 0, l}, {x}];
th[i] = ParametricNDSolveValue[{t'[x] +
bh/L (t[x] - U[i][x, y, w]) == 0, t[0] == thi},
t, {x, 0, L}, {y}];
TC = Interpolation[
Flatten[Table[{{x, y}, tc[i][x][y]}, {x, 0, L, .02 L}, {y, 0, l,
0.02 l}], 1]];
TH = Interpolation[
Flatten[Table[{{x, y}, th[i][y][x]}, {x, 0, L, .02 L}, {y, 0, l,
0.02 l}], 1]];, {i, 1, n}]


Now we can visualize numerical solution for tc,th in 2 points on every iteration to check how fast solution converges:

Plot[Evaluate[Table[tc[i][L][y], {i, 1, n}]], {y, 0, l},
PlotLegends -> Automatic, AxesLabel -> {"y", "tc(L,y)"}]

Plot[Evaluate[Table[th[i][l][x], {i, 1, n}]], {x, 0, L},
PlotLegends -> Automatic, PlotRange -> All,
AxesLabel -> {"x", "th(x,l)"}]


We see that solution converges fast in 10 steps. Now we can visualize T in 3 slice on z and tc, th on the last iteration

{DensityPlot[U[n][x, y, 0], {x, 0, L}, {y, 0, l},
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
PlotRange -> All],
DensityPlot[U[n][x, y, w/2], {x, 0, L}, {y, 0, l},
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
PlotRange -> All],
DensityPlot[U[n][x, y, w], {x, 0, L}, {y, 0, l},
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
PlotRange -> All]}

{DensityPlot[TC[x, y], {x, 0, L}, {y, 0, l},
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
PlotRange -> All, FrameLabel -> Automatic, PlotLabel -> "tc"],
DensityPlot[TH[x, y], {x, 0, L}, {y, 0, l},
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
PlotRange -> All, FrameLabel -> Automatic, PlotLabel -> "th"]}


Finally we calculate average temperature

tcoldAv = NIntegrate[TC[x, l], {x, 0, L}]/L

Out[]= 381.931

thotAv = NIntegrate[TH[L, y], {y, 0, l}]/l

Out[]= 377.481


Now we can try to improve code for analytical solution. First part of code I just take as it is, but delete two lines and extend number of parameters of functions c1,c2 :

T[x_, y_,
z_] = (C1*E^(\[Gamma] z) + C2 E^(-\[Gamma] z))*
Sin[(\[Alpha] x/L) + \[Beta]]*Sin[(\[Delta] y/l) + \[Theta]] + Ta
tc[x_, y_] =
E^(-NTUC*y/l)*{tci + (NTUC/l)*
Integrate[E^(NTUC*s/l)*T[x, s, 0], {s, 0, y}]};
(*tc[x_,y_]=tc[x,y][[1]];*)
bc1 = (D[T[x, y, z], z] /. z -> 0) == pc (T[x, y, 0] - tc[x, y]);
ortheq1 =
Integrate[(bc1[[1]] - bc1[[2]])*Sin[(\[Alpha] x/L) + \[Beta]]*
Sin[(\[Delta] y/l) + \[Theta]], {x, 0, L}, {y, 0, l},
Assumptions -> {C1 > 0, C2 > 0, L > 0,
l > 0, \[Alpha] > 0, \[Beta] > 0, \[Gamma] > 0, \[Delta] >
0, \[Theta] > 0, NTUC > 0, pc > 0, Ta > 0, tci > 0}] == 0;
(*ortheq1=ortheq1//Simplify;*)
th[x_, y_] =
E^(-NTUH*x/L)*{thi + (NTUH/L)*
Integrate[E^(NTUH*s/L)*T[s, y, w], {s, 0, x}]};
(*th[x_,y_]=th[x,y][[1]];*)
bc2 = (D[T[x, y, z], z] /. z -> w) == ph (th[x, y] - T[x, y, w]);
ortheq2 =
Integrate[(bc2[[1]] - bc2[[2]])*Sin[(\[Alpha] x/L) + \[Beta]]*
Sin[(\[Delta] y/l) + \[Theta]], {x, 0, L}, {y, 0, l},
Assumptions -> {C1 > 0, C2 > 0, L > 0,
l > 0, \[Alpha] > 0, \[Beta] > 0, \[Gamma] > 0, \[Delta] >
0, \[Theta] > 0, NTUC > 0, pc > 0, Ta > 0, thi > 0}] == 0;
(*ortheq2=ortheq2//Simplify;*)
soln = Solve[{ortheq1, ortheq2}, {C1, C2}];
CC1 = C1 /. soln[[1, 1]];
CC2 = C2 /. soln[[1, 2]];
expression1 := CC1;
c1[α_, β_, δ_, θ_, γ_, L_, l_, NTUC_, pc_, Ta_, tci_, NTUH_, ph_, thi_, w_] := Evaluate[expression1];
expression2 := CC2;
c2[α_, β_, δ_, θ_, γ_, L_, l_, NTUC_, pc_, Ta_, tci_, NTUH_, ph_, thi_, w_] := Evaluate[expression2];


Now we run the very fast code for numerical solution

 \[Gamma]1[\[Alpha]_, \[Delta]_] :=
Sqrt[(\[Alpha]/L)^2 + (\[Delta]/l)^2]; m0 = 30; n0 = 30;
L = 0.9; l = 1.8; w = 0.0003; NTUH = 17.394; NTUC = 22.151; ph = 8.6; \
pc = 13.93;
\[Alpha]0 = 0.01095439637; \[Delta]0 = 0.0154917784; \[Beta]0 = \
1.56532; \[Theta]0 = 1.56305;
thi = 460; tci = 300; Ta = 380;
b[n_] := Evaluate[ArcTan[1.66 10^4 (\[Alpha]0 + n Pi)]];
tt[m_] := Evaluate[ArcTan[8.33 10^3 (\[Delta]0 + m*\[Pi])]];
Vn = Sum[(c1[\[Alpha]0 + n*\[Pi], b[n], \[Delta]0 + m*\[Pi],
tt[m], \[Gamma]1[\[Alpha]0 + n*\[Pi], \[Delta]0 + m*\[Pi]], L,
l, pc, pc, Ta, tci, ph, ph, thi, w]*
E^(\[Gamma]1[\[Alpha]0 + n*\[Pi], \[Delta]0 + m*\[Pi]]*z) +
c2[\[Alpha]0 + n*\[Pi], b[n], \[Delta]0 + m*\[Pi],
tt[m], \[Gamma]1[\[Alpha]0 + n*\[Pi], \[Delta]0 + m*\[Pi]], L,
l, pc, pc, Ta, tci, ph, ph, thi, w]*
E^(-\[Gamma]1[\[Alpha]0 + n*\[Pi], \[Delta]0 + m*\[Pi]]*z))*
Sin[(\[Delta]0 + m*\[Pi])*y/l + tt[m]]*
Sin[(\[Alpha]0 + n*\[Pi])*x/L + b[n]], {n, 0, n0}, {m, 0, m0}];
Vnet = Vn/2 + Ta;

tc = ParametricNDSolveValue[{t'[y] + pc/l (t[y] - Vnet /. z -> 0) ==
0, t[0] == tci}, t, {y, 0, l}, {x}]; th =
ParametricNDSolveValue[{t'[x] + ph/L (t[x] - Vnet /. z -> w) == 0,
t[0] == thi}, t, {x, 0, L}, {y}]; TC =
Interpolation[
Flatten[Table[{{x, y}, tc[x][y]}, {x, 0, L, .01 L}, {y, 0, l,
0.01 l}], 1]]; TH =
Interpolation[
Flatten[Table[{{x, y}, th[y][x]}, {x, 0, L, .01 L}, {y, 0, l,
0.01 l}], 1]];


Note, I am using Vn/2 to limit low and high temperature. And finally we visualize solution

{DensityPlot[TC[x, y], {x, 0, L}, {y, 0, l},
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
PlotRange -> All, FrameLabel -> Automatic, PlotLabel -> "tc"],
DensityPlot[TH[x, y], {x, 0, L}, {y, 0, l},
ColorFunction -> "Rainbow", PlotLegends -> Automatic,
PlotRange -> All, FrameLabel -> Automatic, PlotLabel -> "th"]}


• Really appreciate this effort. Thanks. From the time the last answer (by Steffen) was posted on this question, I have managed to write a script to carry out the summation $(1)$.Should I post the code here if you would like to take a look ? – Indrasis Mitra Jun 13 '20 at 10:57
• I have posted my code to sum $(1)$ as an edit to the original question. I do not know why but the average temperatures I get are far off. If you are interested further I will post the equations to the V0,Vn terms in my code. NOTE: The term bh,bc in your code is NTUH,NTUC in mine.Also the functions c1,c2 defined in the attempt come from the first code section which I have edited now to reflect these definitions. – Indrasis Mitra Jun 13 '20 at 12:21
• @IndrasisMitra Could you please post parameter $h_a$? I put it as pa=10. – Alex Trounev Jun 13 '20 at 14:05
• @IndrasisMitra Code been speed up due to replacement of symbolic Integrate with numerical ParametricNDSolveValue. – Alex Trounev Jun 19 '20 at 11:42
• @IndrasisMitra This normalization is due to definition of c1,c2. Please, check how you get this coefficients. – Alex Trounev Jun 19 '20 at 15:58

Your T function is on the left-hand side dependent on {x,y,z} but on the right-hand side not a y in the MathML code. You got confused by the name of the functions in special states of the solutions process and forget to use them consequently. The solution of the Subscript[C,1], Subscript[C,2] depends in length on the given parameters but all are not set in the definitions above. It is a deviation from the solution path not to name the solution special at the end of the first Mathematica code section.

T[x_, y_, z_] = (Subscript[C, 1] E^(γ z) + Subscript[C, 2] E^(- γ z))*Sin[(α x/L) + β]*Sin[(δ y/l) + θ] + Subscript[T, a]
tc[x_, y_] = E^(-Subscript[β, c] y/l)*{tci + (Subscript[β, c]/l)*Integrate[E^(Subscript[β, c] s/l)*T[x, s, 0], {s, 0, y}]};
tc[x_, y_] = tc[x, y][[1]];
bc1 = (D[T[x, y, z], z] /. z -> 0) == Subscript[p, c] (T[x, y, 0] - tc[x, y]);
ortheq1 = Integrate[bc1[[1]]*Sin[(α x/L) + β]*Sin[(δ y/l) + θ], {x, 0, L}, {y, 0, l}] == Integrate[bc1[[2]]*Sin[(α x/L) + β]*Sin[(δ y/l) + θ], {x, 0, L}, {y, 0, l}];
ortheq1 = ortheq1 // Simplify;
th[x_, y_] = E^(-Subscript[β, h] x/L)*{thi + (Subscript[β, h]/L)*Integrate[E^(Subscript[β, h] s/L)*T[s, y, w], {s, 0, x}]};
th[x_, y_] = th[x, y][[1]];
bc2 = (D[T[x, y, z], z] /. z -> w) == Subscript[p, h] (th[x, y] - T[x, y, w]);
ortheq2 = Integrate[bc2[[1]]*Sin[(α x/L) + β]*Sin[(δ y/l) + θ], {x, 0, L}, {y, 0, l}] == Integrate[bc2[[2]]*Sin[(α x/L) + β]*Sin[(δ y/l) + θ], {x, 0, L}, {y, 0, l}];
ortheq2 = ortheq2 // Simplify;
soln = Solve[{ortheq1, ortheq2}, {Subscript[C, 1], Subscript[C, 2]}];
Subscript[Csol, 1] = Subscript[C, 1] /. soln[[1, 1]];
Subscript[Csol, 2] = Subscript[C, 2] /. soln[[1, 2]];


From that plug into the definition:

Tsol[x_, y_, z_] = (Subscript[Csol, 1] E^(γ z) + Subscript[Csol, 2] E^(- γ z))*Sin[(α x/L) + β]*Sin[(δ y/l) + θ] + Subscript[T, a]


This Tsol is Your Twnet the variables and parameters plugged in correctly.

It is much better to define:

T[x_, y_, z_,γ_,α_,β_,δ_,θ_,L_,l_,Subscript[T_, a]]

so that another source of confusion. Might be a good idea to name such complicated variable parameters as Subscript[T_, a] shorter such as T_.

Doing so the second part of Your Mathematica code does take a long time too.

α = 0.01095; δ = 0.1549;
β = ArcTan[1.66*10^4 α]; θ =
Tan[δ/(10^3 * 8.33)];

TWnet = (Subscript[Csol, 1] E^(γ z) +
Subscript[Csol, 2] E^(-γ z))*
Sin[(α x/L) + β]*Sin[(δ y/l) + θ] +
Subscript[T, a];

L = 0.9; l = 1.8; w = 0.0003; Subscript[β, h] = 17.394;
Subscript[β, c] = 22.151; Subscript[p, h] = 8.6;
Subscript[p, c] = 13.93;
γ = Sqrt[(α/L)^2 + (δ/l)^2];
thi = 460; tci = 300; Subscript[T, a] = 380;
tc1[x_, y_] =
E^(-Subscript[β, c] y/l)*{tci + (Subscript[β, c]/l)*
Integrate[
E^(Subscript[β, c] s/l)*(TWnet /. {y -> s, z -> 0}), {s,
0, y}]};
th1[x_, y_] =
E^(-Subscript[β, h] x/L)*{thi + (Subscript[β, h]/L)*
Integrate[
E^(Subscript[β, h] s/L)*(TWnet /. {x -> s, z -> w}), {s,
0, x}]};
Plot[tc1[x, l], {x, 0, L}]
Plot[th1[L, y], {y, 0, l}]
THotAvg = Integrate[th1[x, y]/l, {y, 0, l}];
TColdAvg = Integrate[tc1[x, y]/L, {x, 0, L}];
THotAvg /. x -> L
TColdAvg /. y -> l
Plot[THotAvg, {x, 0, L}]
Plot[TColdAvg, {y, 0, l}]


{408.044}

{433.444}


This is as close stuck to the given information and independent of n and m.

A start is

nmax = 3; mmax = 3;

T[x_, y_, z_,γ_,α_,β_,δ_,θ_,L_,l_,Subscript[T_, a]] =
Sum[(Subscript[C, 1] E^(γ z) +
Subscript[C, 2] E^(-γ z))*
Sin[(Subscript[α, n] x/L) + Subscript[β, n]]*
Sin[(Subscript[δ, m] y/l) + Subscript[θ, m]] +
Subscript[T, a], {n, 0, nmax}, {m, 0, mmax}]


And solve for each n and m.

• How long this code run? – Alex Trounev Jun 8 '20 at 19:31
• @user2432923 Thanks for this. But as I have said in the original question $\alpha_0, \delta_0$ are the values that I want to supply. From here i.e. n=1` the subsequent values are defined by $\alpha=n\pi$ and $\delta=m\pi$. This has not been probably used in your answer or am I misunderstanding something ? – Indrasis Mitra Jun 8 '20 at 20:23
• @AlexTrounev: I did not time that, but it was about an hour, just the first part of the code. The second worked fast after inserting the values for the coefficients. – Steffen Jaeschke Jun 9 '20 at 8:29
• @Indrasis Mitra: I am the opinion that this will take me a too long time. Such Fourier Expansions are usually convergent very rapidly, so it suffices to make a plot of the amplitude functions c1 and c2 that are not dependent on the expansion parameter n and m or Your T function have to be corrected for that dependencies. – Steffen Jaeschke Jun 9 '20 at 8:38