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.
NDSolvecan handle integral differential equation or use $T(0, H_a - r_a, z)$ in the equation. $\endgroup$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. $\endgroup$