Coupled heat transfer validation: Heat transfer between concentric cylinders
1. Objective
The objective of this project is to validate a coupled heat transfer case including the new ‘Radiation’ feature in the Convection Heat Transfer analysis of SimScale. A project including conduction, convection and radiation has been solved analytically, and using the platform.
Import this project to your workbench
2. Case description
The geometry consists of two horizontal concentric cylinders with air in between. The inner cylinder (diameter \(D_i\)) has higher temperature (\(T_i\)) than the outer cylinder (diameter \(D_o\)), with temperature (\(T_o\)). Both have the same axial length \(L\) and have an emissivity of \(\varepsilon=0.85\). Side rings are considered adiabatic and black surfaces (\(\varepsilon=1.0\)).
\(D_o={1}{m}\)
\(D_i={0.9}{m}\)
\(L={1}{m}\)
\(T_o={20}{ºC}={293.15}{K}\)
\(T_i={50}{ºC}={323.15}{K}\)
3. Analytical Solution
3.1 Conduction
Conduction is the heat transfer due to a gradient on the particles vibration. Regarding most of the newtonian fluids, it is usually much smaller than convection (this ratio is represented by the Rayleigh Number). For steady-state, conduction is quantified by the Fourier’s (or Newton’s heat) law. Considering a constant conductivity k, the amount of conductive heat \(q_{cond}\) that goes through a surface of area A is:
\(\frac{q_{cond}}{A}=-k\nabla T\)
Taking into account the geometry, the above equation can be expressed in cylindrical coordinates (\(r,\theta,z\)) as follows,
\(\frac{q_{cond}}{A}=-k\left[\frac{\partial T}{\partial r} + \frac{1}{r}\frac{\partial T}{\partial \theta} + \frac{\partial T}{\partial z}\right]\)
Also, considering the axisymmetry of our problem and neglecting temperature gradients in the \(z\) direction, it can be seen that:
\(\frac{q_{cond}}{A}=-k \frac{dT}{dr}\label{eq:prev}\)
Integrating between both cylinders (\(D_o=2r_o and D_i=2r_i\)):
\(\int_{r_i}^{r_o} \frac{dr}{r}= -\frac{2\pi k (T_o-T_i)}{q_{cond}}\)
\(q_{cond}=-\frac{2\pi k (T_o-T_i)}{\ln \frac{r_o}{r_i}}\)
Substituting the values from Appendix A, we have,
\(q_{cond}={46.963}{W}\)
3.2 Convection
Convection is the heat transfer caused by bulk flows when a temperature gradient is applied to a material. Analytical results for convection are cumbersome to carry out due to its transient nature.
The method used here can be found in various references and is based on the conduction equations . Thermal conductivity k is substituted by an ”effective” thermal conductivity \(k_{eff}\). This value is the conductivity of a stationary fluid that will transfer the same amount of heat as the moving fluid:
\(q_{conv}=-\frac{2\pi k_{eff} (T_o-T_i)}{\ln \frac{r_o}{r_i}}\)
The ratio between both conductivities is given by,
\(\frac{k_{eff}}{k}=0.386 \left(\frac{Pr}{0.861+Pr} \right)^{1/4}Ra_c^{1/4}\)
where \(Pr\) is the Prandtl number of the fluid (\(Pr=\mu c_p / k\)) and \(Ra_c\) is a modified version of the Rayleigh number (\(Ra=\frac{g\beta (T_i-T_o)L_c^3}{\nu \alpha}\)) using a characteristic length given by:
\(L_c=\frac{2\left[\ln (r_o/r_i)\right]^{4/3}}{\left(r_i^{-3/5}+r_o^{-3/5}\right)^{5/3}}\)
This method is valid if and only if \(0.7\leq Pr \leq 6000\) and \(Ra_c \leq 10^7\). In addition, the minimum heat transfer rate between the cylinders cannot fall below the conduction limit (\(k_{eff}/k\geq 1\))
Using the air properties specified in Appendix A, we have:
\(q_{conv}={142.586}{W}\)
3.3 Radiation
Radiation is the heat transfer mean caused by the emission of electromagnetic waves. The method used by SimScale is a diffuse model based on view factors. The best way to calculate this heat exchange between surfaces analytically is to build the radiation network of the system.
Here subscript \(i\) represents the inner cylinder, \(o\) refers to the outer cylinder, and \(R\) is the sum of both side rings. Thus, \(q_{i \rightarrow R}\) is the amount of radiative heat that goes from \(i\) to \(R\). Also, it is assumed that \(T_R=(T_o+T_i)/2\), which can be checked with an Area Average result control item in SimScale.
We can calculate the radiosity \(J_i\) by calculating the power balance at the central node:
\(q_{rad}=q_{i \rightarrow R}+q_{i \rightarrow o}\)
\(\frac{E_{bi}-J_i}{\frac{1-\varepsilon_i}{\varepsilon_i A_i}}=\frac{J_i-E_{bR}}{\frac{1}{A_i F_{i,R}}}+\frac{J_i-E_{bo}}{\frac{1}{A_i F_{i,o}}+\frac{1-\varepsilon_o}{\varepsilon_o A_o}}\)
Complexity arises when we calculate the view factors. Analytical solutions for these values are easily calculated for infinite cylinders. However, in our case the cylinders are finite in length. As a result, plots from literature were used to estimate the view factors.
For our case we have, \((r_i/r_o=0.9, L/r_o=2)\). Using figure 3, the view factors for the outer cylinder are:
\(F_{o,i}=0.88\) and \(F_{o,o}=0.08\)
As a result:
\(F_{i,o}=\frac{A_o}{A_i}F_{o,i}=0.9777\)
\(F_{i,R}=1-F_{i,o}=0.0223\)
With above values in hand we calculate the radiosity \(J_i\):
\(J_i={592.079}{W/m^2}\)
Now, we can calculate the net radiative heat that is getting into the domain through the inner cylinder:
\(q_{rad}=\frac{E_{bi}-J_i}{\frac{1-\varepsilon_i}{\varepsilon_i A_i}}={420.111}{W}\)
3.4 Total heat flux
The analytical total heat flux \(q_{net}\) getting into the domain through the inner cylinder can be expressed as:
\(q_{net}=q_{cond}+q_{conv}+q_{rad}={609.66}{W}\)
4. SimScale solution
4.1 Simulation setup
A CAD model of the fluid contained between the cylinders was created. In order to validate the SimScale results against the analytical ones, the following Boundary Conditions were applied:
- \(q_{net}={609.66}{W}\) was assigned as a Turbulent Heat Flux to the inner cylinder.
- Fixed temperature \(T_o={293.15}{K}\) was assigned to the outer cylinder.
- The side rings were assigned a no slip adiabatic wall BC with \(\varepsilon=1.0\).
This set of Boundary Conditions is the closest one to the analytical problem described previously. The objective is to obtain \(T_i={323.15}{K}\) on the surface of the inner cylinder (average value).
5. Result comparison
Results are presented in the table below for different models and setups:
6 Appendix
A. Air properties
For conduction and convection, the properties used for the air must be the ones at the film Temperature \(T_f (T_f=(T_o+T_i)/2={308.15}{K}={35}{ºC}\)). Besides, the value of \(Pr\) has been taken constant for all the \(Re\) regimes (only for validation purposes):
\(ρ={1.145}{kg/m^3}\) \(\mu={1.891e -05 }{kg/ms}\) \(\nu={1.655e -05 }{m^2/s}\) \(k={2.622e -02 }{W/mK}\) \(c_p={1007}{J/kgK}\) \(\alpha={2.277e -05 }{m^2/s}\) \(\beta={0.003245}{1/K}\) \(Pr=0.7268\)B. A deeper look. Heat fluxes
There is another variable that can be used to validate the SimScale results. A surface integral can be made for \(Q_r\) in SimScale along the inner cylinder. This leads to \(q_{rad}=\iint_i Q_rdS\). The results are presented in the table below:
Due to the inevitable agglomeration process, the variation of \(T_i\) along the inner cylinder surface is much more smooth than the \(Q_r\) one.
Presented below are the post-processing results from the SimScale post-processor.
C. Conduction profile validation
By solving the conduction differential equation for a generic radius \(r_i<r<r_o\), we can get the temperature profile \(T(r)\) :
\(\frac{q_{cond}}{A}=-k \frac{dT}{dr}\)
\(\int_{r}^{r_o} \frac{dr}{r}= -\frac{2\pi k (T-T_i)}{q_{cond}}\)
\(T(r)=T_i-\frac{q_{cond}}{2\pi k}\ln\left(\frac{r}{r_i}\right)\)
This is the function that the Temperature profile would have in a case with only conduction. We can set this project up in SimScale by disabling radiation and assigning a zero value to all the gravity components. Using \(q_{cond}\) and the \(k\) specified in Appendix A, results obtained can be seen in the figure below:
With a Root Mean Square Error of \(RMSE={0.145}{K} \)