# 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][];
bc1 = (D[T[x, y, z], z] /. z -> 0) == pc (T[x, y, 0] - tc[x, y]);
ortheq1 = Integrate[bc1[]*Sin[(α x/L) + β]*Sin[(δ y/l) + θ], {x, 0, L}, {y, 0, l}] == Integrate[bc1[]*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][];
bc2 = (D[T[x, y, z], z] /. z -> w) == ph (th[x, y] - T[x, y, w]);
ortheq2 = Integrate[bc2[]*Sin[(α x/L) + β]*Sin[(δ y/l) + θ], {x, 0, L}, {y, 0, l}] == Integrate[bc2[]*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. Jun 8, 2020 at 17:22
• @user2432923 Sorry, that was a typo. Corrected now. Jun 8, 2020 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 Jun 8, 2020 at 19:04
• @IndrasisMitra It looks like you try to solve 3D Laplace equation with Dirichlet boundary condition? Jun 8, 2020 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 ? Jun 8, 2020 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 == tci},
t, {y, 0, l}, {x}];
th[i] = ParametricNDSolveValue[{t'[x] +
bh/L (t[x] - U[i][x, y, w]) == 0, t == 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][];*)
bc1 = (D[T[x, y, z], z] /. z -> 0) == pc (T[x, y, 0] - tc[x, y]);
ortheq1 =
Integrate[(bc1[] - bc1[])*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][];*)
bc2 = (D[T[x, y, z], z] /. z -> w) == ph (th[x, y] - T[x, y, w]);
ortheq2 =
Integrate[(bc2[] - bc2[])*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 == tci}, t, {y, 0, l}, {x}]; th =
ParametricNDSolveValue[{t'[x] + ph/L (t[x] - Vnet /. z -> w) == 0,
t == 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 ? Jun 13, 2020 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. Jun 13, 2020 at 12:21
• @IndrasisMitra Could you please post parameter $h_a$? I put it as pa=10. Jun 13, 2020 at 14:05
• @IndrasisMitra Code been speed up due to replacement of symbolic Integrate with numerical ParametricNDSolveValue. Jun 19, 2020 at 11:42
• @IndrasisMitra This normalization is due to definition of c1,c2. Please, check how you get this coefficients. Jun 19, 2020 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][];
bc1 = (D[T[x, y, z], z] /. z -> 0) == Subscript[p, c] (T[x, y, 0] - tc[x, y]);
ortheq1 = Integrate[bc1[]*Sin[(α x/L) + β]*Sin[(δ y/l) + θ], {x, 0, L}, {y, 0, l}] == Integrate[bc1[]*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][];
bc2 = (D[T[x, y, z], z] /. z -> w) == Subscript[p, h] (th[x, y] - T[x, y, w]);
ortheq2 = Integrate[bc2[]*Sin[(α x/L) + β]*Sin[(δ y/l) + θ], {x, 0, L}, {y, 0, l}] == Integrate[bc2[]*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? Jun 8, 2020 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 ? Jun 8, 2020 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. Jun 9, 2020 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. Jun 9, 2020 at 8:38