Approach to analytically solve a Coupled system of PDE in Mathematica [Heat transfer in cylindrical coordinates]

I have the following two PDEs, which describe steady-state coupled heat transport between a externally heated axi-symmetric solid body (Eq. 1, $$T(r,z)$$) and a fluid (Eq. 2, $$t(z)$$) flowing inside it

$$\frac{\partial^2 T}{\partial r^2}+\frac{1}{r}\frac{\partial T}{\partial r}+\frac{\partial^2 T}{\partial z^2}=0 \tag1$$

$$\frac{\partial t}{\partial z}+\alpha(t-T(r_1,z))=0 \tag2$$

Eq. (1) is defined in the domain $$r\in[r_1,r_2]$$ where $$r_1$$ and $$r_2$$ describe the inner and outer radii of the cylinder and $$z\in[0,L]$$ where $$L$$ is the length of the cylinder. The boundary conditions for Eq. (1) are $$\frac{\partial T(r,0)}{\partial z}=\frac{\partial T(r,L)}{\partial z}=0 \tag3$$

$$\frac{\partial T(r_2,z)}{\partial r}=\gamma \tag4$$

$$\frac{\partial T(r_1,z)}{\partial r}=\beta(T(r_1,z)-t) \tag5$$

For Eq. (2) it is known that $$t(z=0)=t_{in}$$

$$\alpha,\beta,\gamma,t_{in}$$ are known constants. It seems the solid and fluid temperatures are coupled through the B.C. at $$r=r_1$$ (solid-fluid interface, Robin condition).

Any suggestion on how to approach this problem analytically in Mathematica is appreciated. I get that this is not a Mathematica related question but I have had some excellent feedback on my earlier questions which have helped me find better solution methodologies.

Following Bill Watts' answer, I took some realistic parameters.

These constants correspond to a copper circular channel (thermal conductivity = 390 W/mK) with inner and outer radii of $$1 mm$$ and $$2 mm$$ respectively in which fluid enters with a velocity of $$0.0333 m/s$$. The cylinder is heated externally by a heat flux of $$8000 W/m^2$$ and the heat transfer coefficient is $$2000 W/m^2 K$$

which give

\[Alpha] = 28.852; \[Beta] = 5.128; \[Gamma] = 20.51; tin = 300; L = 0.03; r1 = 0.001; r2 = 0.002;

and on plotting the boundary condition $$(5)$$, the discrepancy seems to have reduced

For the same set of parameters except with $$r_2 = 5 mm$$, the discrepancy almost vanishes

• What have you tried? Where is Mathematica code? Apr 28, 2020 at 15:49
• @MariuszIwaniuk I am sorry to say that I have not tried anything still in Mathematica. The maximum I did was express $t$ as an integral function of $T$ from $(2)$ and substituted in the boundary condition $(5)$ which transforms it into $$\frac{\partial T}{\partial r} \vert_{r=r_1}=\beta\Bigg[T-\alpha e^{-\alpha z}\bigg(\int_0^z e^{\alpha s} T(r,s) \mathrm{d}s + \frac{T_{fi}}{\alpha} \bigg)\Bigg]$$ I could not figure out how to proceed after this. Apr 28, 2020 at 16:26
• I guess in Eq(2) the $T$ should be $T(r_1,z)$? May 1, 2020 at 1:36
• Just to be sure, a numeric approach would not be helpful right? Second question: Can any of the constants assumed to be small? It would be very helpful to start with the assumption that $\alpha$ is small. In that case a pertubational expansion in $\alpha$ might be feasible. A fully analytic treatment of a PDE is typically the exception not the default case. One often starts by analyzing certain special cases analytically and then proceeds by using those to verify the more general numeric solutions.
– Max1
May 1, 2020 at 6:03
• The new system is harder to solve analytically compared to the previous one mentioned in this and this post, because $(2)$ and $(5)$ is equivalent to a b.c. involving first order derivative of $z$, thus finite Fourier transform won't help here. There might exists other integral transform that's suitable for this type of b.c., but searching such transform is beyond my reach. May 1, 2020 at 6:37

This solution is not perfect, but I will throw it out there anyway in case anyone has an interest to improve it.
Use separation of variables

Clear["Global*"]


Work on the T equation first

pde = D[T[r, z], r, r] + (1/r)*D[T[r, z], r] + D[T[r, z], z, z] == 0


Separation by multiples

T[r_, z_] = R[r] Z[z]

pde/T[r, z] // Expand
(*R''[r]/R[r] + R''[r]/(r R[r]) + Z''[z]/Z[z] == 0*)


Choose the z equation such that it is sinusoidal in z due to the given boundary conditions.

zeq = Z''[z]/Z[z] == -a^2

DSolve[zeq, Z[z], z] // Flatten

Z[z_] = Z[z] /. % /. {C[1] -> c1, C[2] -> c2}
(*c1 Cos[a z] + c2 Sin[a z]*)


Now the R equation

req = R''[r]/R[r] + R'[r]/(r R[r]) == a^2

DSolve[req, R[r], r] // Flatten

R[r_] = (R[r] /. % /. {C[1] -> c3, C[2] -> c4})
(*c3 BesselJ[0, I a r] + c4 BesselY[0, -I a r]*)


I don't know why Mathematica always insists on complex solutions for this equation. Convert by:

FullSimplify[FunctionExpand[R[r], r > 0]] // Collect[#, BesselI[0, a r]] &


Consolidate constants

R[r_] = % /. {Coefficient[%, BesselI[0, a r]] -> c3, Coefficient[%, BesselK[0, a r]] -> c4}
(*c3 BesselI[0, a r] + c4 BesselK[0, a r]*)


As usual with the diffusivity equation we don't have enough pieces with separation by multiplication. Now separate by addition.

T[r_, z_] = Rp[r] + Zp[z]

pde
(*Rp''[r] + Rp'[r]/r + Zp''[z] == 0*)

zpeq = Zp''[z] == b

DSolve[zpeq, Zp[z], z] // Flatten

Zp[z_] = Zp[z] /. % /. {C[1] -> c5, C[2] -> c6}
(*(b z^2)/2 + c5 + c6 z*)

rpeq = Rp''[r] + Rp'[r]/r + b == 0

DSolve[rpeq, Rp[r], r] // Flatten

Rp[r_] = Rp[r] /. % /. {C[1] -> c7, C[2] -> 0}
(*c7 Log[r] - (b r^2)/4*)


I chose C[1] to be zero because we don't need two constant terms. Put it all together:

T[r_, z_] = R[r] Z[z] + Rp[r] + Zp[z]
(c1 Cos[a z] + c2 Sin[a z]) (c3 BesselI[0, a r] + c4 BesselK[0, a r]) - (b r^2)/4 + (b z^2)/2 + c5 + c6 z + c7 Log[r]


Check

pde // FullSimplify
(*True*)


Apply the boundary conditions

(D[T[r, z], z] /. z -> 0) == 0
(*a c2 (c3 BesselI[0, a r] + c4 BesselK[0, a r]) + c6 == 0*)


so

c2 = 0
c6 = 0


and consolidate constants

c1 = 1

(D[T[r, z], z] /. z -> L) == 0
(*b L - a Sin[a L] (c3 BesselI[0, a r] + c4 BesselK[0, a r]) == 0*)


from which

b = 0


and to make the Sin zero:

a = (n π)/L


with

$Assumptions = n ∈ Integers  T becomes an infinite series in n, but we will leave off the sum for now so MMa won't constantly try to evaluate it. (D[T[r, z], r] /. r -> r2) == γ (*Cos[(π n z)/L] ((π c3 n BesselI[1, (n π r2)/L])/L - (π c4 n BesselK[1, (n π r2)/L])/L) + c7/r2 == γ*)  We can satisfy by c4 = c4 /. Solve[Coefficient[%[[1]], Cos[(\[Pi] n z)/L]] == 0, c4][[1]] (*(c3 BesselI[1, (n π r2)/L])/BesselK[1, (n π r2)/L]*)  and c7 = c7 /. Solve[c7/r2 == γ, c7][[1]] (*γ r2*) T[r, z] // Collect[#, c3] &  Check out the solution when n = 0. BesselK is unbounded with zero arguments, so take the limit. Limit[T[r, z], n -> 0] (*c3 + c5 + γ r2 Log[r]*)  Note that c5 is the equivalent c3 constant when n = 0 in the Fourier series. We only need to keep one of them, so for n = 0 T0[r_, z_] = % /. c3 -> 0  For general n Tn[r_, z_] = T[r, z] - T0[r, z] // Simplify  Now work on the differential equation for t. pdet = (t'[z] + α (t[z] - T[r1, z]) == 0)  General n pde2 = (tn'[z] + α (tn[z] - Tn[r1, z]) == 0) (DSolve[pde2, tn[z], z] // Flatten) tn[z_] = (tn[z] /. % /. C[1] -> c8)  The outputs are getting a little long to show here. For n = 0 pde20 = t0'[z] + α (t0[z] - T0[r1, z]) == 0 DSolve[pde20, t0[z], z] // Flatten t0[z_] = t0[z] /. % /. C[1] -> c80 (*c5 + c80 E^(α (-z)) + γ r2 Log[r1]*)  Now apply the initial condition t[0] == tin Do this by setting the part contain n to zero, and set the constant part to tin. c8 = c8 /. Solve[tn[0] == 0, c8][[1]] c80 = c80 /. Solve[t0[0] == tin, c80][[1]] tn[z_] = tn[z] // Simplify t0[z] // Simplify t[z_] = t0[z] + tn[z]  where it is understood that the part containing n is the sum over n from 1 to infinity. Check the t solution. pdet // Simplify (*True*)  Apply the final bc on general n and n = 0 separately using the orthogonality of Cos[(π n z)/L]. The final boundary condition. bcf = (D[T[r, z], r] /. r -> r1) == β (T[r1, z] - t[z])  For n = 0 Limit[bcf[[1]], n -> 0] (*(γ r2)/r1*) Limit[bcf[[2]], n -> 0] (*β E^(α (-z)) (c3 + c5 + γ r2 Log[r1] - tin)*)  Again, c5 is just the constant term in the fourier series when n = 0, so we don't need both it and c3. bcfn0 = % == %% /. c5 + c3 -> c30 (*β E^(α (-z)) (c30 + γ r2 Log[r1] - tin) == (γ r2)/r1*)  Use orthogonality Integrate[bcfn0[[1]], {z, 0, L}] == Integrate[bcfn0[[2]], {z, 0, L}] c5 = c30 /. Solve[%, c30][[1]] // Simplify  General n ortheq = Integrate[bcf[[1]]*Cos[(n*Pi*z)/L], {z, 0, L}] == Integrate[bcf[[2]]*Cos[(n*Pi*z)/L], {z, 0, L}] c3 = c3 /. Solve[%, c3][[1]] // Simplify  Simplify everything. t0[z_] = t0[z] // Simplify tn[z_] = tn[z] // Simplify T0[r_, z_] = T0[r, z] // Simplify Tn[r_, z] = Tn[r, z] // Simplify  Plug in numbers α = 1/10; β = 1/10; γ = 1; tin = 1; L = 10; r1 = 1; r2 = 2;  I am using exact numbers so I can use lots of terms in the Fourier series if necessary. For calculation, add an additional argument used for the number of terms in the series. T[r_, z_, mm_] := T0[r, z] + Sum[Tn[r, z], {n, 1, mm}] t[z_, mm_] := t0[z] + Sum[tn[z], {n, 1, mm}]  Of course mm should actually be infinity, but we will use a finite series for calculation. And the derivatives dtdz[Z_, mm_] := (D[t0[z], z] /. z -> Z) + Sum[D[tn[z], z] /. z -> Z, {n, 1, mm}] dTdr[R_, z_, mm_] := (D[T0[r, z], r] /. r -> R) + Sum[D[Tn[r, z], r] /. r -> R, {n, 1, mm}] dTdz[r_, Z_, mm_] := (D[T0[r, z], z] /. z -> Z) + Sum[D[Tn[r, z], z] /. z -> Z, {n, 1, mm}]  Compiling the expressions will speed up the calculations, but compiling is limited to machine precision values. For checking I don't want that restriction. Make some plots. T at a few values of z Plot[{Evaluate[T[r, 0, 50]], Evaluate[T[r, L/2, 50]], Evaluate[T[r, L, 50]]}, {r, r1, r2}]  Plot3D[Evaluate[T[r, z, 50]], {r, r1, r2}, {z, 0, L}, PlotRange -> All]  Check t[0] == tin (*True*)  Plot of t Plot[Evaluate[t[z, 50]], {z, 0, L}]  The t pde Steps = 200 Plot[Evaluate[dtdz[z, Steps] + α (t[z, Steps] - T[r1, z, Steps])], {z, 0, L}, PlotRange -> All]  Pretty close to zero. The boundary at r2. Plot[Evaluate[dTdr[r, z, 20] /. r -> r2] - γ, {z, 0, L}]  The final boundary condition. Plot[{Evaluate[dTdr[r, z, 50] /. r -> r1], Evaluate[β (T[r1, z, 50] - t[z, 50])]}, {z, 0, L}, PlotRange -> {1.5, 2.8}] ` All the other checks are good, but these two plots should lie on top of each other. And while they are not way off, I think the difference is too large to just be numerical error. I invite anyone with an interest in this type of problem to review this solution for improvement. • In my experience the diffusivity equation generally requires both multiplication and addition forms of variable separation. You can't get there with just multiplication. I don't have a reference, but I've been doing it since the 70's. Zeroing$c2$and$c6$is necessary to make the first boundary condition$True$. May 8, 2020 at 5:02 • That comes directly from the contribution of$t0(z)\$ in bcf. You shouldn't have lost that. May 8, 2020 at 7:37
• It is possible that your LHS and RHS of the equations could be slightly different from mine. All the same terms should be there though, and they would solve to the same final results. May 8, 2020 at 7:41
• Interesting that different parameters would get a closer match. I'm glad it worked out. May 8, 2020 at 17:35
• I originally learned from taking electrodynamics in grad school out of J.D. Jackson's "Classical Electrodynamics" book. It is an excellent reference, but be prepared to spend a lot of time understanding it. Otherwise, I worked with the diffusivity equation a lot in the oil field, since the reservoir pressure due to fluid flow in porous media is governed by that equation. This is the first time I worked the diffusivity equation simultaneously with another differential equation, however. May 8, 2020 at 20:39