# Coupled heat transfer equations using collocation method

Using the collocation method proposed here, recently this problem has been solved. I am trying to solve a similar problem described by the equations given below. My attempt in Mathematica is following.

$$T$$ is in the domain $$x\in[0,1], y\in[0,1]$$

eq1 = \[Lambda]x D[T[x, y], x, x] + \[Lambda]y D[T[x, y], y, y] ==
0; bc1 = {(D[T[x, y], y] + rh (Tfh[x] - T[x, y]) == 0) /.
y -> 1, (D[T[x, y], y] + rc (T[x, y] - Tfc[x]) == 0) /.
y -> 0}; eq2 = D[T[x], x] + bh (Tfh[x] - T[x]) == 0;
bc2 = Tfh == 1;
eq3 = D[Tfc[x], x] + bc (T[x] - Tfc[x]) == 0;
bc3 = Tfc == 0;

UE[m_, t_] := Cos[m t] Exp[-m t]
nn = 3;
dx = 1/(nn); xl = Table[l*dx, {l, 0, nn}]; ycol =
xcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, nn + 1}]; Psijk =
Table[UE[n, t1], {n, 0, nn - 1}]; Int1 = Integrate[Psijk, t1];
Int2 = Integrate[Int1, t1];
Psi[y_] := Psijk /. t1 -> y;
int1[y_] := Int1 /. t1 -> y;
int2[y_] := Int2 /. t1 -> y; M = nn;
M = nn; U1 = Array[a1, {M, M}]; U2 = Array[a2, {M, M}]; G1 =
Array[g1, {M}]; G2 = Array[g2, {M}]; G4 = Array[g4, {M}]; G5 =
Array[g5, {M}];

Tfhx[x_] := (Psi[x].G5);
Tfcy[x_] := (Psi[x].G4);
Tfh[x_] := (int1[x].G5);
Tfc[x_] := (int2[x].G4);

u1[x_, y_] := (int2[x].U1.Psi[y]) + x Psi[y].G1;
u2[x_, y_] := (Psi[x].U2.int2[y]) + y Psi[x].G2;
uy[x_, y_] := (Psi[x].U2.int1[y]) + Psi[x].G2;
ux[x_, y_] := (int1[x].U1.Psi[y]) + Psi[y].G1;
uxx[x_, y_] := (Psi[x].U1.Psi[y]);
uyy[x_, y_] := (Psi[x].U2.Psi[y]);

\[Lambda]x = 1/0.025^2; \[Lambda]y = 1/0.002^2; bh = 0.625; bc =
4 bh; rh = 3000 0.002/390; rc = 3000 0.002/390;

eqn = Join[
Flatten[Table[(\[Lambda]x uxx[xcol[[i]],
ycol[[j]]] + \[Lambda]y uyy[xcol[[i]], ycol[[j]]]) == 0, {i,
M}, {j, M}]],
Flatten[Table[
u1[xcol[[i]], ycol[[j]]] - u2[xcol[[i]], ycol[[j]]] == 0, {i,
M}, {j, M}]], Flatten[Table[ux[1, ycol[[i]]] == 0, {i, M}]],
Flatten[Table[uy[xcol[[i]], 1] == 0, {i, M}]],
Flatten[Table[ux[0, ycol[[i]]] == 0, {i, M}]],
Flatten[Table[uy[xcol[[i]], 0] == 0, {i, M}]],
Flatten[Table[
uy[xcol[[i]], 1] + rh (Tfh[xcol[[i]]] - u2[xcol[[i]], 1]) ==
0, {i, M}]],
Flatten[Table[
uy[xcol[[i]], 0] + rc (u2[xcol[[i]], 0] - Tfc[xcol[[i]]]) ==
0, {i, M}]],
Flatten[Table[
Tfhx[xcol[[i]]] + bh (Tfh[xcol[[i]]] - u2[xcol[[i]], 1]) ==
0, {i, M}]],
Flatten[Table[
tcy[xcol[[i]]] + bc (u2[xcol[[i]], 0] - Tfc[xcol[[i]]]) == 0, {i,
M}]], Table[Tfh[xcol[[i]]] == 1., {i, M}],
Table[Tfc[xcol[[i]], 0] == 0., {i, M}]];
var = Join[Flatten[U1], Flatten[U2], Flatten[G1], Flatten[G2],
Flatten[G4], Flatten[G5]];

{v, mat} = CoefficientArrays[eqn, var];

sol = LinearSolve[mat, -v];

rule = Table[var[[i]] -> sol[[i]], {i, Length[var]}];

{Plot[Evaluate[Tfc[x] /. rule], {x, 0, 1}, ColorFunction -> "Rainbow",
MeshStyle -> White, PlotLegends -> Automatic, PlotLabel -> Tfc],
Plot3D[Evaluate[Tfh[x] /. rule], {x, 0, 1},
ColorFunction -> "Rainbow", MeshStyle -> White,
PlotLegends -> Automatic, PlotLabel -> Tfh],
Table[Plot3D[Evaluate[u1[x, y] /. rule], {x, 0, 1}, {y, 0, 1},
ColorFunction -> "Rainbow", MeshStyle -> White,
PlotLegends -> Automatic, PlotLabel -> T[y]], {y, 0, 1, .5}]}


The code above tries to solve the following equations:

$$\lambda_x \frac{\partial^2 T}{\partial x^2}+\lambda_y \frac{\partial^2 T}{\partial y^2}=0 \tag1$$

Zero gradient of $$T$$ at $$x=0,1$$.

At $$y=0,1$$ the solid ($$T$$) is exposed to two different fluids following opposite to each other (separated by the solid).

$$\frac{\partial T(x,0)}{\partial y}+r_c(T(x,0)-t_c)=0\tag2$$ $$\frac{\partial T(x,1)}{\partial y}+r_h(t_h-T(x,1))=0\tag3$$

The fluid temperatures $$t_h,t_c$$ are governed by:

$$\frac{\partial t_h}{\partial x}+\beta_h(t_h-T)=0\tag4$$ $$\frac{\partial t_c}{\partial x}+\beta_c(T-t_c)=0\tag5$$

A reduced order model

The above problem can be simplified if one averages the solid temperature along the $$y$$ direction to get the following solid governing equation instead of $$(1)$$:

$$\kappa \frac{\mathrm{d}^2 T}{\mathrm{d} x^2} + \mu b_h(t_h-T) - \nu b_c(T-t_c)=0 \tag6$$

with $$T'(0)=T'(1)=0$$. Equation $$(6)$$ is still coupled with Equation $$(4)$$ and Equation $$(5)$$. Some parameter values are bc=12.38, bh=25.32, mu=1.143, nu=1, kappa=2.16

Some cases behaving wierdly

Although the wavelet method works really fine, but I am facing problems with some particular set of flow configurations, which are physically important. One of them being the following:

λx = 1/0.025^2; λy = 1/0.001^2; bh = 173.6539; bc = 355.1724; rh = 134.31 0.001/16; rc = 305.2252 0.001/16;

Upon executing these parameters, with the wavelet method such that nn=16, I obtain the following plots of the fluid temperature profiles For nn=40, the plot is like the following: Lastly, the nn=96 gives: As it is clear that the fluid temperature profiles experience spikes. Although, they do decrease in magnitude with increasing nn. I would like to know the reason as to why this is happening and if they can be removed. I have checked that the insulation boundary conditions on x=0,1 are satisfied to the order of 10^-12.

• There is a typo in eq2 = D[T[x], x] + bh (Tfh[x] - T[x]) == 0;. Feb 13, 2022 at 3:00
• If fluid flows at y=0, y=1are opposite then it should be eq2 = D[Tfh[x], x] + bh (Tfh[x] - T[x, 1]) == 0; bc2 = Tfh == 1; eq3 = D[Tfc[x], x] + bc (T[x, 0] - Tfc[x]) == 0; bc3 = Tfc == 0; Feb 13, 2022 at 3:32
• @AlexTrounev yes you are right about the typos in eq2 and eq3, but they are not being used while solving the system. The implementation of these in eqn, I have used bh (Tfh[xcol[[i]]] - u2[xcol[[i]], 1]) and similarly for the bc equation. Feb 13, 2022 at 4:04
• There is also typo in eqn in the line tcy[xcol[[i]]] + bc (u2[xcol[[i]], 0] - Tfc[xcol[[i]]]) == 0. tcy is not defined. Feb 13, 2022 at 4:53
• In the last example there is numerical instability due to large bc, bh. To avoid this instability we need take dx<Min[{1/bc,1/bh}], and nn=Round[1/dx ] consequently. Therefore, nn>355 for this case. So, we need supercomputer or other method to compute cases with large bc, bh. Feb 15, 2022 at 16:28

Solution with Euler wavelets:

eq1 = \[Lambda]x D[T[x, y], x, x] + \[Lambda]y D[T[x, y], y, y] ==
0; bc1 = {(D[T[x, y], y] + rh (Tfh[x] - T[x, y]) == 0) /.
y -> 1, (D[T[x, y], y] + rc (T[x, y] - Tfc[x]) == 0) /.
y -> 0}; eq2 = D[Tfh[x], x] + bh (Tfh[x] - T[x, 1]) == 0;
bc2 = Tfh == 1;
eq3 = D[Tfc[x], x] + bc (T[x, 0] - Tfc[x]) == 0;
bc3 = Tfc == 0;

UE[m_, t_] := EulerE[m, t]
psi[k_, n_, m_, t_] :=
Piecewise[{{2^(k/2) Sqrt[2/Pi] UE[m, 2^k t - 2 n + 1], (n - 1)/
2^(k - 1) <= t < n/2^(k - 1)}, {0, True}}]
PsiE[k_, M_, t_] :=
Flatten[Table[psi[k, n, m, t], {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]
k0 = 3; M0 = 4; With[{k = k0, M = M0},
var1 = Flatten[Table[c[n, m], {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]];
nn = Length[var1]
dx = 1/(nn);  xl = Table[ l*dx, {l, 0, nn}]; ycol =
xcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, nn + 1}]; Psijk =
With[{k = k0, M = M0}, PsiE[k, M, t1]]; Int1 =
With[{k = k0, M = M0}, Integrate[PsiE[k, M, t1], t1]];
Int2 = Integrate[Int1, t1];
Psi[y_] := Psijk /. t1 -> y; int1[y_] := Int1 /. t1 -> y;
int2[y_] := Int2 /. t1 -> y;
M = nn; U1 = Array[a1, {M, M}]; U2 = Array[a2, {M, M}]; G1 =
Array[g1, {M}]; G2 = Array[g2, {M}]; G3 = Array[g3, {M}]; G4 =
Array[g4, {M}]; F1 = Array[f1, {M}]; F2 = Array[f2, {M}];

u1[x_, y_] := int2[x] . U1 . Psi[y] + x G1 . Psi[y] + F1 . Psi[y];
u2[x_, y_] := Psi[x] . U2 . int2[y] + y G2 . Psi[x] + F2 . Psi[x];
uy[x_, y_] := Psi[x] . U2 . int1[y] + G2 . Psi[x];
ux[x_, y_] := int1[x] . U1 . Psi[y] + G1 . Psi[y];
uxx[x_, y_] := Psi[x] . U1 . Psi[y];
uyy[x_, y_] := Psi[x] . U2 . Psi[y];
Tfhx[x_] := Psi[x] . G3; Tfcx[x_] := Psi[x] . G4;
Tfh[x_] := int1[x] . G3 + th0; Tfc[x_] := int1[x] . G4 + tc0;

\[Lambda]x = 1/0.025^2; \[Lambda]y = 1/0.002^2; bh = 0.625; bc =
4 bh; rh = 3000 0.002/390; rc = 3000 0.002/390;

eqn = Join[
Flatten[Table[(\[Lambda]x uxx[xcol[[i]],
ycol[[j]]] + \[Lambda]y uyy[xcol[[i]], ycol[[j]]]) == 0, {i,
M}, {j, M}]],
Flatten[Table[
u1[xcol[[i]], ycol[[j]]] - u2[xcol[[i]], ycol[[j]]] == 0, {i,
M}, {j, M}]], Table[ux[1., ycol[[i]]] == 0, {i, M}],
Table[ux[0., ycol[[i]]] == 0, {i, M}],
Table[uy[xcol[[i]], 1.] + rh (Tfh[xcol[[i]]] - u2[xcol[[i]], 1.]) ==
0, {i, M}],
Table[uy[xcol[[i]], 0.] + rc (u2[xcol[[i]], 0.] - Tfc[xcol[[i]]]) ==
0, {i, M}],
Table[Tfhx[xcol[[i]]] + bh (Tfh[xcol[[i]]] - u2[xcol[[i]], 1.]) ==
0, {i, M}],
Table[Tfcx[xcol[[i]]] + bc (u2[xcol[[i]], 0] - Tfc[xcol[[i]]]) ==
0, {i, M}], {Tfh[0.] == 1.}, {Tfc[1.] == 0.}];
var = Join[Flatten[U1], Flatten[U2], G1, G2, G3, G4, F1,
F2, {th0, tc0}];

{v, mat} = CoefficientArrays[eqn, var];

sol = LinearSolve[mat // N, -v];

rule = Table[var[[i]] -> sol[[i]], {i, Length[var]}];


Visualization

{Plot[Evaluate[Tfc[x] /. rule], {x, 0, 1}, PlotLabel -> Tfc],
Plot[Evaluate[Tfh[x] /. rule], {x, 0, 1}, PlotLabel -> Tfh],
Plot3D[Evaluate[u1[x, y] /. rule], {x, 0, 1}, {y, 0, 1},
ColorFunction -> Hue, MeshStyle -> Black, PlotLegends -> Automatic,
PlotLabel -> T, PlotTheme -> "Marketing", AxesLabel -> Automatic]} This code can be testified on antisymmetric solution with bc=bh, rc=rh, Tfh==1, Tfc==-1, then we should have solution satisfies equation tfc[x]+tfh[1-x]==0. Example 1:

eq1 = \[Lambda]x D[T[x, y], x, x] + \[Lambda]y D[T[x, y], y, y] ==
0; bc1 = {(D[T[x, y], y] + rh (Tfh[x] - T[x, y]) == 0) /.
y -> 1, (D[T[x, y], y] + rc (T[x, y] - Tfc[x]) == 0) /.
y -> 0}; eq2 = D[Tfh[x], x] + bh (Tfh[x] - T[x, 1]) == 0;
bc2 = Tfh == 1;
eq3 = D[Tfc[x], x] + bc (T[x, 0] - Tfc[x]) == 0;
bc3 = Tfc == -1;

UE[m_, t_] := EulerE[m, t]
psi[k_, n_, m_, t_] :=
Piecewise[{{2^(k/2) Sqrt[2/Pi] UE[m, 2^k t - 2 n + 1], (n - 1)/
2^(k - 1) <= t < n/2^(k - 1)}, {0, True}}]
PsiE[k_, M_, t_] :=
Flatten[Table[psi[k, n, m, t], {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]
k0 = 2; M0 = 7; With[{k = k0, M = M0},
var1 = Flatten[Table[c[n, m], {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]];
nn = Length[var1]
dx = 1/(nn);  xl = Table[ l*dx, {l, 0, nn}]; ycol =
xcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, nn + 1}]; Psijk =
With[{k = k0, M = M0}, PsiE[k, M, t1]]; Int1 =
With[{k = k0, M = M0}, Integrate[PsiE[k, M, t1], t1]];
Int2 = Integrate[Int1, t1];
Psi[y_] := Psijk /. t1 -> y; int1[y_] := Int1 /. t1 -> y;
int2[y_] := Int2 /. t1 -> y;
M = nn; U1 = Array[a1, {M, M}]; U2 = Array[a2, {M, M}]; G1 =
Array[g1, {M}]; G2 = Array[g2, {M}]; G3 = Array[g3, {M}]; G4 =
Array[g4, {M}]; F1 = Array[f1, {M}]; F2 = Array[f2, {M}];

u1[x_, y_] := int2[x] . U1 . Psi[y] + x G1 . Psi[y] + F1 . Psi[y];
u2[x_, y_] := Psi[x] . U2 . int2[y] + y G2 . Psi[x] + F2 . Psi[x];
uy[x_, y_] := Psi[x] . U2 . int1[y] + G2 . Psi[x];
ux[x_, y_] := int1[x] . U1 . Psi[y] + G1 . Psi[y];
uxx[x_, y_] := Psi[x] . U1 . Psi[y];
uyy[x_, y_] := Psi[x] . U2 . Psi[y];
Tfhx[x_] := Psi[x] . G3; Tfcx[x_] := Psi[x] . G4;
Tfh[x_] := int1[x] . G3 + th0; Tfc[x_] := int1[x] . G4 + tc0;

\[Lambda]x = 1/0.025^2; \[Lambda]y =
1/0.002^2; bh = bc = 0.625; rh = 3000 0.002/390; rc = 3000 0.002/390;

eqn = Join[
Flatten[Table[(\[Lambda]x uxx[xcol[[i]],
ycol[[j]]] + \[Lambda]y uyy[xcol[[i]], ycol[[j]]]) == 0, {i,
M}, {j, M}]],
Flatten[Table[
u1[xcol[[i]], ycol[[j]]] - u2[xcol[[i]], ycol[[j]]] == 0, {i,
M}, {j, M}]], Table[ux[1., ycol[[i]]] == 0, {i, M}],
Table[ux[0., ycol[[i]]] == 0, {i, M}],
Table[uy[xcol[[i]], 1.] + rh (Tfh[xcol[[i]]] - u2[xcol[[i]], 1.]) ==
0, {i, M}],
Table[uy[xcol[[i]], 0.] + rc (u2[xcol[[i]], 0.] - Tfc[xcol[[i]]]) ==
0, {i, M}],
Table[Tfhx[xcol[[i]]] + bh (Tfh[xcol[[i]]] - u2[xcol[[i]], 1.]) ==
0, {i, M}],
Table[Tfcx[xcol[[i]]] + bc (u2[xcol[[i]], 0] - Tfc[xcol[[i]]]) ==
0, {i, M}], {Tfh[0.] == 1.}, {Tfc[1.] == -1.}];
var = Join[Flatten[U1], Flatten[U2], G1, G2, G3, G4, F1,
F2, {th0, tc0}];

{v, mat} = CoefficientArrays[eqn, var];

sol = LinearSolve[mat // N, -v];

rule = Table[var[[i]] -> sol[[i]], {i, Length[var]}];


Visualization

{Plot[Evaluate[Tfc[x] /. rule], {x, 0, 1}, PlotLabel -> Tfc],
Plot[Evaluate[Tfh[x] /. rule], {x, 0, 1}, PlotLabel -> Tfh,
PlotRange -> All],
Plot[Evaluate[Tfc[x] + Tfh[1 - x]] /. rule, {x, 0, 1},
PlotLabel -> "Tfc[x]+Tfh[1-x]"],
Plot3D[Evaluate[u1[x, y] /. rule], {x, 0, 1}, {y, 0, 1},
ColorFunction -> Hue, MeshStyle -> Black, PlotLabel -> T,
PlotTheme -> "Marketing", AxesLabel -> Automatic]} Therefore the max error of sum Tfc[x] + Tfh[1 - x]] is about $$6\times 10^{-14}$$ for M=14.

• Thankyou for the answer. However, I have a doubt. Can you inform where the boundary conditions of zero gradient of $T$ at $x=0,1$ or $\frac{\partial T(0,y)}{\partial x}=\frac{\partial T(1,y)}{\partial x}=0$are being implemented in this methodology. Do the Euler wavelets implicitly assume so? Feb 13, 2022 at 13:40
• I am doubtful about this because when I run this system for bh=bc, which corresponds to equal flow rates of both the hot and cold fluid, the temperature rise in the cold fluid is not equal to the temperature drop in the hot fluid. However, in the case of equal flow they should be same. In other words for equal flow or bc=bh and when the solid faces at x=0 and x=1 are insulated, one is supposed to get Tfc+Tfh=1. Feb 13, 2022 at 13:44
• However, when I use the UE[m_, t_] := Cos[m t] Exp[-m t], the problem of energy balance is taken care of. I think the cos eigen-functions ensure zero $T$ gradients at $x=0,1$. For example for bh=bc=0.625, Tfc+Tfh=1. with this function as the basis. Feb 13, 2022 at 13:58
• I have also posted the governing equations for a reduced-order model for this problem. The solution to this model is partially available in literature for diffferent set of b.c.(s), i.e. $T'(0)=T'(1)=0$ and $T'(0)=1,T(1)=0$. Please give it a try as it might be helpful for validation of the present system. Feb 13, 2022 at 14:19
• Also, although I have been able to achieve some success with UE[m_, t_] := Cos[m t] Exp[-m t], the temperature profiles I recieve with this set of parameter values \[Lambda]x = 1/0.025^2; \[Lambda]y = 1/0.002^2; bh = 25.32; bc = 12.38; rh = 305.2252 0.002/16; rc = 134.31 0.002/16;` are very unphysical. The fluid temperatures basically drop and then rise again, which should not happen. Have a look if time permits. Feb 13, 2022 at 14:31