4
$\begingroup$

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 enter image description here

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

enter image description here

Improve this question
$\endgroup$
9
  • 2
    $\begingroup$ What have you tried? Where is Mathematica code? $\endgroup$
    – Mariusz Iwaniuk
    Apr 28, 2020 at 15:49
  • $\begingroup$ @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. $\endgroup$
    – Avrana
    Apr 28, 2020 at 16:26
  • $\begingroup$ I guess in Eq(2) the $T$ should be $T(r_1,z)$? $\endgroup$
    – xzczd
    May 1, 2020 at 1:36
  • 3
    $\begingroup$ 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. $\endgroup$
    – Max1
    May 1, 2020 at 6:03
  • 1
    $\begingroup$ 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. $\endgroup$
    – xzczd
    May 1, 2020 at 6:37

1 Answer 1

Reset to default
3
+50
$\begingroup$

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

enter image description here

Plot3D[Evaluate[T[r, z, 50]], {r, r1, r2}, {z, 0, L}, PlotRange -> All]

enter image description here

Check

t[0] == tin
(*True*)

Plot of t

Plot[Evaluate[t[z, 50]], {z, 0, L}]

enter image description here

The t pde

Steps = 200

Plot[Evaluate[dtdz[z, Steps] + α (t[z, Steps] - T[r1, z, Steps])], {z, 0, L}, PlotRange -> All]

enter image description here

Pretty close to zero.

The boundary at r2.

Plot[Evaluate[dTdr[r, z, 20] /. r -> r2] - γ, {z, 0, L}]

enter image description here

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

enter image description here

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.

Improve this answer
$\endgroup$
16
  • 1
    $\begingroup$ 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$. $\endgroup$
    – Bill Watts
    May 8, 2020 at 5:02
  • 1
    $\begingroup$ That comes directly from the contribution of $t0(z)$ in bcf. You shouldn't have lost that. $\endgroup$
    – Bill Watts
    May 8, 2020 at 7:37
  • 1
    $\begingroup$ 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. $\endgroup$
    – Bill Watts
    May 8, 2020 at 7:41
  • 1
    $\begingroup$ Interesting that different parameters would get a closer match. I'm glad it worked out. $\endgroup$
    – Bill Watts
    May 8, 2020 at 17:35
  • 1
    $\begingroup$ 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. $\endgroup$
    – Bill Watts
    May 8, 2020 at 20:39

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.

Not the answer you're looking for? Browse other questions tagged or ask your own question.