- Email: [email protected]
An efficient numerical method for the solution of pure convective transport problems with split boundary conditions
Chemical Engineering Science, 1966, Vol. 21, pp. 337-342.
Pergamon Press Ltd., Oxford.
Printed in Great Britain
An efficient numerical method for the solution of pure convective transport problems with split boundary conditions E. H. HERRON* and D. U. VON ROSENBERG Tulane University, New Orleans, Louisianna (Received 29 April 1965; in revisedform 31 August 1965) Abstract-The equations describing the transient behavior of a counter-current heat exchanger have been solved numerically by using finite differences centered between the grid points in both space and time. The space index is shifted for one dependent variable so the equations take a form for solution by an available algorithm for a bi-tri diagonal system. With an increment ratio chosen to minimize truncation error, the solution agreed almost exactly with an analytic solution which applied until the input discontinuity leaves the exchanger, even for a finite discontinuity in temperature at the inlet. The method has also been used for three dependent variables; the resulting equations were solved by an algorithm for the tri-tri diagonal system developed by one of the authors. INTRODUCTION
of obtaining dynamic mathematical simulations of the systems encountered in chemical engineering is rapidly becoming appreciated. When considering such simulations we may profitably separate all systems into two broad classifications: conductive transport processes and convective transport processes. We classify as conductive transport those processes in which a significant diffusion, or dispersion term is present. This means mathematically that the process is described by parabolic partial differential equations. Processes without significant diffusion we classify as convective transport problems, whose mathematical models are generally hyperbolic partial differential equations. The techniques for solving parabolic equations by numerical methods are fairly well established. However, the solutions to hyperbolic equations, with their inherent first order discontinuities are considerably more difficult to obtain. This paper applies an implicit differencing scheme designed specifically to treat hyperbolic equations and presents a method of efficiently solving the resulting difference equations without iteration. THE VALUE
EXAMPLEOF A CONVECTIVE TRANSPORT PROCESS
The numerical method to be described here can perhaps be most easily understood by studying its * Present Address, Esso Production
application to a particular system. As an example of a pure convective transport process, let us consider the dynamic operation of a countercurrent, double-pipe heat exchanger. Neglecting all conduction, the heat capacitance of the tube walls, and all radial gradients, assuming that both flow rates and all physical properties are constant, we obtain the equations describing the exchanger’s operation, which are
au - = at
au
at-
-
au z
av y2(u-u)+ IQ&-
(lb)
The boundary conditions are u(x,O)=
u(x)
u(0, t) = u(t) v(x, 0) = v(x) u(L, t) = v(t) where t = time
x = u= v= V, =
axial distance temperature of inner fluid temperature of outer fluid linear velocity
Research Company, Houston, Texas. 337
Y,(v - u) - VI
w
E. H. HERRON and D. U. VON ROSENBERG ,%ib.
Yl
=
i+M
UP/(pAC)i
U = heat transfer coefficient P = perimeter of the tube separating the two fluids p = density A = cross-sectional area to flow c = heat capacity L = exchanger length 1, 2 as subscripts refer to the inner and outer fluids respectively. To provide a check on the numerical solutions to follow, Eqs. (1) were solved analytically [l] for the inner fluid temperature above steady state temperature after a step change of u, at the exchanger inlet (x = 0). The solution is
se-prJn’?;px)2] t-7.X
u = uoemax
x
BX
- @x)‘]} x 4C~JC~”
n-l
j-l
FIG. 1.
Physical
j
representation encing.
of centered
differ-
and Thomee have called this analog the centered difference method. In this method difference approximations are made assuming that all variables are centered in both time and distance, as illustrated in Fig. 1. Thus,
dz + ewppX
for
j+r/z,n+1/2 x
andu=OforO
1
a
(2)
au 0ax
where c1= (V,Y, - V,Y,)/2T/,V,
lj2
uj+t n+1
- Uj+l.n
* At
Ujn+l
+
’
-
uj,.
At
(3)
j+l/Z,n+l/Z
B = (Vi + Frz)/2V1’,2 Y = (6 - V,)/2V,v, P = (v, y, - Vl y2)/(Vl + v,)
a 1/4(uj+l,n+i
CT2= 4F/,V,Y,Y,/(‘v, + V2)2
+ Uj+l,n +
uj,,+l
+ uj,.)
t5)
Let us number the mesh points in Fig. 1 from 0 to J (J = L/Ax) when referring to the inner fluid of first kind of order 1 and from 1 to J + 1 when referring to the outer It is apparent from (2) that there is a first order fluid. Since the boundary conditions give uO,”and discontinuity with t at every x. Conversely, there v,+i ,“, the points at which the temperatures are is a first order discontinuity in the “temperature to be calculated are numbered from 1 through J in profile”, the temperature distribution with length, both grids. With this numbering uj_112,n+l12 and at every t < L/V,. ~)~+t,~,~+~,~are at the same physical point, and Eqs. (1) become the following implicit difference equations when (3), (4) and (5) are applied to u, v, NUMERICAL SOLUTIONBYTHECENTERED DIFFERENCE and their partial derivatives : METHOD II = modified Bessel function
Equations (1) may be efficiently solved numerically using finite difference equations based on a differencing analog first proposed by WENDROFF [2] and later extended by THOMEE[3, 41. Wendroff 338
An efficient numerical method for the solution of pure convective transport
-
problems
with split boundary
conditions
r,
-“j+l,n+l
2
= & C”j,n+ Uj-l,n) + + (“j+l,n + uj.n - Uj,n-
ui- 1,“)-~(ui,i
-
(64
uj-l,n)
i _
O-
j=l,2
r,
-20 0
-- Y2uj.+1+
--Uj-l,n+l
2
2
-I-‘+ At
( (
, ***,J
y2 ---
- uj-
1.3+
Ax)
$
vj+1
1 +
CPU. J
+
/j!3)u. J J
,+I
+ c!~‘u.,+I I
I
I
It
6
7
8
9
lo
ft
(6b)
"j.0)
by’vj +
+ay'oj_l+
+ c% J a!3)uj_1 J
I 5
Table 1 shows that as few as five distance increments can be used to effect even faster simulations with only a modest sacrifice in the accuracy of the results. The results in Table 1 are from the problem illustrated in Fig. 2.
n+1 '
Cvj+ 1,n -
b!')u. + 1 J
I 4
FIG. 2. Temperature profile of inner fluid for 100°F step change at exchanger inlet after 8 sec. Implicit numerical solution vs. analytic solution.
with uo,“+r and ~.,+.r,~+rknown. Thus there are 2 x J linear equations in 2 x J unknowns at every time level at $ 1 of the general form a!‘)uj_ .I
I
3
(
v2
2
I
Length
"j,n+ 1 +
1 + t+
I
I2
j+l-
+ajP)~j_,
+ c% , j+1=
- d(J) J +
Distance
5 Increments
0.0 f:Z
10000 67.04 81.87
6.0 8.0 10.0
53.93 55.80 -11.39
10 Increments
20 Increments
Analytic
lOO+IO 67.06 81.91
lOOGO 67.08 81.91
lOOGO 67.08 81.92
54.95 56.97 1.34
54.91 56.97 0.01
54.91 44.93 0.00
(7a)
by’vj+ di2’
Table 1. Numerical results using different numbers of space increments
(7b)
where in (6) cy) = oy) = cj3) = ay) = 0, The bi-&i-diagonal system (6) can be solved simultaneously using the algorithm of DOUGLAS PEACEMAN and RACHFORD [5] to provide a temperature profile at the new time level. Except for some minor distortion at the discontinuity, numerical solutions obtained in this manner are in excellent agreement with the analytic results, as Fig. 2 indicates. For t > L/VI the solution is continuous and no distortion at all is observed. Furthermore, since no iteration is required, the implicit calculations proceed extremely rapidly. 339
In the presence of discontinuous boundary conditions, such as the step changes just considered, the difference Eqs. (6) converge most efficiently to the differential Eqs. (1) when At = AX/V, where V is the velocity of the fluid with the discontinuous boundary condition. In the absence of a discontinuity this requirement is unnecessary. Why this is so can be seen by considering (la) with Yr = 0; i.e. du
at’
-V$
E. H. HERRONand D. U. VONROSENBERG
We can show that centered differencing applied to (8) results in the expression ;
+
where w = Tube wall temperature y1.2 = (WP441.2
jl ~o($)‘*“-2”(~)” x (9)
&Zn+ ljU
1
y,
= (Wll(P&,
y,
=
W2lb44,
=
heat transfer film coefficient
h
’ (2r + 1)!(2n - 2r)! ( &x(*~-~~)&(*‘+~))
3 as a subscript refers to the tube wall
= -
v[$ +il ~o~)2r(~) (*n-2*)XEquation 1
$2n+lju
x (2r + 1)!(2n - 2r)! ( dx(2r+1)df(2”-2r)
)I
and further that
so that when At = Ax/V, (9) reduces to (8) exactly. Thus there is no truncation error when (8) is approximated by centered differencing. In more practical problems Yi is small compared to Vi but not zero, so some small truncation error exists. In a series approximation of a discontinuity the higher derivatives are significant, requiring that truncation error be minimized. Hence the convergence criterion that At = Ax/V for problems with discontinuous boundary conditions arises.
(1 Ic) contains no convective term, so we classify it as describing a conductive transport process. We will not therefore apply centered differencing to the MTvariable but instead will center w in time only. To accomplish this let us number u and u as before and number w from 3 to J + 3. Then Uj-l/Z,n, wj,,, and vj+l/2,n represent the same physical point on the x axis. With this numbering and with u and tr still centered in both time and distance, (1 lc) is approximated as
wj.n+1 -
wj
At
n ' =
-
y3C1/2(wj,n+l
-
1/4(uj,s + 1 +
-
y4C1/2(wj,n+l - wj,J -
-
1/4(0j+
Difference
l,n+ 1
+
Uj,n + Uj-l,n+l
+
vi,"+
equations
Wj,J-
+ Uj-l,n)l
-
(12)
*
1 + uj,n + “j+ 1,&l
similar
to
those
approxi-
mating (1) are written for (1 la) and (1 lb). three difference equations form a 3 x J
These
tri-triEXTENSIONSOF THE CENTERED DIFFERENCE METHOD diagonal system at each time step, which can be solved efficiently with an algorithm developed by Extension of the centered difference method to one of the present authors [I]. more complex systems is illustrated by considering Some results of a simulation obtained solving again the heat exchanger of Eqs. (1) without Eqs. (11) numerically are illustrated in Fig. 3, neglecting the heat capacitance of the tube wall where the outlet temperature of the inner fluid separating the two fluids. The mathematical model in response to a step change in its inlet temperature now becomes is plotted. Similar results from the model neglecting l3U au wall heat capacitance are also included. The model 5 = Yl(W - u) - v, ax (lla) including wall heat capacitance is seen to give a more realistic simulation, predicting an abrupt au - = Y*(w - u) -i V*$ jump in temperature when the hot fluid reaches the (lib) at outlet of the exchanger, followed by a gradual further increase in temperature as the exchanger aw -= - Y,(w- u) - Yd(W- u) (llc) wall heats up. at The centered difference equations can also be 4.6 0) = W(X) in addition to the initial condi- applied to hyperbolic equations with non-linear tions and boundary conditions of (1) coefficients, i.e. Yi and F’i are functions of the 340
An efficient numerical
method for the solution of pure convective transport
problems
with split boundary
conditions
strated, and it has been further shown that the centered difference method can be combined with differencing schemes which center variables in time only. Acknowledgements--The fellowshipsupportof the Humble Oil and Refining Company and the National Science
wall capacitance
Foundation and the enthusiastic contributions of Dr. R. E. C. WEAVERare acknowledged with appreciation.
NOTATION
6o0
I I
I 2
I 3
I
I
I
I
I
I
I
4
5
6
7
6
9
10
Time,
set
FIG. 3. Comparison of the mathematical models including and neglecting wall capacitance. Outlet temperature of the inner fluid in response to a step change in its temperature at its inlet.
dependent variables. Two methods similar to those proposed by DOUGLASand JONES [6] for parabolic systems have been developed for hyperbolic systems. These methods require a two-step iteration at each time step. Details of these methods will be presented in a future publication. In summary, this paper has shown the effectiveness of centered differencing in numerically solving convective transport problems when the difference equations generated are solved using &i-diagonal matrix inversion algorithms. Simulation of systems with up to three unknowns has been demon-
A = c = h= Zr = L = P = t = u= U= v= Yi,a = w= x = Yl,a = Y3 = Y4 = c4= /3 =
Cross-sectional area to flow Heat capacity Heat transfer coefficient Modified Bessel function of first kind of order 1 Exchanger length Perimeter of tube separating fluids Time Temperature of inner fluid Overall heat transfer coefficient Temperature of outer fluid Velocity of inner, outer fluid respectively Tube wall temperature Axial distance UPl(pA4l.a or (hZ7pAc)l.s W%/(pAc)s (hP)zl@Ac)a (KY1 - vlYz~/2vlva
(Vl f y = (Va -
VZ)/2VlV2 Vl)/2VlVZ
p = (VZYI -
VI Ya)/(Vl + VZ) in Eq.
p = Density in Yc cr2 = 4 Vl
va Yl Yz/( Vl +
(2)
vzy
Subscripts r = Inner fluid a = Outer fluid a = Tube wall
REFERENCES HERRON E. H. Distributed Parameter System Dynamics Of The Double-Pipe Heat Exchanger, Ph.D. Dissertation Chemical Engineering, Tulane University, 1964. PI WENDROFFB. J. Sot. ind. appl. Math. 1960 8 549. TIiOhfEEV. J. Sot. ind. appl. Math. 1962 10 229. :.i; THOMEEV. J. Sot. ind. appl. Math. 1963 11 964. [51 DOUGLASJ., PEACEMAND. W. and RACHFORDH. H. Petrol, Tram., 1959 216 297-308. WI DOUGLASJ. Jr., and JONESB. F. Jr. J. Sot. ind. appl. Math. 1963 11 195.
Ul
Rt%u&--Les auteurs r&solvent mn&riquement
les equations definissant le r&rime transitoire d’un echangeur a contre courant en utilisant des intervalles finis centres entre les points de la grille simultanement dans le temps et l’espace. La variable espace est remplacee par une variable dependante afin que les equations prennent une forme pouvant etre resolue par un algorithme applicable dans le systeme bi-tridiagonal. Avec un rapport d’increment choisi pour minimiser les erreurs graphiques, la solution est presque exactement en accord avec une resolution analytique appliquee jusqu’a ce que
341
in
E. H. HERRON and D. U. VON ROSENBERG la discontinuite introduite sorte du domaine de Nchangeur, meme pour une discontinuite finie de la temperature a l’entree. La methode a aussi et& utiliste pcur trois variables dependantes; les resultats des equations sont obtenues par un algorithme sur le systeme bi-tridiagonal, developpe par l’un des auteurs. Zusammenfassung-Die Gleichungen zur Beschreibung des Ubergangsverhaltens eines GegenstromWarmeaustauschers werden numerisch gel&t; dazu werden endliche ditferentielle Rechenschritte im Raum- und Zeitnetz durchgefiihrt. Der Raum-Index wird urn eine abhlngige Variable verschoben, so da13 die Gleichungsform mit Hilfe eines zuganglichen Algorithmus ftir das bi-tri-diagonale System l&bar wird. Wenn man das Inkrementen-Verhliltnis so wlhlt, dag der Abrundungsfehler minimiert wird, stimmt die so vorgenommene LGsung exakt tiberein mit einer analytischen Losung, die bis zum Austritt der Eingangs-Storung gilt. Diese Methode wurde such bei drei abhangigen Variablen angewandt; die erhaltenen Gleichungen wurden (von einem der Autoren) mit Hilfe eines Algorithmus fur das tri-tri-diagonale System gelost.
342
Pergamon Press Ltd., Oxford.
Printed in Great Britain
An efficient numerical method for the solution of pure convective transport problems with split boundary conditions E. H. HERRON* and D. U. VON ROSENBERG Tulane University, New Orleans, Louisianna (Received 29 April 1965; in revisedform 31 August 1965) Abstract-The equations describing the transient behavior of a counter-current heat exchanger have been solved numerically by using finite differences centered between the grid points in both space and time. The space index is shifted for one dependent variable so the equations take a form for solution by an available algorithm for a bi-tri diagonal system. With an increment ratio chosen to minimize truncation error, the solution agreed almost exactly with an analytic solution which applied until the input discontinuity leaves the exchanger, even for a finite discontinuity in temperature at the inlet. The method has also been used for three dependent variables; the resulting equations were solved by an algorithm for the tri-tri diagonal system developed by one of the authors. INTRODUCTION
of obtaining dynamic mathematical simulations of the systems encountered in chemical engineering is rapidly becoming appreciated. When considering such simulations we may profitably separate all systems into two broad classifications: conductive transport processes and convective transport processes. We classify as conductive transport those processes in which a significant diffusion, or dispersion term is present. This means mathematically that the process is described by parabolic partial differential equations. Processes without significant diffusion we classify as convective transport problems, whose mathematical models are generally hyperbolic partial differential equations. The techniques for solving parabolic equations by numerical methods are fairly well established. However, the solutions to hyperbolic equations, with their inherent first order discontinuities are considerably more difficult to obtain. This paper applies an implicit differencing scheme designed specifically to treat hyperbolic equations and presents a method of efficiently solving the resulting difference equations without iteration. THE VALUE
EXAMPLEOF A CONVECTIVE TRANSPORT PROCESS
The numerical method to be described here can perhaps be most easily understood by studying its * Present Address, Esso Production
application to a particular system. As an example of a pure convective transport process, let us consider the dynamic operation of a countercurrent, double-pipe heat exchanger. Neglecting all conduction, the heat capacitance of the tube walls, and all radial gradients, assuming that both flow rates and all physical properties are constant, we obtain the equations describing the exchanger’s operation, which are
au - = at
au
at-
-
au z
av y2(u-u)+ IQ&-
(lb)
The boundary conditions are u(x,O)=
u(x)
u(0, t) = u(t) v(x, 0) = v(x) u(L, t) = v(t) where t = time
x = u= v= V, =
axial distance temperature of inner fluid temperature of outer fluid linear velocity
Research Company, Houston, Texas. 337
Y,(v - u) - VI
w
E. H. HERRON and D. U. VON ROSENBERG ,%ib.
Yl
=
i+M
UP/(pAC)i
U = heat transfer coefficient P = perimeter of the tube separating the two fluids p = density A = cross-sectional area to flow c = heat capacity L = exchanger length 1, 2 as subscripts refer to the inner and outer fluids respectively. To provide a check on the numerical solutions to follow, Eqs. (1) were solved analytically [l] for the inner fluid temperature above steady state temperature after a step change of u, at the exchanger inlet (x = 0). The solution is
se-prJn’?;px)2] t-7.X
u = uoemax
x
BX
- @x)‘]} x 4C~JC~”
n-l
j-l
FIG. 1.
Physical
j
representation encing.
of centered
differ-
and Thomee have called this analog the centered difference method. In this method difference approximations are made assuming that all variables are centered in both time and distance, as illustrated in Fig. 1. Thus,
dz + ewppX
for
j+r/z,n+1/2 x
andu=OforO
1
a
(2)
au 0ax
where c1= (V,Y, - V,Y,)/2T/,V,
lj2
uj+t n+1
- Uj+l.n
* At
Ujn+l
+
’
-
uj,.
At
(3)
j+l/Z,n+l/Z
B = (Vi + Frz)/2V1’,2 Y = (6 - V,)/2V,v, P = (v, y, - Vl y2)/(Vl + v,)
a 1/4(uj+l,n+i
CT2= 4F/,V,Y,Y,/(‘v, + V2)2
+ Uj+l,n +
uj,,+l
+ uj,.)
t5)
Let us number the mesh points in Fig. 1 from 0 to J (J = L/Ax) when referring to the inner fluid of first kind of order 1 and from 1 to J + 1 when referring to the outer It is apparent from (2) that there is a first order fluid. Since the boundary conditions give uO,”and discontinuity with t at every x. Conversely, there v,+i ,“, the points at which the temperatures are is a first order discontinuity in the “temperature to be calculated are numbered from 1 through J in profile”, the temperature distribution with length, both grids. With this numbering uj_112,n+l12 and at every t < L/V,. ~)~+t,~,~+~,~are at the same physical point, and Eqs. (1) become the following implicit difference equations when (3), (4) and (5) are applied to u, v, NUMERICAL SOLUTIONBYTHECENTERED DIFFERENCE and their partial derivatives : METHOD II = modified Bessel function
Equations (1) may be efficiently solved numerically using finite difference equations based on a differencing analog first proposed by WENDROFF [2] and later extended by THOMEE[3, 41. Wendroff 338
An efficient numerical method for the solution of pure convective transport
-
problems
with split boundary
conditions
r,
-“j+l,n+l
2
= & C”j,n+ Uj-l,n) + + (“j+l,n + uj.n - Uj,n-
ui- 1,“)-~(ui,i
-
(64
uj-l,n)
i _
O-
j=l,2
r,
-20 0
-- Y2uj.+1+
--Uj-l,n+l
2
2
-I-‘+ At
( (
, ***,J
y2 ---
- uj-
1.3+
Ax)
$
vj+1
1 +
CPU. J
+
/j!3)u. J J
,+I
+ c!~‘u.,+I I
I
I
It
6
7
8
9
lo
ft
(6b)
"j.0)
by’vj +
+ay'oj_l+
+ c% J a!3)uj_1 J
I 5
Table 1 shows that as few as five distance increments can be used to effect even faster simulations with only a modest sacrifice in the accuracy of the results. The results in Table 1 are from the problem illustrated in Fig. 2.
n+1 '
Cvj+ 1,n -
b!')u. + 1 J
I 4
FIG. 2. Temperature profile of inner fluid for 100°F step change at exchanger inlet after 8 sec. Implicit numerical solution vs. analytic solution.
with uo,“+r and ~.,+.r,~+rknown. Thus there are 2 x J linear equations in 2 x J unknowns at every time level at $ 1 of the general form a!‘)uj_ .I
I
3
(
v2
2
I
Length
"j,n+ 1 +
1 + t+
I
I2
j+l-
+ajP)~j_,
+ c% , j+1=
- d(J) J +
Distance
5 Increments
0.0 f:Z
10000 67.04 81.87
6.0 8.0 10.0
53.93 55.80 -11.39
10 Increments
20 Increments
Analytic
lOO+IO 67.06 81.91
lOOGO 67.08 81.91
lOOGO 67.08 81.92
54.95 56.97 1.34
54.91 56.97 0.01
54.91 44.93 0.00
(7a)
by’vj+ di2’
Table 1. Numerical results using different numbers of space increments
(7b)
where in (6) cy) = oy) = cj3) = ay) = 0, The bi-&i-diagonal system (6) can be solved simultaneously using the algorithm of DOUGLAS PEACEMAN and RACHFORD [5] to provide a temperature profile at the new time level. Except for some minor distortion at the discontinuity, numerical solutions obtained in this manner are in excellent agreement with the analytic results, as Fig. 2 indicates. For t > L/VI the solution is continuous and no distortion at all is observed. Furthermore, since no iteration is required, the implicit calculations proceed extremely rapidly. 339
In the presence of discontinuous boundary conditions, such as the step changes just considered, the difference Eqs. (6) converge most efficiently to the differential Eqs. (1) when At = AX/V, where V is the velocity of the fluid with the discontinuous boundary condition. In the absence of a discontinuity this requirement is unnecessary. Why this is so can be seen by considering (la) with Yr = 0; i.e. du
at’
-V$
E. H. HERRONand D. U. VONROSENBERG
We can show that centered differencing applied to (8) results in the expression ;
+
where w = Tube wall temperature y1.2 = (WP441.2
jl ~o($)‘*“-2”(~)” x (9)
&Zn+ ljU
1
y,
= (Wll(P&,
y,
=
W2lb44,
=
heat transfer film coefficient
h
’ (2r + 1)!(2n - 2r)! ( &x(*~-~~)&(*‘+~))
3 as a subscript refers to the tube wall
= -
v[$ +il ~o~)2r(~) (*n-2*)XEquation 1
$2n+lju
x (2r + 1)!(2n - 2r)! ( dx(2r+1)df(2”-2r)
)I
and further that
so that when At = Ax/V, (9) reduces to (8) exactly. Thus there is no truncation error when (8) is approximated by centered differencing. In more practical problems Yi is small compared to Vi but not zero, so some small truncation error exists. In a series approximation of a discontinuity the higher derivatives are significant, requiring that truncation error be minimized. Hence the convergence criterion that At = Ax/V for problems with discontinuous boundary conditions arises.
(1 Ic) contains no convective term, so we classify it as describing a conductive transport process. We will not therefore apply centered differencing to the MTvariable but instead will center w in time only. To accomplish this let us number u and u as before and number w from 3 to J + 3. Then Uj-l/Z,n, wj,,, and vj+l/2,n represent the same physical point on the x axis. With this numbering and with u and tr still centered in both time and distance, (1 lc) is approximated as
wj.n+1 -
wj
At
n ' =
-
y3C1/2(wj,n+l
-
1/4(uj,s + 1 +
-
y4C1/2(wj,n+l - wj,J -
-
1/4(0j+
Difference
l,n+ 1
+
Uj,n + Uj-l,n+l
+
vi,"+
equations
Wj,J-
+ Uj-l,n)l
-
(12)
*
1 + uj,n + “j+ 1,&l
similar
to
those
approxi-
mating (1) are written for (1 la) and (1 lb). three difference equations form a 3 x J
These
tri-triEXTENSIONSOF THE CENTERED DIFFERENCE METHOD diagonal system at each time step, which can be solved efficiently with an algorithm developed by Extension of the centered difference method to one of the present authors [I]. more complex systems is illustrated by considering Some results of a simulation obtained solving again the heat exchanger of Eqs. (1) without Eqs. (11) numerically are illustrated in Fig. 3, neglecting the heat capacitance of the tube wall where the outlet temperature of the inner fluid separating the two fluids. The mathematical model in response to a step change in its inlet temperature now becomes is plotted. Similar results from the model neglecting l3U au wall heat capacitance are also included. The model 5 = Yl(W - u) - v, ax (lla) including wall heat capacitance is seen to give a more realistic simulation, predicting an abrupt au - = Y*(w - u) -i V*$ jump in temperature when the hot fluid reaches the (lib) at outlet of the exchanger, followed by a gradual further increase in temperature as the exchanger aw -= - Y,(w- u) - Yd(W- u) (llc) wall heats up. at The centered difference equations can also be 4.6 0) = W(X) in addition to the initial condi- applied to hyperbolic equations with non-linear tions and boundary conditions of (1) coefficients, i.e. Yi and F’i are functions of the 340
An efficient numerical
method for the solution of pure convective transport
problems
with split boundary
conditions
strated, and it has been further shown that the centered difference method can be combined with differencing schemes which center variables in time only. Acknowledgements--The fellowshipsupportof the Humble Oil and Refining Company and the National Science
wall capacitance
Foundation and the enthusiastic contributions of Dr. R. E. C. WEAVERare acknowledged with appreciation.
NOTATION
6o0
I I
I 2
I 3
I
I
I
I
I
I
I
4
5
6
7
6
9
10
Time,
set
FIG. 3. Comparison of the mathematical models including and neglecting wall capacitance. Outlet temperature of the inner fluid in response to a step change in its temperature at its inlet.
dependent variables. Two methods similar to those proposed by DOUGLASand JONES [6] for parabolic systems have been developed for hyperbolic systems. These methods require a two-step iteration at each time step. Details of these methods will be presented in a future publication. In summary, this paper has shown the effectiveness of centered differencing in numerically solving convective transport problems when the difference equations generated are solved using &i-diagonal matrix inversion algorithms. Simulation of systems with up to three unknowns has been demon-
A = c = h= Zr = L = P = t = u= U= v= Yi,a = w= x = Yl,a = Y3 = Y4 = c4= /3 =
Cross-sectional area to flow Heat capacity Heat transfer coefficient Modified Bessel function of first kind of order 1 Exchanger length Perimeter of tube separating fluids Time Temperature of inner fluid Overall heat transfer coefficient Temperature of outer fluid Velocity of inner, outer fluid respectively Tube wall temperature Axial distance UPl(pA4l.a or (hZ7pAc)l.s W%/(pAc)s (hP)zl@Ac)a (KY1 - vlYz~/2vlva
(Vl f y = (Va -
VZ)/2VlV2 Vl)/2VlVZ
p = (VZYI -
VI Ya)/(Vl + VZ) in Eq.
p = Density in Yc cr2 = 4 Vl
va Yl Yz/( Vl +
(2)
vzy
Subscripts r = Inner fluid a = Outer fluid a = Tube wall
REFERENCES HERRON E. H. Distributed Parameter System Dynamics Of The Double-Pipe Heat Exchanger, Ph.D. Dissertation Chemical Engineering, Tulane University, 1964. PI WENDROFFB. J. Sot. ind. appl. Math. 1960 8 549. TIiOhfEEV. J. Sot. ind. appl. Math. 1962 10 229. :.i; THOMEEV. J. Sot. ind. appl. Math. 1963 11 964. [51 DOUGLASJ., PEACEMAND. W. and RACHFORDH. H. Petrol, Tram., 1959 216 297-308. WI DOUGLASJ. Jr., and JONESB. F. Jr. J. Sot. ind. appl. Math. 1963 11 195.
Ul
Rt%u&--Les auteurs r&solvent mn&riquement
les equations definissant le r&rime transitoire d’un echangeur a contre courant en utilisant des intervalles finis centres entre les points de la grille simultanement dans le temps et l’espace. La variable espace est remplacee par une variable dependante afin que les equations prennent une forme pouvant etre resolue par un algorithme applicable dans le systeme bi-tridiagonal. Avec un rapport d’increment choisi pour minimiser les erreurs graphiques, la solution est presque exactement en accord avec une resolution analytique appliquee jusqu’a ce que
341
in
E. H. HERRON and D. U. VON ROSENBERG la discontinuite introduite sorte du domaine de Nchangeur, meme pour une discontinuite finie de la temperature a l’entree. La methode a aussi et& utiliste pcur trois variables dependantes; les resultats des equations sont obtenues par un algorithme sur le systeme bi-tridiagonal, developpe par l’un des auteurs. Zusammenfassung-Die Gleichungen zur Beschreibung des Ubergangsverhaltens eines GegenstromWarmeaustauschers werden numerisch gel&t; dazu werden endliche ditferentielle Rechenschritte im Raum- und Zeitnetz durchgefiihrt. Der Raum-Index wird urn eine abhlngige Variable verschoben, so da13 die Gleichungsform mit Hilfe eines zuganglichen Algorithmus ftir das bi-tri-diagonale System l&bar wird. Wenn man das Inkrementen-Verhliltnis so wlhlt, dag der Abrundungsfehler minimiert wird, stimmt die so vorgenommene LGsung exakt tiberein mit einer analytischen Losung, die bis zum Austritt der Eingangs-Storung gilt. Diese Methode wurde such bei drei abhangigen Variablen angewandt; die erhaltenen Gleichungen wurden (von einem der Autoren) mit Hilfe eines Algorithmus fur das tri-tri-diagonale System gelost.
342