Modelling Heat Transfer for Pulsating Laminar Duct Flow

[rafisondi]

Good evening everybody.
I am new here I hope it is okay to ask this question here.

As part of a project I am trying to examine the effect of a pulsating laminar duct flow on the heat transfer behavior for constant heat flux at the wall for a Newtonian incompressible fluid.

As such I am assuming an oscillating pressure gradient dp/dz = A*exp( i*ω*t), A being a real constant denoting the pressure amplitude, ω being the radial velocity and time t > 0. Using conservation of mass and the Navier Stokes equations, I get a quasi steady state flow profile that is only dependent on time and the radius w(r,t) in axial direction x which only depends on time and radius r. I attached a plot graph of my velocity profile vs radius (Radius R of duct = 0.01m). The different colored graphes are ploted for different times t_i.

I essentially only evaluated the profile 1dimensional over the radius since we are hydrodynamically developed and can neglect the azimutal direction bc of symmetry arguments.

My next step would have been to evaluate the temperature profile T(r,x,t) using the energy conservation assuming constant heat flux q over the duct wall area. Assuming no dissipation of energy and neglecting the thermal diffusion in axial direction (∂^2T/∂x^2) to be small compared to radial diffusion I get

ρ*c_p (∂T/∂t + u * ∂T/∂x ) = k (1/r *∂/∂r(r*∂T/∂r ) )

My biggest issue right now lies in the red term ∂T/∂x. I don't know how to "deal" with it... I can't really neglect it since that would knock out my fluid velocity. I somehow need to add some kind of "value" or expression to it, such that I can later numerically evaluate this PDE. I somehow feel like this should be actually very simple to solve but I can't seem can't seem to wrap my head around this problem.

Hopefully I framed my problem clearly. I am still a student and this is my first time working on such a problem so please excuse me possibly not being very exact with my formulations.

Anyways thanks in advance and have a nice day.

 

[Chestermiller]

Are you examining the oscillatory behavior superimposed on a constant pressure gradient, or are you really analyzing a flow with no mean axial velocity and a purely oscillatory pressure gradient?

 

[rafisondi]

Are you examining the oscillatory behavior superimposed on a constant pressure gradient, or are you really analyzing a flow with no mean axial velocity and a purely oscillatory pressure gradient?

The plot I attached in my previous post is done by examining the oscillatory behavior with no mean velocity on the flow. But I am intending to superimpose a constant pressure gradient on the flow.

Just to give you some context, the endgoal would be to analyse heat flow in a blood vessel with a heart causing the periodic pressure behavior.

 

[Chestermiller]

What makes you think that the temperature will vary with axial position? Does the velocity vary with axial position? I assume you are modeling an infinitely long tube, right?

 

[rafisondi]

What makes you think that the temperature will vary with axial position? Does the velocity vary with axial position? I assume you are modeling an infinitely long tube, right?

Well my inital thought comes from my understanding of the steady poiseulle flow temperature profile with constant heat flux on the walls (and as you mentionend assumed infinitly long pipe). In that case we should have a constant heat gradient in axial direction. I understood this to be the fluid heating up along the pipe.

In the case of a steady Poiseulle flow, you could show that axial temperature gradient to be constant and the energy equation I formulated above could be solved analytically.

I don't expect there to be an analytical solution for my PDE but some kind of argument for the gradient, such that I could solve for the temperature profile along the radial axis numerically.

 

[Chestermiller]

Well my inital thought comes from my understanding of the steady poiseulle flow temperature profile with constant heat flux on the walls (and as you mentionend assumed infinitly long pipe). In that case we should have a constant heat gradient in axial direction. I understood this to be the fluid heating up along the pipe.

In the case of a steady Poiseulle flow, you could show that axial temperature gradient to be constant and the energy equation I formulated above could be solved analytically.

I don't expect there to be an analytical solution for my PDE but some kind of argument for the gradient, such that I could solve for the temperature profile along the radial axis numerically.

So you are talking about an entrance flow where, up to x = 0, there is no heating. Because, if you are talking about an infinitely long pipe, the temperature profile will be purely radial, and there will be no axial temperature gradient. What are your BCs on x if you have only oscillation and no underlying steady flow?

 

[rafisondi]

So you are talking about an entrance flow where, up to x = 0, there is no heating.

You are right there is no heating for x<0. As BCs I guess I only have a uniform entrance temperature T(x=0) = T_in.

 

[Chestermiller]

You are right there is no heating for x<0. As BCs I guess I only have a uniform entrance temperature T(x=0) = T_in.

You are supplying very limited piecemeal information.

Are you trying to solve for the stationary oscillatory state only or are you also trying to also include the transient behavior leading up to the stationary oscillatory state? In any case, what does your solution for the stationary oscillatory state look like for the velocity?

 

[rafisondi]

I am sorry for being inaccurate. I really appreciate you for taking your time with my problem.

Are you trying to solve for the stationary oscillatory state only or are you also trying to also include the transient behavior leading up to the stationary oscillatory state?

I am only solving for the steady oscillatory state without transient behavior.

In any case, what does your solution for the stationary oscillatory state look like for the velocity?

I attached a velocity profile for my stationary oscillatory flow velocity profile at my initial question. The plot depicts the velocity profile oscillating along the radial axis. I think this is what you are referring to, right?

 

[Chestermiller]

I am sorry for being inaccurate. I really appreciate you for taking your time with my problem.
I am only solving for the steady oscillatory state without transient behavior.
I attached a velocity profile for my stationary oscillatory flow velocity profile at my initial question. The plot depicts the velocity profile oscillating along the radial axis. I think this is what you are referring to, right?

So you are expressing the axial velocity as $$v=A(\omega, r)\cos{\omega t}+B(\omega, r)\sin{\omega t}$$I'd like to see how you solved for A and B. I'd also like what the equations look like when you have a constant pressure gradient plus an oscillatory pressure gradient superimposed (in component form).

I guess you are also assuming that the viscosity is independent of temperature, correct?

 

[rafisondi]

I guess you are also assuming that the viscosity is independent of temperature, correct?

Yes you are right, I am assuming viscosity and density to be temperature independent.

The deriviation of the velocity profile can be done like this:

Looking at the axial component of the Navier Stokes equation:
\begin{equation}
\rho \frac{\partial}{\partial t}w = -\frac{\partial}{\partial z} p + \mu \frac{1}{r} \frac{\partial }{\partial r} \left(r \frac{\partial}{\partial r} w \right)
\end{equation}

Now assuming \begin{equation}\frac{\partial}{\partial z} p = A \cdot e^{i \omega t}\end{equation} with A being the pressure amplitude we can try to solve this using the ansatz that our veloctiy also oscillates with the same frequency w:

\begin{equation}
w(r,t) = \hat{w} \cdot e^{iwt}
\end{equation}

Plugin all in yields:

\begin{equation}
\mu \frac{1}{r} \frac{\partial }{\partial r} \left(r \frac{\partial}{\partial r} \hat{w} \right) - i \omega \rho \hat{w} = A
\end{equation}

Solving this ODE is a bit stressful but the particular solution is \begin{equation}w_p = i \frac{A}{\omega \rho} \end{equation}, while to solve the homogenous solution one has to notice that the Homogenous ODE takes the form of a Bessel equation of zeroth order. This yields the solution:

\begin{equation}
\hat{w}_h = i\frac{A}{\omega \rho } + c_1 J_0(\hat{r}) +c_2 I_0(\hat{r}) = 0
\end{equation}

Using the no slip at the wall condition and no singularites (C2 = 0) yields

\begin{equation}
w(r,t) = i\frac{A}{\rho \omega} \left[1 - \frac{J_0(i^{3/2}\alpha (r/a))}{J_0(i^{3/2}\alpha)} \right] \cdot e^{i \omega t}
\end{equation}

with
\begin{equation}
\alpha = L \sqrt{\frac{\omega}{\nu}}
\end{equation}
where L is usually the pipe radius.

This is essentially u = A(ω,r)cosωt+B(ω,r)sinωt as you wrote it. I hope this deriviation is clear. This flow profile was essentially derived by John R. Womersley(though I think he derived it differently) and is labeled as the Womersley flow profile in case you want to google it.

 

[Chestermiller]

Yes you are right, I am assuming viscosity and density to be temperature independent.

The deriviation of the velocity profile can be done like this:

Looking at the axial component of the Navier Stokes equation:
\begin{equation}
\rho \frac{\partial}{\partial t}w = -\frac{\partial}{\partial z} p + \mu \frac{1}{r} \frac{\partial }{\partial r} \left(r \frac{\partial}{\partial r} w \right)
\end{equation}

Now assuming \begin{equation}\frac{\partial}{\partial z} p = A \cdot e^{i \omega t}\end{equation} with A being the pressure amplitude we can try to solve this using the ansatz that our veloctiy also oscillates with the same frequency w:

\begin{equation}
w(r,t) = \hat{w} \cdot e^{iwt}
\end{equation}

Plugin all in yields:

\begin{equation}
\mu \frac{1}{r} \frac{\partial }{\partial r} \left(r \frac{\partial}{\partial r} \hat{w} \right) - i \omega \rho \hat{w} = A
\end{equation}

Solving this ODE is a bit stressful but the particular solution is \begin{equation}w_p = i \frac{A}{\omega \rho} \end{equation}, while to solve the homogenous solution one has to notice that the Homogenous ODE takes the form of a Bessel equation of zeroth order. This yields the solution:

\begin{equation}
\hat{w}_h = i\frac{A}{\omega \rho } + c_1 J_0(\hat{r}) +c_2 I_0(\hat{r}) = 0
\end{equation}

Using the no slip at the wall condition and no singularites (C2 = 0) yields

\begin{equation}
w(r,t) = i\frac{A}{\rho \omega} \left[1 - \frac{J_0(i^{3/2}\alpha (r/a))}{J_0(i^{3/2}\alpha)} \right] \cdot e^{i \omega t}
\end{equation}

with
\begin{equation}
\alpha = L \sqrt{\frac{\omega}{\nu}}
\end{equation}
where L is usually the pipe radius.

This is essentially u = A(ω,r)cosωt+B(ω,r)sinωt as you wrote it. I hope this deriviation is clear. This flow profile was essentially derived by John R. Womersley(though I think he derived it differently) and is labeled as the Womersley flow profile in case you want to google it.

Good, now let's see your differential equation(s) for the heat balance, assuming that, superimposed on the pressure gradient $A^0$ there is an oscillatory component of amplitude $A^1$, letting the temperature be $T^0+T^1$ and the velocity be $w^0+w^1$, with the zero superscripts signifying the steady state and the superscript 1 signifying the oscillatory perturbation.

 

[rafisondi]

Firstoff a little correction to my previous post:

This yields the solution:

​

(6)\begin{equation} \hat{w}_h = i\frac{A}{\omega \rho } + c_1 J_0(\hat{r}) +c_2 I_0(\hat{r}) = 0 \end{equation}
Using the no slip at the wall condition and no singularites (C2 = 0) yields

(7)\begin{equation} w(r,t) = i\frac{A}{\rho \omega} \left[1 - \frac{J_0(i^{3/2}\alpha (r/a))}{J_0(i^{3/2}\alpha)} \right] \cdot e^{i \omega t} \end{equation}
with

​

(8)\begin{equation} \alpha = L \sqrt{\frac{\omega}{\nu}} \end{equation} where L is usually the pipe radius.

I skipped several steps here: using the no slip boundary condition( Velocity has to be zero) we find an expression for c1. We set c2 = 0 because of a singularity at r = 0. And then we combine the homogenous and particular solution together to get:

\begin{equation} w(r,t) = w_h + w_p = i\frac{A}{\rho \omega} \left[1 - \frac{J_0(i^{3/2}\alpha (r/a))}{J_0(i^{3/2}\alpha)} \right] \cdot e^{i \omega t} \end{equation}

I wanted to add this correction in case someone wanted to read my deriviation and got confused.

Good, now let's see your differential equation(s) for the heat balance

For a fully developed convection heat transfer the governing equation should look like this:

​

\begin{equation} \rho c_p (\frac{\partial T}{\partial t} w\frac{\partial T}{\partial x} + u_r \frac{\partial T}{\partial r} + u_\theta (\frac{1}{r}\frac{\partial T}{\partial \theta}) )

= \frac{1}{r}\frac{\partial}{\partial r} ( r k \frac{\partial T}{\partial r} ) + \frac{\partial }{\partial x}(k \frac{\partial T}{\partial x}) + \frac{1}{r^2} \frac{\partial}{\partial \theta} ( k \frac{\partial T}{\partial \theta} )

\end{equation}

with fluid density ρ, fluid heat capacity c_p and heat conductivity coefficient k.
u_r (radial velocity) is zero and $$u_\theta (\frac{1}{r}\frac{\partial T}{\partial \theta})$$ (Aximutal-Velocity) can be neglected using symmetry arguments. Same for the last term on the right side.

One can further ague that the temperature diffusion in axial direction can be neglected since it is very small compared to the radial direction:
$$\frac{1}{r}\frac{\partial}{\partial r} ( r k \frac{\partial T}{\partial r} ) >> \frac{\partial }{\partial x}(k \frac{\partial T}{\partial x})$$

further simplifying to

\begin{equation} \rho c_p (\frac{\partial T}{\partial t} +w \frac{\partial T}{\partial x} )= \frac{1}{r}\frac{\partial}{\partial r} ( r k \frac{\partial T}{\partial r} ) \end{equation}

My Boundary conditions read as follows:
\begin{equation}
r = 0: \frac{\partial T}{\partial r} = 0 ;
\end{equation}

and
\begin{equation}
r = R:k \frac{\partial T}{\partial r} = q"_{wall} ;
\end{equation}

R being the Radius of the pipe and $$q"_w$$ the heat flux over the wall.

-------------------

This is how far I got. I have not yet thought of assuming Temperature T to be deconstructed as T0+T1 a steady and an ocillatory component. I would know how to solve the steady part introducing the nondimensionalized temperature Θ

\begin{equation}
\Theta = \frac{T - T_{wall}}{T_{mean} - T_{wall}}
\end{equation}

But I would need some time to actually compute everything. But this would still leave open on how to solve the transient Temperature part.
I hope all of this makes sense, I tried to be as thorough as possible.

 

[Chestermiller]

The is good. Now, factor out k from the rhs of Eqn. 22 of the above post and combine it with $\rho$ and $C_p$ to make the thermal diffusivity $\alpha$ as a multiplier on the rhs. Next, make the substitutions I indicated at the bottom of post #12, and linearize with respect to the oscillatory perturbation by neglecting products of superscript 1 quantities. What do you get for the DE and BC's?

 

[rafisondi]

combine it with ρ\rho and CpC_p to make the thermal diffusivity α\alpha as a multiplier on the rhs

Okay so here is the equation written with thermal diffusivity:

\begin{equation}
\frac{\partial T}{\partial t} + w \cdot \frac{\partial T}{\partial x} = \alpha \cdot \frac{1}{r}\frac{\partial}{\partial r} ( r \frac{\partial}{\partial r} T )
\end{equation}

Next, make the substitutions I indicated at the bottom of post #12

Inserting into the equation above and rearranging\begin{equation}
\frac{\partial (T^1)}{\partial t} = \alpha \cdot \frac{1}{r}\frac{\partial}{\partial r} ( r \frac{\partial}{\partial r} (T^0+T^1) ) - (w^0+w^1) \cdot \frac{\partial (T^0+T^1)}{\partial x}
\end{equation}

-----

linearize with respect to the oscillatory perturbation

Honestly, I am not really sure on how to do so... If I understand correctly you want to linearize the PDE at the steady state $$T^0$$

Do you want me to partially deriviate the LHS for x and r evaluated at steady state and expanded by the coordinates?

\begin{align}
\frac{\partial (T^1)}{\partial t} =
\begin{pmatrix}
\begin{bmatrix}
\frac{\partial f_{LHS}}{\partial r}_{ss} , \frac{\partial f_{LHS}}{\partial r}_{ss} \\
\end{bmatrix} \cdot
\begin{bmatrix}
r \\
x
\end{bmatrix}
\end{pmatrix}
\end{align}

Again thanks in advance

 

[Chestermiller]

Okay so here is the equation written with thermal diffusivity:

\begin{equation}
\frac{\partial T}{\partial t} + w \cdot \frac{\partial T}{\partial x} = \alpha \cdot \frac{1}{r}\frac{\partial}{\partial r} ( r \frac{\partial}{\partial r} T )
\end{equation}
Inserting into the equation above and rearranging\begin{equation}
\frac{\partial (T^1)}{\partial t} = \alpha \cdot \frac{1}{r}\frac{\partial}{\partial r} ( r \frac{\partial}{\partial r} (T^0+T^1) ) - (w^0+w^1) \cdot \frac{\partial (T^0+T^1)}{\partial x}
\end{equation}

-----
Honestly, I am not really sure on how to do so... If I understand correctly you want to linearize the PDE at the steady state $$T^0$$

Do you want me to partially deriviate the LHS for x and r evaluated at steady state and expanded by the coordinates?

\begin{align}
\frac{\partial (T^1)}{\partial t} =
\begin{pmatrix}
\begin{bmatrix}
\frac{\partial f_{LHS}}{\partial r}_{ss} , \frac{\partial f_{LHS}}{\partial r}_{ss} \\
\end{bmatrix} \cdot
\begin{bmatrix}
r \\
x
\end{bmatrix}
\end{pmatrix}
\end{align}

Again thanks in advance

If I linearize Eqn. 27 with respect to the transient perturbation, I get:

$$\frac{\partial T^1}{\partial t}+w^0\frac{\partial T^0}{\partial x}+w^0\frac{\partial T^1}{\partial x}+w^1\frac{\partial T^0}{\partial x}=\alpha \frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial T^0}{\partial r}\right)+\alpha \frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial T^1}{\partial r}\right)$$ But the steady state solution satisfies:
$$w^0\frac{\partial T^0}{\partial x}=\alpha \frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial T^0}{\partial r}\right)$$Therefore, the transient perturbation must satisfy:
$$\frac{\partial T^1}{\partial t}+w^0\frac{\partial T^1}{\partial x}=-w^1\frac{\partial T^0}{\partial x}+\alpha \frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial T^1}{\partial r}\right)$$where $w^0$, $w^1$, and $T_0$ are already known.

OK so far? If so, now do the same thing with the boundary conditions.

 

[rafisondi]

Sorry about this late reply, I was lately a bit short on time.

OK so far? If so, now do the same thing with the boundary conditions.

I am a bit confused, to be honest...

If I linearize Eqn. 27 with respect to the transient perturbation, I get:

​

\begin{equation}
\frac{\partial T^1}{\partial t}+w^0\frac{\partial T^0}{\partial x}+w^0\frac{\partial T^1}{\partial x}+w^1\frac{\partial T^0}{\partial x}=\alpha \frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial T^0}{\partial r}\right)+\alpha \frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial T^1}{\partial r}\right)
\end{equation}

But the steady-state solution satisfies:

​

\begin{equation}
w^0\frac{\partial T^0}{\partial x}=\alpha \frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial T^0}{\partial r}\right)

\end{equation}

Pls, correct me if I am wrong but if I understand everything correctly you expanded my equation (27) with your argument, neglected the term with transient temperature, and used the steady-state solution as an additional condition that equation (27) has to fulfill.

Unfortunately, I do not really understand how to apply this method to my boundary conditions. I substituted BCs (23) and (24) with your arguments:\begin{equation} r = 0: \frac{\partial T^1 }{\partial r} + \frac{\partial T^0 }{\partial r} = 0 ; \end{equation}

and

\begin{equation} r = R:k \frac{\partial T^1 }{\partial r} + \frac{\partial T^0 }{\partial r} = q"_{wall} ; \end{equation}

Though I fail to see how to linearize these conditions.

Again thanks a lot for your help (and patience ;) ).

 

[Chestermiller]

Sorry about this late reply, I was lately a bit short on time.
I am a bit confused, to be honest...
Pls, correct me if I am wrong but if I understand everything correctly you expanded my equation (27) with your argument, neglected the term with transient temperature,

I neglected the term involving the product of transient quantities as small compared to zero-order and first-order terms in transient quantities.

\begin{equation} r = 0: \frac{\partial T^1 }{\partial r} + \frac{\partial T^0 }{\partial r} = 0 ; \end{equation}

But, $$\frac{\partial T^0 }{\partial r} = 0$$ Therefore, at r = 0, $$ \frac{\partial T^1 }{\partial r}=0$$

and

\begin{equation} r = R:k \frac{\partial T^1 }{\partial r} + \frac{\partial T^0 }{\partial r} = q"_{wall} ; \end{equation}

But, $$k \frac{\partial T^0 }{\partial r} = q"_{wall}$$ Therefore, at r = R, $$\frac{\partial T^1 }{\partial r} =0$$

 

[rafisondi]

I think I finally understand what you were going for. You are linearizing with respect to the steady-state solution using the steady-state equations as "critical points". So now using your solution

\begin{equation}

\frac{\partial T^1}{\partial t}+w^0\frac{\partial T^1}{\partial x}=-w^1\frac{\partial T^0}{\partial x}+\alpha \frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial T^1}{\partial r}\right)

\end{equation}

And the linearized BC:

\begin{equation}\frac{\partial T^1 }{\partial r}=0\end{equation} at r=0 and
\begin{equation}\frac{\partial T^1 }{\partial r} = 0\end{equation} at r=RAs for the solution of this PDE, I will probably have to assume an oscillating temperature profile as an Ansatz:

\begin{equation}\frac{\partial T^1}{\partial t}=A \cdot e^{i \omega \cdot t}\end{equation}

where the frequency is the same as for the velocity.

 

[Chestermiller]

Not exactly, You do the exact same thing as for the velocity: $$T^1=A(x,r)e^{i\omega t}$$

 

[rafisondi]

Sure no problem. It might take me a little while to derive a reasonable solution but now I finally got something to go on. I'll post as soon as I get a somewhat closed solution. Again thanks a lot for your help.