5
$\begingroup$

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 tcfavg thfavg And the outlet temperatures are

tcfavg /. y -> l // Chop
401.984
thfavg /. x -> L // Chop
344.348
$\endgroup$
10
  • $\begingroup$ Your T function is on the left hand side dependent on {x,y,z} but on the right hand side not a y. $\endgroup$ Jun 8, 2020 at 17:22
  • $\begingroup$ @user2432923 Sorry, that was a typo. Corrected now. $\endgroup$
    – Avrana
    Jun 8, 2020 at 17:55
  • $\begingroup$ @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 $\endgroup$
    – Avrana
    Jun 8, 2020 at 19:04
  • $\begingroup$ @IndrasisMitra It looks like you try to solve 3D Laplace equation with Dirichlet boundary condition? $\endgroup$ Jun 8, 2020 at 19:24
  • $\begingroup$ @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 ? $\endgroup$
    – Avrana
    Jun 8, 2020 at 19:27

2 Answers 2

7
+50
$\begingroup$

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["NDSolve`FEM`"];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)"}] 

Figure 1 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"]}

Figure 2

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"]}

Figure 2

$\endgroup$
15
  • $\begingroup$ 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 ? $\endgroup$
    – Avrana
    Jun 13, 2020 at 10:57
  • $\begingroup$ 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. $\endgroup$
    – Avrana
    Jun 13, 2020 at 12:21
  • $\begingroup$ @IndrasisMitra Could you please post parameter $h_a$? I put it as pa=10. $\endgroup$ Jun 13, 2020 at 14:05
  • 1
    $\begingroup$ @IndrasisMitra Code been speed up due to replacement of symbolic Integrate with numerical ParametricNDSolveValue. $\endgroup$ Jun 19, 2020 at 11:42
  • 1
    $\begingroup$ @IndrasisMitra This normalization is due to definition of c1,c2. Please, check how you get this coefficients. $\endgroup$ Jun 19, 2020 at 15:58
5
$\begingroup$

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}]

Graphics

Graphics

{408.044}

{433.444}

Graphics

Graphics

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.

$\endgroup$
4
  • 1
    $\begingroup$ How long this code run? $\endgroup$ Jun 8, 2020 at 19:31
  • 2
    $\begingroup$ @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 ? $\endgroup$
    – Avrana
    Jun 8, 2020 at 20:23
  • $\begingroup$ @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. $\endgroup$ Jun 9, 2020 at 8:29
  • $\begingroup$ @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. $\endgroup$ Jun 9, 2020 at 8:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.