Numerical solution of internal convection

From Thermal-FluidsPedia

Jump to: navigation, search
Typical boundary conditions for internal convection
Typical boundary conditions for internal convection.

The algorithms to solve the convection-diffusion equation and flow field that were addressed in computational methodologies for forced convection are still applicable to the internal convection problems. However, special attention must be paid to proper treatment of the boundary conditions. For an internal flow and heat transfer problem, there are several different types of boundary conditions (see figure to the right): (1) inflow condition, (2) axisymmetric condition, (3) impermeable solid surface, and (4) outflow condition.[1]

The inflow conditions are normally specified by a given distribution of the inlet velocity and the general variable, φ, at the inlet (x = 0). The axisymmetric condition can be implemented by setting the gradient of the velocity component in the y-direction, \partial u/\partial y, and the gradient of the general variable, \partial \phi /\partial y, equal to zero along the axisymmetric line (y = 0). In addition, the velocity component in the y-direction, v, at the axisymmetric boundary should also be zero. For the impermeable solid surface, all three kinds of boundary conditions for the general variables are possible: specified φ (the first kind), specified gradient of φ in the direction perpendicular to the impermeable surface, \partial \phi /\partial y, (the second kind), or specified relation between φ and \partial \phi /\partial y (the third kind). The velocity components in both the tangential and normal direction of the impermeable surface are zero, i.e., u = v = 0, due to no-slip and impermeable conditions. No special treatment for the boundary conditions at inflow, axisymmetric, and impermeable surfaces is required. On the contrary, the outflow boundary condition requires special treatment as outlined below.

Since the momentum equation in the internal flow direction (the x-direction in above figure) and the conservation equation for the general variable are elliptic, the second order derivatives with respect to x appeared in the partial difference equation. Mathematically, boundary conditions at both the inflow and outflow boundary need to be specified. While the inflow boundary conditions are always known, the outflow boundary conditions are usually unknown. Unless the experimentally measured distribution at the outflow boundary is available, we cannot directly use the algorithms introduced in the previous chapter to solve the internal convection problem.

A coordinate (such as x in above figure, or time for a transient problem) can be either two-way or one-way depending on the nature of the problem [2][3]. If the condition at a given location, x, is influenced by changes of conditions at either side of the given location, the coordinate is said to be two-way. On the other hand, if the conditions at the given location, x, are influenced by changes of conditions from only one side, the coordinate is said to be one-way. Mathematically, if the second order derivation with respect to x appeared in the partial differential equation, the coordinate in the x-direction is two-way. For example, the spatial coordinate in the one-dimensional steady-state heat conduction in a fin is two-way because the temperature at any point is influenced by the temperatures at either side, and the second order derivative appears in the energy equation that describes heat conduction in the fin. On the other hand, time is a one way coordinate because the conditions at a given time are only influenced by what happened before that time, and what happen afterward will not affect the conditions at the current time. While the space coordinate in heat and mass transfer is normally always two-way, it can become one-way for some special cases. For internal flow, the change of condition at a given point will have a more profound effect on the conditions of the points downstream, while its influence on the conditions of the points upstream will be relatively weak. For the case that convection overpowers diffusion, one can assume that the condition changes at one point can only propagate downstream and the coordinate becomes one-way.

Control volume for continuity equation at boundary
Control volume for continuity equation at boundary.

The most commonly used approach to handle the outflow boundary condition is to assume the coordinate at the outflow boundary is locally one-way. Thus, the value of the general variable, φ, at the inner point, P, shown in the figure to the right is not affected by the value of φ at the outflow boundary grid point, E, which is located downstream. Computationally, one can set aE = 0 in the discretized equation for the point P. This simple treatment effectively avoids the problem associated with an unknown value ofφat the outflow boundary. It follows from computational methodologies for forced convection that this treatment will be more accurate for cases with a high Peclet number.

Another situation that can be handled by a similar approach is the case that the unknown variable is fully developed at the outflow boundary, in which case the outflow boundary condition becomes \partial \phi /\partial x=0. The implication of the fully developed condition at the outflow is exactly the same as the local one-way behavior discussed above, and can be treated using the same approach. However, it should be pointed out that for the case of fully developed heat transfer, one must not use \partial T/\partial x=0, and the correct condition for fully developed heat transfer should be \partial [(T-{{T}_{m}})/({{T}_{w}}-{{T}_{m}})]/\partial x=0 (see fully-developed flow and heat transfer).

Many practical problems are not similar to the cases corresponding to simple internal channel flow. In many practical applications for internal flow, one needs to deal with conjugate effects, compressibility, multi-domain, and porous media, as well as non-conventional boundary conditions.

Heat pipe configuration
Heat pipe configuration.

To show the power of numerical simulation, the modeling of a conventional wicked heat pipe with variable heat flux (see figure to the right) will be presented here, compared to the non-conventional internal flow case presented in previous sections. The problem is modeled as a two-dimensional conjugate, with compressible flow including the effect of porous media and coupling between various regions.

Since the heat pipe in the figure is closed at both ends, it is required that the vapor which flows out of the evaporator segment enters into the condenser section. The conservation of mass, momentum, and energy equations for the compressible vapor flow region including viscous dissipation are

\frac{1}{r}\frac{\partial }{\partial r}({{\rho }_{v}}r{{v}_{v}})+\frac{\partial }{\partial z}({{\rho }_{v}}{{w}_{v}})=0


\begin{align}  & {{\rho }_{v}}\left( {{w}_{v}}\frac{\partial {{w}_{v}}}{\partial z}+{{v}_{v}}\frac{\partial {{w}_{v}}}{\partial r} \right)=-\frac{\partial {{p}_{v}}}{\partial z} \\  & +{{\mu }_{v}}\left[ \frac{4}{3}\frac{{{\partial }^{2}}{{w}_{v}}}{\partial {{z}^{2}}}+\frac{1}{r}\frac{\partial }{\partial r}\left( r\frac{\partial {{w}_{v}}}{\partial r} \right)+\frac{1}{r}\frac{\partial }{\partial r}\left( r\frac{\partial {{v}_{v}}}{\partial z} \right)-\frac{2}{3}\frac{\partial }{\partial z}\left( \frac{1}{r}\frac{\partial }{\partial r}(r{{v}_{v}}) \right) \right] \\ \end{align}


\begin{align}  & {{\rho }_{v}}\left( {{w}_{v}}\frac{\partial {{v}_{v}}}{\partial z}+{{v}_{v}}\frac{\partial {{v}_{v}}}{\partial r} \right) \\  & =-\frac{\partial {{p}_{v}}}{\partial r}+{{\mu }_{v}}\left[ \frac{{{\partial }^{2}}{{v}_{v}}}{\partial {{z}^{2}}}+\frac{4}{3r}\frac{\partial }{\partial r}\left( r\frac{\partial {{v}_{v}}}{\partial r} \right)-\frac{4}{3}\frac{{{v}_{v}}}{{{r}^{2}}}+\frac{1}{3}\frac{{{\partial }^{2}}{{w}_{v}}}{\partial z\partial r} \right] \\ \end{align}


\begin{align}  & {{\rho }_{v}}{{c}_{pv}}\left( {{w}_{v}}\frac{\partial {{T}_{v}}}{\partial z}+{{v}_{v}}\frac{\partial {{T}_{v}}}{\partial r} \right) \\  & =\frac{{{k}_{v}}}{r}\left[ \frac{\partial }{\partial r}\left( r\frac{\partial {{T}_{v}}}{\partial r} \right)+r\frac{{{\partial }^{2}}{{T}_{v}}}{\partial {{z}^{2}}} \right]+{{v}_{v}}\frac{\partial {{p}_{v}}}{\partial r}+{{w}_{v}}\frac{\partial {{p}_{v}}}{\partial z}+{{\mu }_{v}}\Phi  \\ \end{align}


where the viscous dissipation is

\Phi =2\left[ {{\left( \frac{\partial {{v}_{v}}}{\partial r} \right)}^{2}}+{{\left( \frac{{{v}_{v}}}{r} \right)}^{2}}+{{\left( \frac{\partial {{w}_{v}}}{\partial z} \right)}^{2}}+\frac{1}{2}{{\left( \frac{\partial {{v}_{v}}}{\partial z}+\frac{\partial {{w}_{v}}}{\partial r} \right)}^{2}}-\frac{1}{3}{{(\nabla \cdot {{\mathbf{V}}_{v}})}^{2}} \right]


and \nabla \cdot {{\text{V}}_{v}} is given by

\nabla \cdot {{\mathbf{V}}_{v}}=\frac{1}{r}\frac{\partial }{\partial r}(r{{v}_{v}})+\frac{\partial {{w}_{v}}}{\partial z}


The ideal gas law (p = ρvRgTv) is employed to account for the compressibility of the vapor. The use of liquid capillary action is a unique feature of the heat pipe. From a fundamental point of view, the liquid capillary flow in heat pipes with screen wicks should be modeled as a flow through a porous media. It is assumed that the wicking material is isotropic and of constant thickness. In addition, the wick is saturated with liquid, and the vapor condenses and the liquid evaporates at the liquid-vapor interface. The averaging technique has been applied by many investigators to obtain the general equation which describes the conservation of momentum in a porous structure. Since the development of Darcy’s semi-empirical relation, which characterizes the fluid motion under certain conditions, many researchers have tried to develop and extend Darcy’s law in order to see the effect of the inertia terms. In this respect, those who have tried to model the flow with the Navier-Stokes equations were the most successful. The general equations of continuity, momentum and energy for steady state laminar incompressible liquid flow in porous media in terms of the volume-averaged velocities are:

\frac{\partial }{\partial z}({{\rho }_{\ell }}{{w}_{\ell }})+\frac{1}{r}\frac{\partial }{\partial r}({{\rho }_{\ell }}r{{v}_{\ell }})=0


\frac{1}{{{\varepsilon }^{2}}}\left( {{w}_{\ell }}\frac{\partial {{w}_{\ell }}}{\partial z}+{{v}_{\ell }}\frac{\partial {{w}_{\ell }}}{\partial r} \right)=-\frac{1}{{{\rho }_{\ell }}}\frac{\partial {{p}_{\ell }}}{\partial z}-\frac{{{\nu }_{\ell }}{{w}_{\ell }}}{K}+\frac{{{\nu }_{\ell }}}{\varepsilon }\left[ \frac{1}{r}\frac{\partial }{\partial r}\left( r\frac{\partial {{w}_{\ell }}}{\partial z} \right)+\frac{{{\partial }^{2}}{{w}_{v}}}{\partial {{z}^{2}}} \right]


\frac{1}{{{\varepsilon }^{2}}}\left( {{w}_{\ell }}\frac{\partial {{v}_{\ell }}}{\partial z}+{{v}_{\ell }}\frac{\partial {{v}_{\ell }}}{\partial r} \right)=-\frac{1}{{{\rho }_{\ell }}}\frac{\partial {{p}_{\ell }}}{\partial r}-\frac{{{\nu }_{\ell }}{{v}_{\ell }}}{K}+\frac{{{\nu }_{\ell }}}{\varepsilon }\left[ \frac{1}{r}\frac{\partial }{\partial r}\left( r\frac{\partial {{v}_{\ell }}}{\partial z} \right)-\frac{{{v}_{\ell }}}{{{r}^{2}}}+\frac{{{\partial }^{2}}{{v}_{v}}}{\partial {{z}^{2}}} \right]


{{\rho }_{\ell }}{{c}_{p\ell }}\left( {{w}_{\ell }}\frac{\partial {{T}_{\ell }}}{\partial z}+{{v}_{\ell }}\frac{\partial {{T}_{\ell }}}{\partial r} \right)=\frac{1}{r}\frac{\partial }{\partial r}\left( r{{k}_{eff}}\frac{\partial {{T}_{\ell }}}{\partial r} \right)+\frac{\partial }{\partial z}\left( {{k}_{eff}}\frac{\partial {{T}_{\ell }}}{\partial z} \right)


where ε is the volume fraction or porosity of the wick, and K is the permeability of the wick structure.

The effective thermal conductivity of the wick, keff, is related to the thermal conductivity of the solid and liquid phases

{{k}_{eff}}=\frac{{{k}_{\ell }}[({{k}_{\ell }}+{{k}_{s}})-(1-\varepsilon )({{k}_{\ell }}-{{k}_{s}})]}{[({{k}_{\ell }}+{{k}_{s}})+(1-\varepsilon )({{k}_{\ell }}-{{k}_{s}})]}


The steady state energy equation that describes the temperature in the heat pipe wall is

\frac{\partial }{\partial z}\left( {{k}_{w}}\frac{\partial {{T}_{w}}}{\partial z} \right)+\frac{1}{r}\frac{\partial }{\partial r}\left( {{k}_{w}}r\frac{\partial {{T}_{w}}}{\partial r} \right)=0


where kw is the local thermal conductivity of the heat pipe wall. The boundary conditions are

\begin{matrix}   {} & {}  \\\end{matrix}w(0,r)=v(0,r)=0


\begin{matrix}   {} & {}  \\\end{matrix}w({{L}_{t}},r)=v({{L}_{t}},r)=0


\begin{matrix}   {} & {}  \\\end{matrix}w(z,{{R}_{v}})=0


\begin{matrix}   {} & {}  \\\end{matrix}w(z,{{R}_{i}})=0


\frac{\partial w}{\partial r}(z,0)=\frac{\partial v}{\partial r}(z,0)=0


v(z,{{R}_{v}})=\frac{1}{{{\rho }_{v}}{{h}_{\ell v}}}\left( {{\left. -{{k}_{eff}}\frac{\partial {{T}_{\ell }}}{\partial r} \right|}_{r={{R}_{v}}}}+{{\left. {{k}_{v}}\frac{\partial {{T}_{v}}}{\partial r} \right|}_{r={{R}_{v}}}} \right)


v(z,Ri) = 0


\frac{\partial T}{\partial z}(0,r)=\frac{\partial T}{\partial z}({{L}_{t}},r)=0


T(z,{{R}_{v}})={{\left[ \frac{1}{{{T}_{0}}}-\frac{{{R}_{g}}}{{{h}_{\ell v}}}\ln \left( \frac{{{p}_{v}}}{{{p}_{0}}} \right) \right]}^{-1}}


{{k}_{w}}\frac{\partial {{T}_{w}}}{\partial r}(z,{{R}_{0}})=\pm \frac{Q(z)}{A}


In the boundary conditions specified by eqs. (13)-(22), the heat flux (Q(z)/A) was uniform and positive at the active evaporators, and zero in all adiabatic and transport sections as well as in the inactive evaporators. In the condenser section, the heat flux was assumed to be uniform and negative based on the total heat input through the evaporators. The temperature along the liquid-vapor interface was taken as the equilibrium saturation temperature corresponding to its equilibrium pressure condition. Through the above procedure, both the effects of energy conduction in the heat pipe wall and the liquid flow in the wick are taken into account.

For simplicity and generality, the problem should be solved as a single domain problem using conjugate heat transfer analysis. To achieve this, the conservation equations for mass, momentum, and energy are generalized such that the energy equations with temperature as a dependent variable in the three regions (i.e., wall, wick, and vapor) have the same source term. In addition, the continuity of temperature, heat flux, and mass flux, as well as some special boundary conditions such as thermodynamic equilibrium, should be satisfied at each of the interfaces. The energy equation can be written in terms of temperature as

\rho \frac{DT}{Dt}=\frac{k}{{{c}_{p}}}{{\nabla }^{2}}T+\frac{s}{{{c}_{p}}}


where s is the source term which includes viscous dissipation and pressure work. It should be noted that solving the problem in terms of enthalpy does not preserve the condition of temperature continuity at an interface when harmonic averaging is employed, and will lead to an incorrect calculation of diffusion. For the vapor region, the energy equation is

{{\rho }_{v}}\frac{DT}{Dt}=\frac{{{k}_{v}}}{{{c}_{p,v}}}{{\nabla }^{2}}T+\frac{s}{{{c}_{p,v}}}


Multiplying both sides of the energy equation for the liquid-wick region by {{c}_{p,\ell }}/{{c}_{p,v}} results in

{{\rho }_{l}}\frac{{{c}_{p,l}}}{{{c}_{p,v}}}\frac{DT}{Dt}=\frac{{{k}_{l}}}{{{c}_{p,v}}}{{\nabla }^{2}}T+\frac{s}{{{c}_{p,v}}}


Similarly, the energy equation for the heat pipe wall is

{{\rho }_{w}}\frac{{{c}_{p,w}}}{{{c}_{p,v}}}\frac{DT}{Dt}=\frac{{{k}_{w}}}{{{c}_{p,v}}}{{\nabla }^{2}}T+\frac{s}{{{c}_{p,v}}}


From this transformation, the following can be observed:
1. The source term for the energy equation is divided by cp,v for all three regions.
2. The density is equal to the modified density, ρ*, for different regions, i.e., ρv for vapor, {{\rho }_{\ell }}{{c}_{p,\ell }}/{{c}_{p,v}}for liquid, and ρwcp,w/cp,v for solid.
3. The diffusion coefficients contain only one cp, namely, cp,v. The transformation makes possible an exact representation of diffusion across an interface when the harmonic average is used for the diffusion coefficient at the interfaces.
The momentum equation for the vapor flow does not need any special treatment since ρ*= ρv in the vapor region. In the wick region, the momentum equation can be transformed in the same manner:

\rho *\frac{DV}{Dt}=-\nabla p*+{{v}_{l}}\rho *V-\frac{{{v}_{l}}\rho *\varepsilon }{K}V


It should be noted that:
1. A source term was added in the momentum equation for the porosity effect in the wick region.
2. The pressure solved for is the modified pressure, p*, which is proportional to the actual pressure. The actual pressure drop between two points in the flow fluid can be calculated as

\begin{matrix}   {} & {}  \\\end{matrix}\Delta p={{c}_{p,v}}/{{c}_{p,l}}\Delta p*


The transformed continuity equation can be written as

\nabla \cdot (\rho V)=0


The finite volume method discussed above was employed by Faghri and Buchko [4] in the solution of the elliptic conservation eqs. (1)-(4) and (12) subject to the boundary conditions (13)-(22). The source terms due to viscous dissipation and pressure work in the energy equation, and the source term due to the porous matrix in the momentum equation are linearized, and the SIMPLEST algorithm was employed for the momentum equations. The solution procedure is based on a line-by-line iteration method in the axial direction and the Jacobi point-by-point procedure in the radial direction.

The energy equation is not continuous at the liquid-vapor interface due to the latent heat of evaporation and condensation. The term {{h}_{\ell v}}must be added as a source term in the energy equation at the liquid-vapor interface. For the first few iterations, it was assumed the heat flux at the liquid-vapor interface is equal to the outer wall heat flux. An exact energy balance given by eq. (18) is satisfied after these initial iterations.


  1. Faghri, A., Zhang, Y., and Howell, J. R., 2010, Advanced Heat and Mass Transfer, Global Digital Press, Columbia, MO.
  2. Patankar, S.V., 1980, Numerical Heat Transfer and Fluid Flow, Hemisphere, Washington, DC.
  3. Patankar, S.V., 1991, Computation of Conduction and Duct Flow Heat Transfer, Innovative Research.
  4. Faghri, A., and Buchko, B., 1991, “Experimental and Numerical Analysis of Low Temperature Heat Pipes with Multiple Heat Sources”, ASME J. Heat Transfer, Vol. 113, pp. 728-734.