4.4.2 Cooling rib with variable cross-section
Figure 56: A cooling rib of length L, width b and a varying cross-sectional area \( A(X) \) (trapezoidal).
In Figure 56 trapezoidal cooling rib with length L and width b is illustrated. The thickness varies from d at \( X=0 \) to D at \( X=L \). The cooling rib has a heat loss at \( X=0 \) and a prescribed temperature \( T_L \) at \( X_L \), with a constant surrounding temperature \( T_{\infty} \) .
In Figure 57 the upper half of the cooling rib is illustrated for a better quantitative demonstration of the geometry, where \( d \) and \( D \) denote the diameters of the cooling rib at the inlet and outlet, respectively, whereas \( L \) denotes the length of the rib.
Figure 57: A cooling rib with a trapezoidal cross-section.
The following geometrical relations may be deducted: $$ \begin{equation*} \tan(\phi) = \frac{D/2 - d/2}{L} = \frac{d^*/2-d/2}{X} \end{equation*} $$
Further, the spatial dependent diameter \( d^* \) may be isolated and presented as: $$ \begin{equation} d^*(X) = d + \left(\frac{D-d}{L}\right)\;X \tag{4.75} \end{equation} $$
We assume that the temperature mostly varies in the X-direction such that the derivation in appendix B in Numeriske beregninger is valid. We assume accordingly that D \( \ll \) L and that \( 0 < \frac{d}{D} < 1 \) .
In Numeriske beregninger the quasi one-dimensional differential equation for heat conduction is derived, and for stationary conditions is reduces to equation (B.0.22): $$ \begin{equation} \frac{d}{dX}\left[A(X) \; \frac{d(T-T_{\infty})}{dX}\right] = \frac{\bar{h}P(X) }{k}(T-T_{\infty}) \tag{4.76} \end{equation} $$
with corresponding boundary conditions: $$ \begin{equation} \tag{4.77} \frac{dT(0)}{dX}= \frac{\bar{h}_0}{k}[T(0)-T_{\infty}] \ , \ T(L) = T_L \end{equation} $$
For the trapezoidal cross-section in Fig. (56): $$ \begin{equation*} P(X) = 2(b+d^*) \approx 2b \qquad \text{ for } d^* \ll b \qquad \text{and} A(X) = b\;d^* \end{equation*} $$
For convenience we introduce dimensionless quantities for temperature \( \theta = \dfrac{T-T_{\infty}}{T_L - T_{\infty}} \), length \( x=\dfrac{X}{L} \), and the diameter ratio \( \alpha = \dfrac{d}{D} \), \( 0 < \alpha < 1 \) and the Biot numbers: $$ \begin{equation} \beta^2 = \frac{2\bar{h} \;L^2}{D\;k} \qquad \text{ and } \qquad \beta_0^2 = \dfrac{\bar{h}_0}{k}\;L \tag{4.78} \end{equation} $$
Equipped with these dimensionless quantities we may render equation (4.76) as: $$ \begin{equation} \frac{d}{dx}\left[\{\alpha + (1-\alpha)\;x\}\frac{d\theta(x)}{dx}\right] - \beta^2\theta(x) = 0 \tag{4.79} \end{equation} $$
and corresponding boundary conditions become: $$ \begin{equation} \tag{4.80} \frac{d\theta}{dx}(0)=\beta_0^2 \;\theta(0) \ , \ \theta(1) = 1 \end{equation} $$
As \( d \rightarrow 0 \) the trapezoidal profile approaches a triangular profile and \( \alpha \rightarrow 0 \), and consequently equation (4.79) takes the form: $$ \begin{equation} \tag{4.81} x\frac{d^2\theta}{dx^2}+\frac{d\theta}{dx} -\beta^2\theta(x) = 0 \end{equation} $$
For the left boundary at \( x = 0 \) we get: $$ \begin{equation} \frac{d\theta}{dx}(0)-\beta^2\theta(0) = 0 \tag{4.82} \end{equation} $$
An analytical solution of equation (4.79) and (4.81) is outlined in G.5 of Numeriske Beregninger.
4.4.3 Example: Numerical solution for specific cooling rib
In the current example we will look at the solution for a specific cooling rib based on the generic representation above. We consider the a rib with the following geometry: $$ \begin{equation} L=0.1m, \quad D = 0.01m, \quad d=0.005m,\ \tag{4.83} \end{equation} $$
and physical constants: $$ \begin{equation} \bar{h} = 80W/m^2/^\circ C, \bar{h}_0 = 200 W/m^2/^\circ C, \text{og }k = 40W/m/^\circ C \tag{4.84} \end{equation} $$
Together, these particular values of geometrical and physical variables correspond to \( \beta^2 = 4.0 \), \( \alpha = \frac{1}{2} \) and \( \beta_0^2 = 0.5 \) and equation (4.79) will consequently take the form: $$ \begin{equation} \frac{d}{dx}\left[ (1+x)\;\frac{d\theta}{dx}\right] - 8\;\theta(x) = 0 \tag{4.85} \end{equation} $$
Analytical solution
An analytical solution of the resulting differential equation (4.85) for this example may be found in appendix G.5 in Numeriske beregninger and a corresponding python snippet is shown below along with the resulting plot:
# coding: utf-8
# src-ch3/trapes.py;
import numpy as np
from matplotlib.pyplot import *
from numpy import sqrt as sqrt
from scipy.special import kv as besselk
from scipy.special import iv as besseli
# change some default values to make plots more readable
LNWDT=3; FNT=11
rcParams['lines.linewidth'] = LNWDT; rcParams['font.size'] = FNT
font = {'size' : 16}; rc('font', **font)
def trapes_analytical(h):
N=int(round(1.0/h)) # number of unknowns, assuming the RHS boundary value is known
x=np.arange(N+1)*h
# x = np.linspace(0, 1, N)
z0 = 4*sqrt(2)
z1 = 8
g = sqrt(2)/8
k1z0, i0z1 = besselk(1, z0), besseli(0, z1)
k0z1, i1z0 = besselk(0, z1), besseli(1, z0)
k0z0, i0z0 = besselk(0, z0), besseli(0, z0)
J = k1z0*i0z1 + k0z1*i1z0 + g*(k0z0*i0z1 - k0z1*i0z0)
A = (k1z0 + g*k0z0)/J
B = (i1z0 - g*i0z0)/J
z = sqrt(32.*(1 + x))
theta = A* besseli(0, z) + B* besselk(0, z)
dtheta = 16*(A*besseli(1, z) - B*besselk(1, z))/z
return x, theta, dtheta
if __name__ == '__main__':
legends=[] # empty list to append legends as plots are generated
h=0.2
x, theta, dtheta = trapes_analytical(h)
plot(x,theta,'b',)
legends.append(r'$\theta$')
plot(x,dtheta,'r')
legends.append(r'$\theta^{\prime}$')
## Add the labels
legend(legends,loc='best',frameon=False) # Add the legends
ylabel(r'$\theta$ and $\theta^{\prime}$')
xlabel('x')
show()
close()
Figure 58: The dimensionless temperature \( \theta \) and the spatial derivative \( \theta^{\prime} \) as a function of the spatial variable \( x \).
Numerical solution
We will solve equation (4.85) numerically by central differences, and first rewrite it as: $$ \begin{equation} (1+x)\frac{d^2\theta}{dx^2}+\frac{d\theta}{dx}-8 \; \theta(x)=0 \tag{4.86} \end{equation} $$
with the boundary conditions: $$ \begin{equation} \frac{d\theta}{dx}(0) = \beta_0^2 \, \theta(0) = \frac{\theta(0)}{2}, \qquad \theta(1) = 1 \tag{4.87} \end{equation} $$
The procedure for this example will follow the along the lines of the procedure in section (4.4.1 Example: Heat exchanger with constant cross-section) for mixed boundary conditions.
Figure 59: Discretization of a bar of unit length.
The x-coordinates are discretized by: \( x_i = i\;h \) , \( i =0,1,2,\ldots,N \) where \( h=\frac{1}{N} \) and a central difference approximation of equation (4.86) leads to the following discretized equation: $$ \begin{equation*} (1+x_i) \;\frac{\theta_{i+1} -2\theta_i + \theta_{i-1}}{h^2}+\frac{\theta_{i+1}-\theta_{i-1}}{2h}-8\theta_i = 0 \end{equation*} $$
which after collection of terms may be presented as for generic point \( x_i \) along the unit length bar in Figure 59: $$ \begin{equation} -(1-\gamma_i)\;\theta_{i-1} + 2 \;(1+8h \;\gamma_i) \;\theta_i - (1+\gamma_i)\;\theta_{i+1}=0 \tag{4.88} \end{equation} $$
where $$ \begin{equation} \gamma_i = \frac{h}{2\;(1+x_i)}=\frac{h}{2\;(1+i\;h)}, \qquad i=1,\ldots,N-2 \tag{4.89} \end{equation} $$
We observe that equation (4.88) form a tridiagonal equation system as only the immediate neighbours of \( \theta_{i} \), namely \( \theta_{i-1} \) and \( \theta_{i-1} \), are involved for the \( i-th \) equation.
At the left boundary we need equation (4.87) to be satisfied, which has the second order discretized equivalent \( \frac{\theta_2 - \theta_0 }{2h} = \frac{\theta_1}{2} \) which yields \( \theta_0 = \theta_2 - \theta_1 \;h \) and after substitution in equation (4.88) for \( i = 0 \) we get the following discretized version at the boundary: $$ \begin{equation} [2+h \;(1+15\;\gamma_0)] \;\theta_0 - 2\; \theta_1 = 0 \tag{4.90} \end{equation} $$
Note that the discretization \( \theta_0 = \theta_2 - \theta_1 \;h \) is equally valid, it is just that the implementation of the boundary condition would then violate the nice tridiagonal structure of the discretized equation system, and thus we prefer equation (4.90) which fits very nicely in the tri-diagonal structure.
For \( i=N-1 \) we make use of the fact that \( \theta_N=1 \) is given : $$ \begin{equation} -(1-\gamma_{N-1}) \;\theta_{N-2} + 2 \;(1 + 8h \;\gamma_{N-1})\;\theta_{N-1} = (1 + \gamma_{N-1}) \; \theta_N= 1 + \gamma_{N-1} \tag{4.91} \end{equation} $$
Together, the equations (4.88), (4.90),and (4.91) form a tridiagonal, linear equation system which we will
solve using the sparse library scipy.sparse from SciPy. An example of how to use the sparse library may
be found in our Introduction to Scientific Python programming.
The python code for the implementation of the algorithm is shown below. We compare the solution with the analytical solution and again we observe the striking difference between the very complex analytical solution (involving Bessel-functions of various kinds) and the relatively simple numerical solution involving only the solution of a tridiagonal linear equation system.
As we have discretized equation (4.86) using central differences, we expect the numerical solution to be second order accurate and have added a code segment to approximate the order of the scheme. You may test the code an see whether our implementation has the expected order. Optionally, we have also provided the possibility to plot the solution as a function of the mesh size.
# coding: utf-8
# src-ch3/trapes_numerical.py;
import numpy as np
from matplotlib.pyplot import *
import scipy as sc
import scipy.sparse
import scipy.sparse.linalg
import time
from trapes import trapes_analytical
# change some default values to make plots more readable
LNWDT=3; FNT=11
rcParams['lines.linewidth'] = LNWDT; rcParams['font.size'] = FNT
font = {'size' : 16}; rc('font', **font)
def solve_trapes_numerically(h):
L=1.0 # length of domain
n=int(round(L/h)) # number of unknowns, assuming the RHS boundary value is known
x=np.arange(n+1)*h # x includes xmin and xmax at boundaries
gamma=h/(2*(1+x[0:-1])) # uses all but the last x-value for gamma
subdiag=gamma[:-1]-1.0 #the subdiagonal has one less element than the diagonal
diag=2.0*(1+8.0*h*gamma)
diag[0]= 2.0+h*(1+15.0*gamma[0])
superdiag=np.zeros(n-1)
superdiag[0]=-2
superdiag[1:]=-(gamma[0:-2]+1.0)
A_diags=[subdiag,diag,superdiag] #note the offdiags are one elt shorter than the diag
A = sc.sparse.diags(A_diags, [-1,0,1], format='csc')
#Create rhs
d=np.zeros(n)
d[n-1]=1+gamma[-1]
#Solve linear problem
tic=time.time()
theta = sc.sparse.linalg.spsolve(A,d) #theta=sc.linalg.solve_triangular(A,d)
toc=time.time()
print('sparse solver time:',toc-tic)
return x,theta
if __name__ == '__main__':
h=0.2
xa, atheta, adtheta = trapes_analytical(h)
x, theta=solve_trapes_numerically(h)
# Plot solutions
figure('Solutions')
plot(xa,atheta,)
plot(x[:-1],theta,'-.')
legends=[] # empty list to append legends as plots are generated
legends.append(r'$\theta$-analytical')
legends.append(r'$\theta$')
## Add the labels
legend(legends,loc='best',frameon=False) # Add the legends
ylabel(r'$\theta$')
xlabel('x')
# Approximate the order of the scheme
h=0.2
h_range=[h/2.**(d) for d in np.arange(0,6)]
error=[]
for h in h_range:
xa, atheta, adtheta = trapes_analytical(h)
x,theta=solve_trapes_numerically(h)
error.append(np.sqrt(np.sum(((atheta[:-1]-theta)/atheta[:-1])**2)))
log_error = np.log2(error)
dlog_err=log_error[0:-1]-log_error[1:]
print(('Approximate order:{0:.3f}'.format(1./np.mean(dlog_err))))
plot_error=0
if (plot_error):
figure('error')
ax = plot(h_range,error)
yscale('log',basey=2)
xlabel(r'meshsize $h$')
ylabel(r'$\frac{\epsilon}{\theta_a}$')
show()
close()