instead oftemperature, on the right-hand side ofEq. (10) takes into account the stabilizing nature oflatent heat in cases with phase transformation. That is, energy ...
FSI-95-TN43
ENHANCEMENTS TO HEAT TRANSFER AND SOLIDIFICATION SHRINKAGE MODELS IN FLOW-3D®
M. R. Barkhudarov Flow Science, Inc. July 1995
I. INTRODUCTION This work outlines recent developments in heat transfer, conduction and solidification algorithms of the general purpose CFD code FLOW-3D, which resulted from efforts to improve its performance and accuracy. Section II describes a new, fully implicit treatment of heat fluxes in the fluid and obstacle energy equations. Section III concerns the Rapid Solidification Shrinkage (RSS) model, which in many cases offers a substantial increase in computational efficiency compared to the existing dynamic shrinkage model without loss of accuracy. Special sub-time-stepping during implicit calculations maintains accuracy even when large time step sizes are employed. Finally, Section IV outlines enhancements to the fluid/obstacle heat transfer algorithm that lead to improved accuracy. All these additions are illustrated by simple examples at the end of each section.
II. IMPLICIT TREATMENT OF HEAT TRANSFER AND CONDUCTION 1. Overview FLOW-3D solves the energy equation which includes terms describing conduction and fluid/obstacle heat transfer. These terms are numerically evaluated explicitly giving a simple and accurate solution procedure but also requiring the time step size to be kept below a certain
1
limiting value to prevent numerical instability. For example, in a uniform one-dimensional mesh with cell size ~, the time step size must satisfy the limiting condition
I1t ~ CONHT· min ( ~:'
2f)
(1)
where a = p~ is the diffusion coefficient and p = :c in which p, C, k and h are the density, specific heat, thermal conductivity and heat transfer coefficient, respectively. CONHT is a safety factor used in FLOW-3D to ensure stability of two- and three-dimensional transient solutions and usually equals to 0.45 [1]. Since the minimum on the right-hand side ofEq. (1) is estimated over the whole mesh, it is obvious that this time step size limitation could be too severe to allow an efficient numerical solution. This may be the case in highly·distorted meshes where the time step size is controlled by the stability limit in the smallest cell, even though the temperature may not be changing in this region. The existing optional locally implicit treatment of the heat fluxes in FLOW-3D, employed by setting IMPHTC=l [1], is based on an implicit representation of only the temperature of the cell for which the equation is currently being solved, while leaving all other temperatures used to compute the fluxes explicit. This algorithm yields a diagonal solution matrix and does not require iterations to obtain the solution. Despite its simplicity and unconditional stability, the locally implicit method is known to be inaccurate. One of its main drawbacks is that truncation errors in the numerical evaluation of the fluxes, which are first order with respect to ~t, accumulate with each tilne step possibly leading to significant errors in the final solution [2]. Therefore it can be used to advantage only in a limited number of situations. In this note we present a modification to FLOW-3D that eliminates the numerical stability condition expressed by Eq. (1) and provides an accurate solution. The modification consists of making the heat fluxes fully implicit functions of the new-time temperatures. Since implicit fluxes are not always needed, this addition has been added to the code as an option that can be selected by the user through input data (IMPHTC=2). A more detailed discussion of the differences between explicit and implicit approximations and iterative convergence of implicit solution methods can be found in [3] which describes an implicit algorithm for viscous stresses in the Navier-Stokes equations. Addition of implicit fluxes in the energy equation means that these terms, which previously were explicit and only needed to be computed once per time step, must now be made part of an iterative solution strategy. The complete solution procedure is outlined below, including some limitations with respect to other options in the FLOW-3D code.
2
2. Implicit Solution Outline Let us consider a one-dimensional linear heat conduction equation for temperature T with a constant coefficient of heat diffusion k (see [3])
(2) where E(T) is the fluid energy, which for an incompressible fluid is a linear function of temperature. Err) may also include heat of transformation and a second fluid component for two fluid problems. An unconditionally stable finite-difference approximation for this second order differential equation is
E(T~+I) == E(T:) + d (T~+1 - 2T~+1 1 1 1+1 I
I + T:+ ) I-I
(3)
where
(4) in which & and !1t are space and time increments [2]. A superscript n denotes the time t = nf1/, while a subscript i denotes the space location x = i&. Although this difference equation could be solved directly by a tri-diagonal solver, such direct solutions are not always possible in two- and three-dimensional cases. Solutions are therefore sought by iterative techniques. A simple way of setting an iteration scheme is given by
(5) Here we have used a superscript k to denote the iteration level. Level k=O corresponds to time level nf1/. As k increases, one hopes the iteration process will converge to a steady state in which the values of T are no longer changing. When this happens, the resulting values are the desired solution at time (n + 1)11/ ofEq. (3). Since Err) is a linear function, Eq. (5) offers a simple way of obtaining temperatures at iteration level k+1 using values at level k. In the simplest case of
E=pCT
(6) (7)
Eq. (7) constitutes a Jacobi iteration process. The new iteration value is entirely based on the previous iteration values. It can be shown that for finite values of p~ this iteration procedure converges to a unique solution [3,4]. 3
3. Convergence Criterion and Time Step Control A solution is considered converged when
max 1T k+1 - Tkl < EPSIHT
(8)
where
EPSIHT= EPSHTC· max(0.001 dTk+1, 0.0001 '1'+.1) mIn dTk+l
= _1 max IE k+1 -Enl pC
yk+.1 = min(Tk+1) mIn
(9) (10) (11)
EPSHTC is a user specified constant (defaults to 1.0). Min and max in Eqs. (8), (10) and (11) are evaluated over the entire mesh including obstacles and boundaries. The use of energy, instead of temperature, on the right-hand side ofEq. (10) takes into account the stabilizing nature of latent heat in cases with phase transformation. That is, energy may change significantly between times (' and ('+1 due to latent heat release during solidification, but at the same time temperature may vary little. For example, for pure material solidification the cell temperature stays constant while the solidification front passes through it. Equation (10) then ensures that fewer iterations will be needed in that case to reach temperature solution convergence. The presence of the Tm;n in the right-hand side of Eq. (9) ensures that the code does not iterate excessively to reach an accuracy too fine compared with the minimum temperature in the domain. Here the cutoff is suggested to be 0.0001· EPSHTC· T min. It is implicitly assumed that the value of T min is somehow related to the actual magnitude of temperature variation in the flow. If, for example, the minimum temperature is 1000 degrees, then a cell temperature variation of one degree will be considered insignificant. If this is not the case, then the user can adjust the cutoff parameter by adjusting the value EPSHTC. Since the value of EPSIHT varies from time step to time step, its value is printed out to the HD30UT.DAT file, together with other diagnostics 'data, at time intervals specified by SPRTDT [1 ]. As in all other iteration processes in FLOW-3D, such as pressure and viscous iterations, a parameter that specifies the maximum number of iterations allowed for the implicit heat transfer algorithm is provided in the input namelist LIMITS. That limit is ITHTMX which defaults to 200. In addition, the input parameter IDTHT (namelist LIMITS, defaults to 50) is used to control the time step size. If the number of temperature iterations in a time cycle, ITHTC, is greater than IDTHT, then the time step size for the next cycle will be decreased proportionally to the ratio of IDTHT and ITHTC. This ensures efficient and accurate solution. 4
4. Examples
A simple one-dimensional heat flow example is described in this section to demonstrate the implicit algorithm. The properties of the fluid were assumed to be constant during the simulation and close to those of liquid aluminum, Density: p=2,500 kg/m 3 Specific Heat: C=1000 J/kgK; Thermal Conductivity: k= 100 W/mK.
;
The problem involved a transient conductive heat flux from a boundary at xb=O.O at a constant temperature of Tb =1000oK into the fluid which was initially at a uniform temperature T;n=300oK. The fluid region extended from x}=O.O to x]=1.0 m. No phase transformation was included in this problem. The solution was sought for 0< t 0 and x >0 is [8] T(x,t)-T;n To-Tin
with X=
= 1-
;4/TV\ - [ erJ\..t'1.)
at) J[ 1 - t(X + -hj(ii) - ]
exp (hx -k + -hk 2 2
er
k
(28)
_x_
2j(ii .
For the numerical simulation the length of the sand domain was set at 200 mm, which was a good approximation of semi-infinity during the first hour of the simulation for Density: Specific Heat: Thermal Conductivity: Heat Transfer Coefficient:
p=1,500 kg/m 3 ; C=1675 J/kgK; k=O.75 W/mK; h=l,OOO W/m 2K. 15
--~
--------------------------------------
A uniform mesh with a 20 mm spacing was employed to solve the energy equation in the sand. The coarseness of the mesh was used to magnify the differences between the results of the old and new methods. Figure lOa shows temperature histories computed at the center of the interfacial cell (at x=IO mm from the interface) compared with the analytical solution. The result of the old heat transfer method is obviously less accurate. More over, it does not even converge to the analytical solution at the end of the simulation at (=3500 s. In contrast, the temperature predicted using the modified algorithm is much closer to the analytical solution and practically coincides with it at the end of the simulation. Figure lOb shows the ratio of predicted to analytical heat fluxes (i.e. normalized heta flow) at the interface. It can be seen that the old method largely overpredicts the heat flux at the initial stage of the simulation, « 150 S, while the new one underestimates it in the same period of time. Finally, predicted, normalized total energies fluxed through the interface are plotted in Fig. IOc. The result of the new method converges to the analytical solution much faster as the simulation progresses than the one computed using the old method. At the end of the simulation, (=3500 S, the energy predicted by the old method is still 30% above the analytical solution, while the new method is practically identical to the theoretical result. The results in this section are computed for a case when the fluid/wall interface lies exactly on a mesh face. What happens if we place the interface inside a mesh cell and keep the same resolution? Figures Iia and lIb show the interfacial heat flux and the total energy fluxed through the interface at (=130 S from the beginning of the simulation, as predicted by the two methods and compared with the analytical solution, as a function of the fractional open volume in the cell, VF• It can be seen that the results of the old method for the heat flux vary by than 50%, compared to the analytical solution, as VF changes between 0 and 1. At the same time the new method predictions vary by about 15%, a significant impr9vement in accuracy. The error is still substantial with the improved method because of the coarseness of the mesh used for this test, tu = 20 mm. These errors can be significantly reduced by using a finer mesh. A more typical grid spacing employed when computing heat fluxes in sand are in the range of 2 to 5 mm. According to Eq. (25), truncation errors in evaluating the interfacial heat flux decrease as tu 2 , therefore, by changing tu from 20 to 5 mm the maximum error should be decreased from 15% to less than 1%. For tu = 2 mm the maximum error should be 0.15%. An actual simulation for the latter case gave an error of 0.45%. It is larger than the estimated 0.15% due to additional truncation errors introduced by the first-order discretization method employed in the code for the conductive terms in the energy equation.
16
a
15.0
0.0 0.0
20.0
40.0
60.0
80.0
x h
Figure 1. Uniform (a) and non-uniform (b) meshes used in one-dimensional heat flux calculations in Section 2.
850
..
..
740 t
e m p e
r a t u r e
:(~.;'"~~.~~.I' .. f
630 ........................ :
,f , ,
520
ppP pP
t : 1 : I
:
t .p
I.: t
/.
410
:/ o
·······1········ ~
··········1··
.
~
E
i ~
.. ......•... t..;..........•i............
K
300
p:
~ .
•
:
•
:
.
.
:
i ~
:
~
"\".
·················r
,
,
.
.
.
.
iii
I
300
400
100
200
500
time, s a
analytical solution
b
implicit solution, y=l
c
implicit solution, y=10
d
'locally' implicit solution, y=l, imphtc=l
Figure 2. Comparison of the 'locally' implicit solution (IMPHTC=1, d) with fully implicit solutions at x=80 mm. The implicit solution for y=10 (c) is more accurate than the locally implicit for y= 1 (d).
850
.
- :' - -
• I;
t e m
740
..
r
~:::.-·::::::::·:.. :.:·.·.·~.·.·.·r:
.
:
..
·..·1
r-
.
~
!;
;!.
:
P",
. · · · . ·. · · · · · · . ·1
:
:'"
:/~l
630
a t
.
~
:
./~).I: 1··
p
e
~
•••••••••••••• 1"
:
.
u r
............. ···..··..····· ··
e K
~
I
410
300
o
100
:
200
300
400
500
time, s
analytical solution implicit solution, implicit solution,
o-=10, a=10,
epshtc=1.0 epshtc=O.l
Figure 3. Influence of the value of the convergence criterion on the implicit solution accuracy at x=80 mm. The convergence criterion is ten times smaller for the more accurate solution (dashed line).
850
t e
740
.
e r
a t u r e
'-'c /
630
. . . . . . . . . . . . ,.4·~·y·
