# Finite element for coupled 1D and 3D equations

This is a renewed post for this question, since the old post is already 2 years old and the problem hasn't been solved completely.

The problem is to solve heat conduction with robin type boundary condition, which is coupled with a 1D equation. The 1D equation describes the heat convection between fluid flow in a pipe and the surface of the pipe. This equation is necessary, because I would like to determine not only the temperature distribution inside the solid domain, but also the fluid flow temperature distribution along the z-axis.

For simplification, let's ignore the pipe wall and only consider steady-state first. The domain is a cuboid with a half cylinder hole. The equations are as follows:

Eq1: $$\Delta T(x,y,z) = 0$$

BC1: $$k \frac{\partial T(x,y,z)}{\partial n} = \alpha_a(T(x,y,z) - T_e)$$ if $$y = 0$$

BC2: $$- k \frac{\partial T(x,y,z)}{\partial n} = \alpha_i(T(x,y,z) - T_e)$$ if $$y = H$$

BC3: $$\frac{\partial T(x,y,z)}{\partial n} = 0$$ if $$x = 0$$

BC4: $$\frac{\partial T(x,y,z)}{\partial n} = 0$$ if $$x = \frac{L_p}{2}$$

BC5: $$\frac{\partial T(x,y,z)}{\partial n} = 0$$ if $$z = 0$$

BC6: $$\frac{\partial T(x,y,z)}{\partial n} = 0$$ if $$z = l$$

BC7: $$- k \frac{\partial T(x,y,z)}{\partial n} = \alpha_f (T(x,y,z) - T_f(z))$$ if $$x^2 + (y - H_a )^2 = r_a^2$$

Eq2: $$- c_W \rho_W \dot{V}_f \frac{dT_f(z)}{dz} = 2 \alpha_f \pi r_a (T_f(z) - T_m(z))$$

where $$T_m(z)$$ should be the average temperature at position z in the internal curve,

$$T_m(z) = \frac{1}{\pi}\int_{-\pi/2}^{\pi/2} T[r_a \cos(\vartheta), H_a + r_a \sin(\vartheta), z] d \vartheta$$

We could also simplify the $$T_m(z)$$ into $$T_m(z) = T(0, H_a - r_a, z)$$.

As we can see the difficulty is the coupling of 1D and 3D equation. @user21 suggested to transform $$T_f$$ into a 3D function and the problem can be bypassed. Although this method gives me results, the solution $$T_f$$ is not strictly 1D, thus the energy might not be conserved. I implemented @user21's suggestion with the following code.

(* constants *)
Lp = 0.25;
da = 0.02;
hi = 0.15;
ha = 0.15;
ks = 2.1;
alphaf = 1000;
alphai = 6;
alphaa = 0.6;
ra = da/2;
Hi = hi + ra;
Ha = ha + ra;
H = Hi + Ha;
cW = 4200;
rhoW = 1000;
kW = 0.6;
l = 1;
Vf = 10 Lp l/3600000;

(* mesh *)
pointsStructure = {{0, 0}, {Lp/2, 0}, {Lp/2, H}, {0, H}};
pointsPipeOuter = Table[ra {Cos[theta Degree], Sin[theta Degree]} + {0, Ha}, {theta, 90, -90, -20}];
{len1, len2} = Length /@ {pointsStructure, pointsPipeOuter};
contour = Table[{i, If[i == len1 + len2, 1, i + 1]}, {i, 1, len1 + len2}];
line1D = MeshRegion[Table[{i}, {i, 0, l, l/50.}], Line /@ Table[{i, i + 1}, {i, 50}]];
bmesh = ToBoundaryMesh["Coordinates" -> Join[pointsStructure, pointsPipeOuter], "BoundaryElements" -> {LineElement[contour]}];
mesh2D = ToElementMesh[bmesh, "MeshOrder" -> 1, MaxCellMeasure -> 5 10^-5];
mesh2D["Wireframe"];
region2D = MeshRegion[mesh2D["Coordinates"], Triangle /@ mesh2D["MeshElements"][[1, 1]]];
region3D = RegionProduct[region2D, line1D];
mesh3D = ToElementMesh[region3D(*,MaxCellMeasure\[Rule]0.05 10^-5*)]

(* equations and NDSolve *)
eq = {-Inactive[Div][{{-ks, 0, 0}, {0, -ks, 0}, {0, 0, -ks}}.
Inactive[Grad][t[x, y, z], {x, y, z}], {x, y, z}] ==
NeumannValue[alphai t[x, y, z], y == H] +
NeumannValue[alphaa t[x, y, z], y == 0] +
NeumannValue[alphaf (t[x, y, z] - tf[x, y, z]),
x^2 + (y - Ha)^2 == ra^2],
-cW rhoW Vf D[tf[x, y, z], z] == Pi da alphaf (tf[x, y, z] - t[x, y, z]),
DirichletCondition[tf[x, y, z] == 10, z == 0]};
{T, Tf} = NDSolveValue[eq, {t, tf}, Element[{x, y, z}, mesh3D]];
SliceContourPlot3D[Tf[x, y, z], "ZStackedPlanes", Element[{x, y, z}, mesh3D]] Any idea to solve this problem is appreciated.

Edit 1

One thing I tried is to force the 3 dimensional $$T_f$$ to be 1D function by modifying the Eq2 into

$$\Big[c_W \rho_W \dot{V}_f \frac{\partial T_f(x, y, z)}{\partial z} + 2 \alpha_f \pi r_a (T_f(x, y, z) - T(x, y, z)) \Big]^2 + (\frac{\partial T_f(x, y, z)}{\partial x})^2 + (\frac{\partial T_f(x, y, z)}{\partial y})^2 = 0$$

eq = {-Inactive[Div][{{-ks, 0, 0}, {0, -ks, 0}, {0, 0, -ks}}.
Inactive[Grad][t[x, y, z], {x, y, z}], {x, y, z}] ==
NeumannValue[alphai t[x, y, z], y == H] +
NeumannValue[alphaa t[x, y, z], y == 0] +
NeumannValue[alphaf (t[x, y, z] - tf[x, y, z]),
x^2 + (y - Ha)^2 == ra^2],
(Pi da alphaf (tf[x, y, z] - t[x, y, z]) + cW rhoW Vf D[tf[x, y, z], z])^2
+ D[tf[x, y, z], x]^2 + D[tf[x, y, z], y]^2 == 0,
DirichletCondition[tf[x, y, z] == 10, z == 0]};
{T, Tf} = NDSolveValue[eq, {t, tf}, Element[{x, y, z}, mesh3D]];


But the error FindRoot::nosol: Linear equation encountered that has no solution. occurs.

• See the solution to a similar problem. Choose the option that suits you. I will help. mathematica.stackexchange.com/questions/193789/… – Alex Trounev Aug 30 '19 at 7:53
• @AlexTrounev thank you for your advice. I'm trying to understand the code you wrote in that post. I found the Michael E2's approach quite interesting and easy. My concern is that my problem is a little bit more complicated, since the 1D equation contains the $T_m(z)$, which is the integral of the 3D solution or the solution at a specific point. I'm not sure if NDSolve can handle integral differential equation or use $T(0, H_a - r_a, z)$ in the equation. – 407PZ Aug 30 '19 at 9:24
• This is a piece of cake. On this site there are many such problems solved. Why do not you want to calculate the flow of fluid in the pipe and heat transfer? – Alex Trounev Aug 30 '19 at 11:29
• @AlexTrounev Because it would take too long for steady-state and transient case, and the temperature distribution is actually more important to me. A simplified method to calculate the convective heat transfer between the fluid and the pipe wall could save a huge amount of time. – 407PZ Aug 30 '19 at 11:34
• @AlexTrounev Could you give me advice on how to handle this integral differential equation? I'm thinking about calculating the $T(x, y, z)$ and $T_f(z)$ iteratively with NDSolve: solve $T(x, y, z)$ first and then solve $T_f(z)$, repeating the procedure until the results converges. I could imagine that it would work nice for the steady-state case. For the transient case however, it would be trickier. – 407PZ Aug 30 '19 at 11:44

We use the method of the false transient and the simple model Tf[z]=T[0,Ha-ra,z] to show the difference from the solution with the temperature Tf[x,y,z] defined in the volume.

Needs["NDSolveFEM"];(*constants*)Lp = 0.25;
da = 0.02;
hi = 0.15;
ha = 0.15;
ks = 2.1;
alphaf = 1000;
alphai = 6;
alphaa = 0.6;
ra = da/2;
Hi = hi + ra;
Ha = ha + ra;
H = Hi + Ha;
cW = 4200;
rhoW = 1000;
kW = 0.6;
l = 1;
Vf = 10 Lp l/3600000;

(*mesh*)
pointsStructure = {{0, 0}, {Lp/2, 0}, {Lp/2, H}, {0, H}};
pointsPipeOuter =
Table[ra {Cos[theta Degree], Sin[theta Degree]} + {0, Ha}, {theta,
90, -90, -20}];
{len1, len2} = Length /@ {pointsStructure, pointsPipeOuter};
contour =
Table[{i, If[i == len1 + len2, 1, i + 1]}, {i, 1, len1 + len2}];
line1D = MeshRegion[Table[{i}, {i, 0, l, l/50.}],
Line /@ Table[{i, i + 1}, {i, 50}]];
bmesh = ToBoundaryMesh[
"Coordinates" -> Join[pointsStructure, pointsPipeOuter],
"BoundaryElements" -> {LineElement[contour]}];
mesh2D = ToElementMesh[bmesh, "MeshOrder" -> 1,
MaxCellMeasure -> 10^-4];
mesh2D["Wireframe"];
region2D =
MeshRegion[mesh2D["Coordinates"],
Triangle /@ mesh2D["MeshElements"][[1, 1]]];
region3D = RegionProduct[region2D, line1D];
mesh3D = ToElementMesh[region3D(*,MaxCellMeasure\[Rule]0.05 10^-5*)]

(*equations and NDSolve*)
eq = {-Inactive[
Div][{{-ks, 0, 0}, {0, -ks, 0}, {0, 0, -ks}}.Inactive[Grad][
t[x, y, z], {x, y, z}], {x, y, z}] ==
NeumannValue[alphai t[x, y, z], y == H] +
NeumannValue[alphaa t[x, y, z], y == 0] +
NeumannValue[alphaf (t[x, y, z] - tf[x, y, z]),
x^2 + (y - Ha)^2 == ra^2], -cW rhoW Vf D[tf[x, y, z], z] ==
Pi da alphaf (tf[x, y, z] - t[x, y, z]),
DirichletCondition[tf[x, y, z] == 10, z == 0]};
{T, Tf} = NDSolveValue[eq, {t, tf}, Element[{x, y, z}, mesh3D]];


The method of the false transient

Tn[x_, y_, z_] := T[x, y, z]
Tfn[z_] := T[0, Ha - ra, z]
n = 3;
Do[Tn[i] =
NDSolveValue[{-Inactive[
Div][{{-ks, 0, 0}, {0, -ks, 0}, {0, 0, -ks}}.Inactive[Grad][
t[x, y, z], {x, y, z}], {x, y, z}] ==
NeumannValue[alphai t[x, y, z], y == H] +
NeumannValue[alphaa t[x, y, z], y == 0] +
NeumannValue[alphaf (t[x, y, z] - Tfn[i - 1][z]),
x^2 + (y - Ha)^2 == ra^2]}, t, Element[{x, y, z}, mesh3D]];
Tfn[i] = NDSolveValue[{-cW rhoW Vf D[tf[z], z] ==
Pi da alphaf (tf[z] - Tn[i][0, Ha - ra, z]), tf == 10},
tf, {z, 0, 1}];, {i, 1, n}]


The solution quickly converges. You can see the difference with the original solution.The difference is not very large as can be seen from the figures.

Plot[Evaluate[Table[Tfn[i][z], {i, 0, n}]], {z, 0, 1},
PlotLegends -> Automatic]

Table[ContourPlot[Tn[i][x, y, .5], {x, 0, .125}, {y, 0, .32},
ColorFunction -> Hue, Contours -> 20, PlotLegends -> Automatic], {i,
0, n}] • Thank you for your answer. This solution seems quite good for the steady-state case but would it also be suitable for the transient case? – 407PZ Aug 30 '19 at 13:18
• For transient, this method works even better than NDSolve[] in the new version. I tested the method on an unsteady flow around a cylinder and on problems of thermal convection, see community.wolfram.com/groups/-/m/t/1433064 – Alex Trounev Aug 30 '19 at 13:25
• I have some difficulty for understanding the code you wrote for the unsteady flow with NS-equation. Could you please give me a brief how you implemented the explicit euler with NDSolve[]? – 407PZ Sep 2 '19 at 7:48
• Mathematica has a good linear FEM solver. I used it to solve non-linear non-stationary problems. For this, I used a simple algorithm similar to an explicit Euler (in time). At each step, linear FEM solver and interpolation are used there. This is due to the fact that the time derivative enters the NSE linearly. Therefore, constructing a numerical algorithm is not difficult. But there are tricks. Therefore, I tested the algorithm very carefully. – Alex Trounev Sep 2 '19 at 13:01