• «
• 1
• ...
• 34
• 35
• 36
• 37
• 38
• 39
• 40
• 41
• 42
• 43
• 44
• 45
• 46
• 47
• 48
• 49
• 50
• 51
• ...
• 79
• »

4.4 Examples

4.4.1 Example: Heat exchanger with constant cross-section

A circular cylindrical rod of length L, environmental temperature \( T_{\infty} \), and a constant temperature \( T_r \) at $ X = L$ is illustrated in Figure 53. The physical parameters relevant for the heat conduction are the heat transfer coefficient \( \bar{h} = 100 W/m^{2 \ \circ}C \) and the thermal conductivity \( k = 200 W/m/^\circ C \).

From equation (B.0.22), appendix B in Numeriske Beregninger we find that the governing equation for the heat transfer in the circular cylindrical rod in Figure 53 is governed by: $$ \begin{equation} \tag{4.44} \frac{d}{dX} \left[A(X) \;\frac{d(T - T_{\infty}}{dX} \right] = \dfrac{\bar{h}P}{k}(T-T_\infty) \end{equation} $$

which for rod with constant cross-section, i.e. \( A(x) = A \) reduces to: $$ \begin{equation} \frac{d^2}{dX^2}(T - T_{\infty}) = \frac{\bar{h}P}{kA}(T- T_{\infty}) \tag{4.45} \end{equation} $$

For convenience we introduce dimensionless variables: $$ \begin{equation} \tag{4.46} x = \frac{X}{L} \ ,\ \theta = \frac{T - T_{\infty}}{T_r - T_{\infty}}\ ,\ \beta^2 = \frac{\bar{h}P}{kA}L^2 \end{equation} $$

where L is a characteristic length (see Figure B.5 i Numeriske Beregninger). A relative, dimensionless temperature may be constructed by using a reference temperature \( T_r \) and temperature and the temperature \( T_\infty \) in combination. The parameter \( \beta^2 \) is frequently referred to as the Biot-number, \( Bi=\beta^2 \).

Given the assumptions and prerequisites above (4.44) may be presented in the more convenient and generic form: $$ \begin{equation} \tag{4.47} \frac{d^2 \theta}{dx^2}-\beta^2\theta = 0 \end{equation} $$

which is a second order ODE with the following generic analytical solution: $$ \begin{equation} \theta(x) = A\sinh(\beta x)= - B\cosh(\beta x) \tag{4.48} \end{equation} $$

where the constants \( A \) and \( B \) must be determined from the boundary conditions.

In the following we will investigate the numerical solution to two different boundary conditions commonly denoted as prescribed boundary conditions (which for this particular corresponds to prescribed temperatures) and mixed boundary conditions.

Prescribed temperatures: $$ \begin{align*} T &= T_{\infty} & \text{ for } X = 0 \\ T & = T_ r & \text{ for } X = L \end{align*} $$

which may be presented on dimensionless form as: $$ \begin{align} \theta & = 0 & \text{ for } x = 0 \tag{4.49} \\ \theta &= 1 & \text{ for } x = 1 \nonumber \end{align} $$

The analytical solution (4.48) for these particular boundary conditions reduces to: $$ \begin{equation} \theta(x) = \frac{\sinh(\beta x)}{\sinh(\beta)} \ , \ \frac{d \theta}{dx} = \beta \;\frac{\cosh (\beta x)}{\sinh(\beta)} \tag{4.50} \end{equation} $$

With \( D = 0.02m \) and \( L = 0.2m \) the Biot-number \( \beta^2 = 4 \). By doubling the length the Biot-number quadruples. $$ \begin{equation} \beta^2 = \frac{\bar{h}P}{kA} \;L^2 = \frac{2}{D} \;L^2 \tag{4.51} \end{equation} $$

Mixed boundary conditions: $$ \begin{align*} Q_x &= 0 = \frac{dT}{dX} & \text{ for } X = 0 \\ T & = T_r & \text{ for } X = L \end{align*} $$

A zero temperature gradient at \( X = 0 \) corresponds to an isolated rod at that location. The dimensionless representation of this boundary conditions is: $$ \begin{align} \frac{d \theta}{dx} & =0 & \text{ for } x = 0 \tag{4.52} \\ \theta &= 1 & \text{ for } x = 1 \tag{4.53} \end{align} $$

The analytical solution (4.48) reduces for the mixed boundary conditions to: $$ \begin{equation} \theta(x) = \frac{\cosh(\beta x)}{\cosh(\beta)} \ , \ \frac{d \theta}{dx} = \beta \;\frac{\sinh(\beta x)}{\cosh(\beta)} \tag{4.54} \end{equation} $$

4.4.1.1 Numerical solution

We will use central differences for the term \( \dfrac{d^2 \theta}{dx^2} \) and with \( \dfrac{d^2 \theta}{dx^2}\bigg|_i \approx \dfrac{\theta_{i-1} - 2\theta_i + \theta_{i+1}}{h^2} \) we get the following difference equation: $$ \begin{equation} \theta_{i-1}-(2+\beta^2 h^2 ) \;\theta_i + \theta_{i+1} = 0 \tag{4.55} \end{equation} $$

Prescribed temperatures

We enumerate the nodes for the unknowns as shown in Figure 54:

The x-coordinates can be denoted in a compact manner by:: \( x_i = i \;h \) , \( i = 0,1,...,N+1 \) where \( h = \frac{1}{N+1} \), i.e. the prescribed temperatures at the boundaries are denoted \( \theta_0 \) and \( \theta_{N+1} \), respectively.

A good practice is to apply the generic scheme (4.55) for nodes in the immediate vicinity of the boundaries, and first we take a look at \( i=1 \) $$ \begin{equation*} \theta_0 - (2 + \beta^2h^2) \;\theta_1 + \theta_2 = 0 \to (2+\beta^2h^2) \; \theta_1 + \theta_2 = 0 \end{equation*} $$

which may be simplified by substitution of the prescribed boundary value \( \theta(0) = \theta_0 = 0 \) $$ \begin{equation*} (2+\beta^2h^2) \; \theta_1 + \theta_2 = 0 \tag{4.56} \end{equation*} $$

For the other boundary at \( i = N \) the generic scheme (4.55) is: $$ \begin{equation*} \theta_{N-1} - (2 + \beta^2h^2) \; \theta_N + \theta_{N+1} = 0 \end{equation*} $$

which by substitution of the prescribed value for \( \theta_{N+1} = 1 \) yields: $$ \begin{equation*} \theta_{N-1} - (2 + \beta^2h^2) \; \theta_N = -1 \tag{4.57} \end{equation*} $$

A complete system of equations may finally be obtained from (4.56), (4.55), and (4.57): $$ \begin{align} i = 1:& \ \ -(2 + \beta^2h^2)\;\theta_1 + \theta_2 = 0 \nonumber \\ i = 2,3,...,N-1:& \ \ \theta_{i-1} - (2 + \beta^2h^2)\;\theta_i + \theta_{i+1} = 0 \tag{4.58}\\ i = N:& \ \ \theta_{N-1} - (2 + \beta^2h^2)\;\theta_N = -1 \nonumber \end{align} $$

See appendix A, section A.5, example A.15 in Numeriske Beregninger, this example is treated in more detail. The following system of coefficients may be obtained by comparison with (4.36): $$ \begin{align} a_i &= 1 \ , \ & i = 2,3,...,N \nonumber \\ b_i &= -(2+ \beta^2h^2) \ , &i = 1,2,...,N \tag{4.59} \\ c_i &= 1 \ , \ &i = 1,2,...,N-1 \nonumber \\ d_i &= 0 \ , \ & i = 1,2,...,N-1 \nonumber \\ d_N& = -1 \nonumber \end{align} $$

In the program ribbe1.py below the linear, tri-diagonal equation system (4.59) resulting from the discretization in (4.55) is solved by using two SciPy modules, the generic scipy.linalg.solve and the computationally more efficient scipy.sparse.linalg.spsolve.

SciPy is built using the optimized ATLAS LAPACK and BLAS libraries and has very fast linear algebra capabilities. Allegedly, all the raw lapack and blas libraries are available for even greater speed. However, the sparse linear algebra module of scipy offers easy-to-use python interfaces to these routines. Numpy has also a linalg-module, however, an advantage of using scipy.linalg over numpy.linalg is that it is always compiled with BLAS/LAPACK support, while this is optional for numpy. Note that the scipy.linalg algorithms do not exploit the sparse nature of the matrix as they use direct solvers for dense matrices.

However, SciPy offers a sparse matrix package scipy.sparse. The spdiags function may be used to construct a sparse matrix from diagonals. Note that all the diagonals must have the same length as the dimension of their sparse matrix - consequently some elements of the diagonals are not used. The first \( k \) elements are not used of the \( k \) super-diagonal, whereas the last \( k \) elements are not used of the \( -k \) sub-diagonal. For a quick tutorial of the usage of these sparse solvers see SciPy sparse examples.

We have implemented a simple means to compare the two solution procedures with a tic-toc statements. You may experiment with various problems sizes (by varying the element size h) to assess the impact on the computational speed of the two procedures. The analytical solution is given by (4.50).

# src-ch3/section321/ribbe1.py

import numpy as np
import scipy as sc
import scipy.linalg
import scipy.sparse
import scipy.sparse.linalg
import time
from math import sinh

#import matplotlib.pyplot as plt
from matplotlib.pyplot import *
# Change some default values to make plots more readable on the screen
LNWDT=3; FNT=20
matplotlib.rcParams['lines.linewidth'] = LNWDT; matplotlib.rcParams['font.size'] = FNT

# Set simulation parameters
beta = 5.0
h = 0.001               # element size
L =1.0                  # length of domain
n = int(round(L/h)) -1  # number of unknowns, assuming known boundary values
x=np.arange(n+2)*h      # x includes min and max at boundaries were bc are imposed.


#Define useful functions

def tri_diag_setup(a, b, c, k1=-1, k2=0, k3=1):
    return np.diag(a, k1) + np.diag(b, k2) + np.diag(c, k3)

def theta_analytical(beta,x):
    return np.sinh(beta*x)/np.sinh(beta)

#Create matrix for linalg solver
a=np.ones(n-1)
b=-np.ones(n)*(2+(beta*h)**2)
c=a
A=tri_diag_setup(a,b,c)

#Create matrix for sparse solver
diagonals=np.zeros((3,n))
diagonals[0,:]= 1                       #all elts in first row is set to 1
diagonals[1,:]= -(2+(beta*h)**2)  
diagonals[2,:]= 1 
A_sparse = sc.sparse.spdiags(diagonals, [-1,0,1], n, n,format='csc') #sparse matrix instance

#Crete rhs array
d=np.zeros(n)
d[n-1]=-1

#Solve linear problems
tic=time.time()
theta = sc.sparse.linalg.spsolve(A_sparse,d) #theta=sc.linalg.solve_triangular(A,d)
toc=time.time()
print('sparse solver time:',toc-tic)

tic=time.time()
theta2=sc.linalg.solve(A,d,)
toc=time.time()
print('linalg solver time:',toc-tic)

# Plot solutions
plot(x[1:-1],theta,x[1:-1],theta2,'-.',x,theta_analytical(beta,x),':')
legend(['sparse','linalg','analytical'])
show()
close()
print('done')

The numerical predictions are presented and compared with analytical values in the table below:

Mixed boundary conditions. Version 1

For mixed boundary conditions the nodes with unknown temperatures are enumerated as shown in Fig. (55), to make the first unknown temperature we compute have index 1. The x-coordinates are given by: \( x_i = (i -1) \;h, \ i = 1,2,...,N+1 \) where \( h = \frac{1}{N} \) and we will use second order central differences as an approximation of the zero-gradient boundary condition \( \frac{d\theta}{dx}=0 \) at the left boundary where \( x = 0 \).

For a generic node 'i' the central difference approximation may be denoted: $$ \begin{equation} \frac{d\theta}{dx}\bigg|_i \approx \frac{\theta_{i+1}-\theta_{i-1}}{2h} \tag{4.60} \end{equation} $$

which for node 1 takes the form: $$ \begin{align} \frac{\theta_2 - \theta_0}{2h}= 0 \tag{4.61} \end{align} $$

and we observe that this strategy requires a value of the temperature \( \theta_0 \) at a node 0 which is not included in Figure 55. However, this is not a problem since \( \theta_0 \) may be eliminated using the identity obtained from (4.61): $$ \begin{align} \theta_0 = \theta_2 \tag{4.62} \end{align} $$

Due to the zero gradient boundary condition, the first of the N equations, represented on generic form by (4.55), takes the particular form: $$ \begin{equation}< -(2+\beta^2h^2)\;\theta_1 + 2\theta_2 = 0 \tag{4.63} \end{equation} $$

This first equation (4.63) is the only equation which differs from the resulting equation system for prescribed boundary conditions in (4.58). All the coefficients \( a_i \), \( b_i \), and \( d_i \) are the same as in for the prescribed temperature version in (4.59), except for \( c_1 \): $$ \begin{align} c_1 &= 2, \qquad c_i = 1, \qquad i = 2,...N-1 \tag{4.64} \end{align} $$

Mixed boundary conditions. Version 2

An alternative version 2 for implementation of the zero-gradient boundary condition may be obtained by using a forward approximation for the gradient as given by (4.55) for a generic node i: $$ \begin{equation} \frac{d\theta}{dx}\bigg|_i \approx \frac{-3\theta_i + 4\theta_{i+1} - \theta_{i+2}}{2h} \tag{4.65} \end{equation} $$

which takes the following form for node 1 where is should evaluate to zero: $$ \begin{equation} \frac{d\theta}{dx}\bigg|_1 \approx \frac{-3\theta_1 + 4\theta_2 - \theta_3}{2h} = 0 \tag{4.66} \end{equation} $$

From equation (4.66) we see that \( \theta_3 \) may be eliminated by: $$ \begin{align} \theta_3 = 4\theta_2 - 3\theta_1 \tag{4.67} \end{align} $$

The first difference equation (4.55) in which \( \theta_3 \) occurs, is the one for node 2 $$ \begin{equation} \theta_1 - (2+\beta^2h^2)\;\theta_2 + \theta_3 = 0 \tag{4.68} \end{equation} $$

and we may eliminate \( \theta_3 \) from equation (4.68) by substitution of (4.67): $$ \begin{equation} 2\, \theta_1 - (2 + \beta^2h^2) \; \theta_2 = 0 \tag{4.69} \end{equation} $$

This is the first equation in the system of equations and the only one which differs from (4.58). Rather than (4.59), we get the following equations for the coefficients: $$ \begin{align} b_1 &= 2 & b_i &= (2+\beta^2h^2) & i &= 2,...,N \tag{4.70} &\\ c_1 &= -(2-\beta^2 h^2) & c_i &= 1 & i &= 2,...,N-1 & \tag{4.71} \end{align} $$

For convenience we summarize the resulting system of equations by: $$ \begin{align} & 2\, \theta_1 -(2-\beta^2h^2)\;\theta_2 = 0, &i & = 1 \nonumber \\ &\theta_{i-1}-(2+\beta^2h^2) \; \theta_i + \theta_{i+1} = 0, &i &= 2,3,...,N-1 \tag{4.72}\\ &\theta_{N-1} - (2+\beta^2h^2)\;\theta_N = -1, & i &= N \nonumber \end{align} $$

Note that we in this case, with a forward difference approximation of a gradient boundary condition, had to combine the approximation of the gradient with the difference equation to get a resulting tri-diagonal system. The the program ribbe2 we solve the mixed boundary conditions with both Version 1 and Version 2.

# src-ch3/section321/ribbe2.py

import numpy as np
import scipy as sc
import scipy.linalg
import scipy.sparse
import scipy.sparse.linalg
import time
from numpy import cosh

#import matplotlib.pyplot as plt
from matplotlib.pyplot import *
# Change some default values to make plots more readable on the screen
LNWDT=3; FNT=20
matplotlib.rcParams['lines.linewidth'] = LNWDT; matplotlib.rcParams['font.size'] = FNT

# Set simulation parameters
beta = 3.0
h = 0.001               # element size
L =1.0               # length of domain
n = int(round(L/h))  # # of unknowns, assuming known bndry values at outlet
x=np.arange(n+1)*h      # x includes min and max at boundaries were bc are imposed.


#Define useful functions

def tri_diag_setup(a, b, c, k1=-1, k2=0, k3=1):
    return np.diag(a, k1) + np.diag(b, k2) + np.diag(c, k3)

def theta_analytical(beta,x):
    return np.cosh(beta*x)/np.cosh(beta)

#Create matrix for linalg solver
a=np.ones(n-1)                  # sub-diagonal
b=-np.ones(n)*(2+(beta*h)**2)   # diagonal
c=np.ones(n-1)                  # sub-diagonal
#c=a.copy()                      # super-diagl, copy as elts are modified later
#c=a
# particular diagonal values due to derivative bc
version=1
if (version==1):
    c[0]=2.0
else:      
    b[0]=2.0
    c[0]=-(2-(beta*h)**2)
    print('version 2') 
                     
A=tri_diag_setup(a,b,c)

#Create matrix for sparse solver
diagonals=np.zeros((3,n))
diagonals[0,:]= 1.0               # all elts in first row is set to 1
diagonals[0,0]= 1.0               
diagonals[1,:]= -(2+(beta*h)**2)  
diagonals[2,:]=1.0

if (version==1):
    diagonals[2,1]= 2.0               # index 1 as the superdiagonal of spdiags is not used,  
else:
    diagonals[1,0]=2.0                # Sets the first element in the main diagonal  
    diagonals[2,1]= -(2-(beta*h)**2)  # index 1 as the superdiagonal of spdiags is not used,   super-diagonal

A_sparse = sc.sparse.spdiags(diagonals, [-1,0,1], n, n,format='csc') #sparse matrix instance

#Crete rhs array
d=np.zeros(n)
d[-1]=-1

#Solve linear problems
tic=time.time()
theta = sc.sparse.linalg.spsolve(A_sparse,d) #theta=sc.linalg.solve_triangular(A,d)
toc=time.time()
print('sparse solver time:',toc-tic)

tic=time.time()
theta2=sc.linalg.solve(A,d,)
toc=time.time()
print('linalg solver time:',toc-tic)

# Plot solutions
plot(x[0:-1],theta,x[0:-1],theta2,'-.',x,theta_analytical(beta,x),':')
xlabel('x')
ylabel(r'Dimensionless temperature $\mathregular{\theta}$')
legend(['sparse','linalg','analytical'])
show()
close()
print('done')

The relative error is computed from \( \varepsilon_{rel} = |(\theta_{num}-\theta_{analyt})/\theta_{analyt}| \), and the results of the computations are given in the table below:

We observe that the two versions for zero-gradient boundary condition yields approximately the same result except close to \( x=0 \), where the forward difference is somewhat better. In section (6.3.1 Neumann boundary conditions) we will take a closer look at the accuracy of gradient boundary conditions.

Mixed boundary conditions. Version 3

Rather than the enumeration in Figure (55), we may use the enumeration in Fig. (54) with \( x_i = i\;h \ , \ i=0,1,...,N+1 \) where \( h=\frac{1}{N+1} \) such that we must take \( i=1 \) in (4.55) to get: $$ \begin{equation*} \theta_0 - (2+ \beta^2h^2) \;\theta_1 + \theta_2 = 0 \end{equation*} $$

The boundary condition in (4.65) the becomes: $$ \begin{equation*} \frac{d\theta}{dx}\bigg|_0 = \frac{-3\theta_0 + 4 \theta_1 - \theta_2}{2h}=0 \end{equation*} $$

from which an elimination equation for \( \theta_0 \) may be obtained: $$ \begin{equation} \theta_0 = 4(\theta_1 - \theta_2)/3 \tag{4.73} \end{equation} $$ Consequently, equation (4.73) may be used to eliminate \( \theta_0 \) from the first difference equation: $$ \begin{equation} \tag{4.74} -(2 + 3\beta^2h^2) \;\theta_1 + 2\theta_2 = 0 \end{equation} $$

With approach \( \theta_0 \) is not solved with the equation system, but retrieved subsequently from equation (4.73).

• «
• 1
• ...
• 34
• 35
• 36
• 37
• 38
• 39
• 40
• 41
• 42
• 43
• 44
• 45
• 46
• 47
• 48
• 49
• 50
• 51
• ...
• 79
• »