Front tracking thermomechanical model for hypoelastic-viscoplastic behavior in a solidifying body

COMPUTER METHODS NORTH-HOLLAND

IN APPLIED

MECHANICS

AND ENGINEERING

81 (1990) 333-364

FRONT TRACKING THERMOMECHANICAL MODEL FOR HYPOELASTIC-VISCOPLASTIC BEHAVIOR IN A SOLIDIFYING BODY Nicholas ZABARAS, Department of Mechanical Engineering,

Yimin RUAN

~n~ver~~t~of Minnesota, Minneapoi~,

MN 55455, U.S. A.

Owen RICHMOND Alcoa Laboratories,

Alcoa Center, PA 15069, U.S.A

Received 26 May 1989 Revised manuscript received 21 October

1989

In assessing the quality of castings, a major consideration is the formation of cracks due to the induced thermal stress field. A means for understanding the casting process is the development and numerical implementation of mathematical models which account for all the heat transfer and deformation phenomena occurring in a solidifying body. In this paper a front tracking finite element method is used for the calculation of temperature and stress field development in a solidifying pure metal. The solid/liquid interface position and its velocity are considered as primary variables of the heat transfer analysis. A rate form of the virtual work principle and a rate dependent viscoplastichypoelastic constitutive model are employed to solve the equilibrium equations and to account for the hydrostatic stress state on the freezing interface and the fact that the material is in a state of residual stress immediately after solidification. Examples of the applicabiIity of the technique are given with the analysis of the solidification of pure aluminium under realistic cooling rates and material representation. The effects of melt pressure, cooling conditions and geometry of a continuously cast metal strand on the residual stress pattern are examined and reported.

1. Iutr~u~tion

A general theme in much of the previous work in the thermomechanical analysis of castings of pure metals is the lack of direct coupling between the stress development and the thermal field. Also limited attention has been given to the proper simulation of the mechanical and thermal conditions at the solid/liquid interface. Furthermore, a need exists to account for more realistic constitutive representation and cooling conditions. A complete thermomechanical description of the solidification of a metal includes the description of the heat flow, the thermal stress field development in the solid shell and the description of the coupling mechanisms, if any. Heat transfer boundary value problems with phase changes have been analyzed extensively. An important element of such an analysis is the proper handling of the solid/liquid front. The central problem is that on this continuously moving interface, velocity and heat transfer conditions have to be met. This means that standard numerical techniques based on a fixed domain (fixed discretization in space) have to be modified. For example, in the so called enthalpy method, one solves the governing equations on a fixed mesh in which interface 0045-7825/90/$03.50

0

1990 - Elsevier Science Publishers

B.V. (North-Holland)

334

N. Zabaras et al., Front tracking thermomechanical

model

conditions are accounted for implicitly via a new form of the governing equations [ 1, 21. Front tracking methods have also been applied extensively in the solution of phase change problems. Moving and deforming finite element front tracking studies include those of Lynch and O’Neill [3], Lynch [4] and Zabaras and Ruan [5, 61. Also moving boundary elements applied to Stefan problems have been employed by O’Neill [7] and Zabaras and Mukherjee [8]. In these methods, the mesh is deformed dynamically so that a continuous track of the solid/liquid front location is achieved. Further references and information on the above techniques and on the so called front fixing methods can be found in the above references and in a recent book by Crank [2]. In contrast to the number of papers on the heat transfer analysis of solidification, there is rather limited literature on the deformation and stress modeling of solidifying bodies. The most critical stages in such an analysis are in the early stage of solidification and when solidification is almost complete. The importance of the early stage is due to large stresses in the solid shell close to the metal/mold interface induced by large temperature gradients. These stresses are later reduced due to inelastic relaxation. For examples, see Richmond and Tien [9] and Tien and Richmond [lo]. The final stage of solidification is also important since poor design can result in partial solidification and cavities in the bulk of the shell. Cavitation is generally due to solidification contraction and usually appears in regions of the melt enclosed by solid before complete solidification has occurred. Over the temperature range considered, it is well-known that there is a large variation of the mechanical and thermal properties of the material [ll]. Further high temperatures require the use of rate dependent viscoplastic models. These constitutive models show time independent elasto-plastic behavior at room temperature, viscoplastic rate dependent behavior at high temperatures and a kind of fluid behavior close to or beyond the melting temperature [12, 131. Thermal stress analysis of solidification problems can be found in several papers including the paper by Kristiansson [14] who used a power creep law, by Weiner and Boley [15] who presented an analytical stress calculation with an elastic/perfectly plastic model, by Inoue and Wang [16] where Perzyna’s model [17] was used, by Thomas et al. [18, 191 where Wray’s model ]20] was adopted and by Richmond and co-workers [9-l 1, 21, 221, where a hyperbolic type of constitutive model was used. Experimental measurements of surface displacements and study of the air-gap formation are also given by Nishida et al. [23]. In addition to correct material modeling one must account for the proper boundary and interface mechanical conditions. At first, a proper physical modeling of the solidification problem requires that the stress state immediately after solidification be the same as that of the liquid before solidification [lo, 111. The deformation of solidifying material is very different from that of standard fixed bodies. It is the usual assumption that in bodies without phase changes the material deforms from an initial state of zero residual stresses. Nevertheless, a solidifying body develops residual stresses immediately after solidification and is never in a state of zero stresses (stress free state). Other complications which will not be analyzed in this paper include convective effects in the melt and their effect on the residual stress pattern and also dynamic coupling mechanisms between the stress and thermal analysis. This coupling can be due to the formation of air-gaps between the solid shell and the mold, for example. Such air-gaps can drastically change the heat transfer mechanisms and the friction boundary conditions on the solid shell/mold surface. Even though coupling mechanisms are important

N. Zabaras et al., Front tracking thermomechanical

model

335

in the stress pattern development, they will not be discussed further, so that attention can be given to thermal and mechanical phenomena on the solid/ liquid freezing front. This paper presents a general moving and deforming FEM analysis which carefully accounts for the solid/liquid front phenomena. The motivation for using a front tracking analysis lies in the belief that it is easy to account for conditions on a moving boundary whose location and velocity are part of the solution of the problem. A one way coupled thermal/stress model (thermal field affecting the stress field but not vice versa) is considered. This model is based on a rate form of the virtual work principle and it employs a hypoelastic-vis~oplastic constitutive model. The plan of this paper is as follows. First, the model problem is defined. Then the heat transfer part of the problem is considered, and a deforming-moving FEM analysis is presented. The paper continues with the FEM analysis of the deformation part and details for the consideration of the interface phenomena are given. The capabilities of the present method are demonstrated by analyzing one-dimensional and two-dimensional square and circular sections of a continuously cast aluminium strand as it is gradually solidifying from the molten state. The effects of different cooling and geometric conditions on the development of residual stresses are examined. Finally, qualitative discussion of the results and further potential applications of the proposed method are presented.

2. The heat conduction problem in a solidifying body 2. I. Governing

equations

A region Q, with boundary ~30~is initially occupied by liquid metal of temperature T,(x, y) equal to or above the melting temperature T,. Solidification is assumed to start at time t > 0 when the boundary aR, is cooled down to a temperature lower than the melting temperature. The situation after some time has elapsed is depicted in Fig. 1. Let diI,(t) be the isothermal freezing interface at time t. In the absence of heat sources, the governing heat conduction

Fig. 1. The solidification

process.

N. Zabaras et al., Front tracking thermomechanical

336

equations

model

are given by Carslaw and Jaeger [l] as

(2) and O,(t) denote the solid and liquid regions, respectively (O,(t) U L!,(t) = ii!*), at the point (x, y) E n,(t) at time t, TL(x, y, t) the temperature at the point (x, y) E O,(t) at time t, and p, c and K are the density, specific heat and conductivity, respectively. Subscripts S and L denote the solid and liquid phases, respectively. The boundary and initial conditions are defined as

where a,(t)

Il”,(x,y, t) the temperature

T(x, Y, f) = ~,(x, y, 0,

(x,

y) f ai&,,

(3)

am, Y, 0 = 4dx, Y, t) , (x, Y) E af&, , ano m, Y, 4 = T, , (x, Y> E an,(t> ,

(5)

W=O)=O,

(6)

K

S

aL(t=o)=Q),

where TO(x, y, t) is the prescribed temperature history on boundary a&, and qO(x, y, t) is the prescribed normal heat flux on boundary aGO, with a6?@,U afl,, = iM20.Also aT(x, y, t) /an, denotes the temperature gradient on the boundary N&,, at (x, y) E aa,, in the direction of the unit normal n,, to aK$, as shown in Fig. 1. Finally, energy balance on the freezing front gives the so called Stefan condition, K

s

aWy an

Y,

4 _

awx,Y, 0

K

I..

an

=pLV-n,

(8)

where n is a unit normal vector on the interface boundary aOr at a point (x, y) e a&?,(t) pointing away from the solid region, V is the velocity vector at the same point on the interface and L denotes the latent heat of fusion. Equation (8) connects the normal fluxes of the two sides of the solid/liquid interface with the normal interface velocity. Due to the solid/liquid interface motion, the problem defined by (l)-(8) is nonlinear even with temperature independent material properties. In summary, the primary unknowns of the above defined problem are the evolution of the solid/liquid interface velocity and position as well as the development of the temperature field in the solid and liquid phases. 2.2. A deforming and moving finite element formulation The regions OS(t) and O,(t) are changing with time and a front tracking finite element technique will require the use of finite elements which can move and deform as time passes.

N. Zabaras et al., Front tracking thermomechanical

337

model

This finite element motion can be introduced by considering the nodal coordinates to be functions of time. More specifically, let the region 0,(t) U O,(t) = f& be divided into E,(t) + EL(t) = E(t) finite elements with M,(t) + ML(f) - M,(t) = M(t) nodes. The E, elements of the solid phase share the same M, nodes on the solid/liquid interface boundary an,(t) with the E, elements of the liquid phase. In general, E,, E,, M, and M, are changing with time. Since the finite element nodes are permitted to move with time, the shape functions have implicit dependence on time through the nodal coordinates. For example, for each finite element e, one can write y, t) ,

2(x, y, t) = T;(t)@;@,

i = 1,2, . . . , M’ ,

(9

where summation on i is implied over the number of nodes in an element M’, T:(t) denote the nodal temperatures and @y(x, y , t) the element shape functions. Moving finite element formulations based on such an interpolation have been proposed earlier in [3-61. For completeness, a rather brief review of the involved ideas is given here. Let us assume the isoparametric transformation x = XF(WF(

5, ‘I) ,

Y = K(Wo>

4 >

(10)

where (Xe, Ye) are the nodal coordinates of element e and @T( 5, n) are the element shape functions on time independent master cordinates (5, r]) with - 1 d (5, n) d + 1. Using (9) and (lo), it has been shown by Zabaras and Ruan [6] that

where the second term on the right-hand side accounts for the fact that the point (x, y) at time is moving with velocity V(x, y, t). This velocity can be easily calculated by time differentiation of (10) assuming that the nodal velocities (XF, I’:) have been specified. Usually, these nodal velocities are calculated based on the velocity of the nodes at the freezing front and a rearrangement of the position of the internal nodes which guarantees a uniform mesh (Zabaras and Ruan [6]). Using (11) and employing a Galerkin type of weak formulation for (1) and (2), the following assembled system of equations is finally obtained as t

I, J = 1,2, ,~+(B,,+K,,)T,=F,,

C

where Z and J are the global nodal numbers (i, j = 1,2, . . . ) M’), respectively, with

. . . ,

M (sum on J) ,

corresponding

to element

(12) nodes

i

and j

(13)

B,, =

5 B; = - 2 jaPC@;(X,

e=l

e=l

c

y, t)V@;(x, y, t)- v(x, y) da

,

(14)

N. Zabaras et al., Front tracking thermomechanical

338

K,, =

5 K; = Ej-O

KV@f(x,

e=l

e=l

y, t) T@;(x,

model

y, t) dL! ,

e

where ani and afig, are the boundary segments of elements which belong to the solid/liquid interface afz, and the fixed boundary CM&,,,respectively. Note that the second term in (16) is due to the heat flux discontinuity on the interface aJIr and the matrix [B] is due to the motion of the finite element nodes. The calculations described in (13)-( 16) can be feasible only when the nodal velocities, front position and the domains L&(t) and L!,(t) are known at each time step. To perform the time integration in (12), an implicit scheme may be used. Denoting the nodal temperature vector at time t = t,_, + At with T”, n = 1,2, . . . , where At is a time step, a stable integration scheme similar to that proposed by Hughes [24] is

cn-I+, At

1

+ Yw,-1+, + Kn-1ty) T”

(17) L

J

where the subscript it - 1 + y indicates the reference time for the calculation of the temperatures and material properties as well as the reference time for the interface position and velocity required to calculate the integrals in (13)-( 16). For example, to determine Bn_r+,,, one calculates the material properties and the temperature field at time t and then performs the integration indicated in (14) based on the interface location and velocity at the same time t”,where F= (1 - y)t,_,

+ yt,

(18)

with 0.5 G y < 1. The location P” and velocity ? of a nodal point at the freezing front at time t^ are given as

p”=(l-y)P,_,+yP,,

(19)

-c;= (I- y)V,_, + yv, ,

(20)

PII = f At + P,_,

(21)

where

and the velocity V, is calculated iteratively. Suppose that fhe nodal temperature vector T”-’ as well as an estimate of the interface velocity at time t are known, one can then calculate T” from (17) by applying the boundary conditions on afiR, (eqs. (3) and (4)) and the temperature boundary conditions on da,(t) (eq.

N. Zabaras et al., Front tracking thermomechanical

339

model

(5)). This last boundary condition is enforced through a penalty method. ture field at time t, is calculated, one can return to (17) and calculate the right-hand side of (16) for each interface node. As it will become clear these terms can directly be used to update the interface velocity at time

Once the temperasecond term on the in the next section, t”.

2.3. A Galerkin weak formulation of the Stefan condition For a global conservation of energy at the freezing front, the following Galerkin type of weak formulation of the Stefan condition (eq. (8)) at time ?is introduced:

i= 1,2,.

. . , M;,

b = 1,2,.

. . ,

E, ,

(22)

where E, and M, are the number of boundary segments and boundary nodes on the solid / liquid interface, respectively, and ?Pq(x, y, t) and Mf are the shape functions and number of nodes for the interface segment b, respectively. The velocity of a point at the solid/liquid interface can be obtained as follows: V=Vflyf(x,y,t),

i=l,2,.

. . ,MF,

(23)

for each segment b, where Vf = (Xi, ji) are the interface nodal velocities on the segment b. After assembling, (22) becomes

i=l,2,...,

Mi , h = 1,2, . . . , Mf (sum on h) ,

(24)

where Nb is the normal to the boundary segment b on the solid/liquid interface and points away from the solid region (see Fig. 25 in Appendix A). Notice that the above equation defines the normal velocities to the solid/liquid interface. If the tangential nodal velocities are known, they can be used together with (24) to form a complete system of equations for the x and y components of the interface nodal velocities. In this work, the nodal tangential velocities are specified as V . ti = a, ,

i = 1,2, . . . , M, (no summation

on i) ,

(25)

where tj is the tangent vector (Fig. 25) to the solid/liquid interface da,(t) on the nodal point i, i=l,2,. . . , M, and ai can be selected a priori such that the geometry of the problem is properly represented. In forming (24) care has to be taken that undesirable spreading or bunching of nodes and tangling of the interface do not occur. After such possibilites have been taken care of, ai is chosen to be zero. Combining (24) and (25), a system of algebraic equations for the x and y components of the interface nodal velocities is obtained as Aijhkahk=Hjj,

i=l,2

,...,

M,,

h=l,2

,...,

M,,

j=l,2,

k=l,2,

(26)

340

N. Zabaras et al., Front tracking thermomechanical

model

where ahk is the kth component of the hth nodal velocity on the interface at time t”.The vector {H} and matrix [A] are given in Appendix A in a non-assembled form for one interface segment. Equation (26) can be used to iteratively calculate the nodal velocities on the solid/liquid interface front, once the evaluation of the right-hand side of (24) has been performed. Fortunately no extra calculations are involved, since these terms have been calculated from (17) as it is discussed at the end of Section 2.2. In other words, (26) provides an accurate representation of the physical phenomena on the freezing front as well as an easy methodology for the calculation of the freezing front motion. 2.4. Finite element nodal motion Continuously moving and deforming finite elements are used to analyze boundary value problems over regions which continuously change geometry. Assuming that the number of finite elements on a,, is fixed, it is quite clear that one should rearrange the positions of the nodes so that a non-distorted mesh is always present. Furthermore, to avoid interpolation of variables from an old mesh to a new mesh, one can assume that the remeshing is occurring because of a properly specified continuous nodal motion. In other words, each discrete variable is identified with a node rather than with a spatial location and of course, with time. A general remeshing technique based on a transfinite mapping with a bilinear projector was given by Zabaras and Ruan [6], and the associated problem of specifying nodal velocities was addressed. In this work, emphasis is given to problems with T, = T, for which a less ,expensive technique with time varying number of elements is used. To demonstrate the ideas, let us concentrate on the solid phase. The liquid is not examined, since T,(x, y, t) = T,,, for (x, y) E L!,(t). The region O,(t) is divided into two subregions. In one of the subregions, the nodes are always stationary and in the other one, the nodes lying on the solid/liquid interface are moving with respect to time (Fig. 2). The moving region includes only the elements of the solid phase next to the interface do,(t). These elements are allowed to change their sizes to adapt to the movement of the interface front. When the sizes of these elements in the deforming region become larger than a prescribed size, a new moving region is generated and the previous moving region becomes part of the fixed region at the same time. The new region created consists of elements of zero area. The new set of elements will not introduce singularities since during the next time step the solid/liquid interface will move and finite sized elements will be created. These elements will continue deforming over the coming time steps until a necessity exists for the introduction of new elements. Obviously, the matrix [B] in (12) is non-zero only on the moving region of O,(t). 2.5. Summary of the algorithm For each time step t, = t,_, + At, n = 1,2, . . . (to = 0) A. Determine if the sizes of the elements at t = t,_, in the deforming region are equal to or less than the prescribed value. 1. If yes, the interface nodes are considered to move with the velocity calculated at time t = t,_, while the non-interface nodes remain stationary. 2. If no, a new set of elements are introduced to the system based on the discussion in Section 2.4 (Fig. 2).

N. Zabaras et al., Front tracking therrnomechanical model

341

Fig. 2. The element configurations before and after new elements are introduced (the shaded area denotes a fixed region and the unshaded region a deforming region). (a) Element ~on~guration at time t,_ t when the size of the unshaded elements is maximum permissible.

Fig. 2. (b) Element configuration

after new elements of zero width are introduced

at time I,_ ,.

16

Fig. 2. (c) Element

configuration

at time t, = t,_, + bt. The width of newly introduced

elements

is now finite.

N. Zabaras et al., Front tracking thermomechanical

342

model

B. For each iteration i = 1,2, . . . , maximum number of iterations 1. Calculate all matrices in (17), then apply boundary conditions on 8Q, and solve for the temperature field at t,. 2. Using the calculated temperatures, return to (17) and calculate the second term on the right-hand side of (16), for e_ach finite element node on the freezing front. 3. Calculate a new interface velocity V using (26). 4. Update the location of the interface nodal points as (x, y)“+’ = ? At + (x, y)” .

(27)

5. Check the relative velocity error as

IKll - W-III IICLII ’

error =

(28)

where 1)- 11is the Euclidean norm of a vector and Q: denotes the calculated vector of the normal nodal velocities on da,(t) at t”at the ith iteration. If error < tolerance, stop iteration. If error > tolerance, continue iteration (return to B). End of iteration End of time stepping. If the temperature history at a specific point is required, an interpolation must be performed using the nodal temperature histories and (9). To start the algorithm, a small finite solid region is assumed at time zero. Usually the algorithm fails just before complete solidification has occurred due to the extremely dense nodal concentration.

3. Thermomechanical modeling of a solidifying body 3.1. Governing equations A quasistatic thermal stress theory as discussed by Boley and Wiener [25] is employed and attention is given only to the solid phase. Body forces and inertia forces are neglected, leading to the simple mechanical equilibrium equations aij,i(x7

_YYt> =

O

7

i, j = 1,2,3

(summation

on i) (n, y) E O,(t) ,

(29)

where aij are the Cartesian components of the Cauchy stress tensor at time t and the usual indicial notation is used. It is emphasized that the equilibrium equations (29) have to be satisfied at any time t in a region which continuously changes (grows) with respect to time. Assuming small strains and rotations, the total strain rate tensor iii is additively decomposed into an elastic, ii:, a thermal, .@: and an inelastic, 2: part, respectively. Therefore, iii =

$ + “; + i;

)

i,j=l,2,3,

N. Zabaras et al., Front tracking thermomechanical

where in terms of the rate of the displacement

model

field tii, Eii can be calculated

343

as

$(z.ij + l&) .

iii =

(31)

As already discussed in the introduction, the deformation of a solidifying body differs from that of a fixed body in that a solidifying body develops residual stresses immediately upon solidification and is never in a stress free state. More specifically, the stress state at the solid/liquid interface must always be purely dilatational, i.e., 0 Ull

=

0 uz2

=

0 uj3

=

-p(t)

)

CT;*

=

c&

=

c&

=

0

)

(32)

where p(t) is the fluid (melt) pressure at time t. In other words, the stress state immediately after solidification is the same as that in the liquid right before solidification. In this work, plane stress conditions are assumed and the constraint ai3 = -p(t) is not enforced explicitly. To account for such behavior, a hypoelastic constitutive equation is assumed to govern the deformation of a solidifying body. In this work, the following hypoelastic model is used: bij=Eijk,(T)iF,,

i, j,k,1=1,2,3,

where the superimposed dot denotes a time derivative and the temperature constants Eijk[(T) are given as

(33)

dependent

elastic

where A(T) and p(T) are the temperature dependent Lame parameters and aij is the Kronecker delta. The ‘initial’ residual stresses can now be treated as an integration constant arising when (33) is integrated in time. Dilatational thermal strains are assumed as (35)

where a(T) is the temperature dependent coefficient of thermal expansion and T, is a reference temperature at which the thermal strains are zero, i.e. T, = T,. Finally, a viscoplastic constitutive model is needed to prescribe the inelastic deformation. The following general form is assumed here ir=fij(oij,

q;, T),

i, j= 1,2,3,

(36)

where q: denote properly defined state variables, if any, for which evolution equations of the following form are given: cjt=gF(a,,q;,

T),

i, j=l,2,3.

(37)

Several viscoplastic constitutive models fall in the unified category of (36) and (37). Examples are the hyperbolic sine law with or without state variables [26], Anand’s model [12],

344

N. Zabaras et al., Front tracking thermomechanical

model

the model due to Hart [13] or that of Bodner and Partom [27]. Most of these models include important effects such as rate sensitivity, strain hardening and recovery, and they are valid for temperatures ranging from room temperature up to the melting point of the solidifying pure metal. Implementation of such models is relatively easy, but usually they are stiff in nature and require special attention in their numerical integration. The complete deformation problem is defined above by (29)-(37) together with proper traction and/or displacement boundary conditions on do,.,. This problem is coupled with the heat transfer problem as described earlier. Nevertheless, the heat transfer problem is considered as not being affected by the deformation problem. So, for example, the heat generated by inelastic deformation is neglected. Also, other coupling mechanisms such as the air gap between the mold and the casting, which is formed as the solidifying metal shell shrinks away from the mold, are neglected and the cooling boundary conditions are applied directly to the boundary afi, which is assumed stress free at all time. The effect of fluid flow in the melt, the induced shrinkage in the vicinity of the freezing front and non-isothermal solidification problems will not be examined here. 3.2. Finite element model The finite element analysis of the stress problem is performed with the same discretization of the solid region l&.(t) as that used for the heat conduction analysis and using the same shape functions @F. Let tii denote the nodal displacement rates in the solid phase, then the field of displacement rate ci can be approximated as ci= z e=l

@pe&f=Nri ,

(summation

on i) i = 1,2, . . . , M” ,

(38)

where ri is the column vector of nodal components tij and N is the matrix of shape functions. Spatial differentiation of the above equation defines the field of strain rate as d=Bti,

(39)

where B is the standard matrix defined by differentiation of the shape functions in N. one can write a weak statement of the Assuming that time stepping is performed, equilibrium equation (29) at time t, as

with a, the stress tensor (a vector form) at time t,, the body forces have been neglected and the boundary a& is assumed traction free. As discussed before, a hypoelastic constitutive law is used to account for a hydrostatic stress state of a material point at the time it solidifies. A rate form of the principle of virtual work such as that proposed by Morgan et al. [28] could be employed here. However, since numerical implementation of such a model leads to solutions which gradually shift away from equiIibrium, (40) will be used to derive a modified weak formulation in a rate form.

N. Zabaras et al., Front tracking thermomechanical

model

345

Let us concentrate on a material point where it is assumed that the stress state, a,,_, , at t,_, is known. Then the stress state, a,, at current time t, can be obtained using Euler’s backward integration scheme as a,, =

ci, At* +

(41)

on_, ,

where a, denotes the stress rate at time t, and At* is a time step defined as follows. If the material point at which the time integration is performed is solid at time t = t,_, , At* is equal to At(= t, - t,_l) and in the case that the material point is liquid at time t = t,_, , At* is defined as At* = t, - t,*_, , where t,*_l denotes the arrival time of the freezing front at the corresponding material point. Generally, in finite element calculations, the integrations in (40) are performed with Gaussian integration, therefore (41) will be evaluated at Gauss points. As discussed previously in the thermal model, the solid phase is divided into two sub-regions, a fixed and a moving region (Fig. 2). In the fixed region, the interface arrival time for all the Gauss points is earlier than time t, _ 1, therefore At* = At = t, - t, _ 1. In the moving region, if the arrival time at a Gauss point is earlier than t, _1, At* is equal to At(=t, - tn_l), and if the arrival time of the point is later than time t,_, , At* is equal to t, - t,*_l. Therefore, one can conclude that, if the front arrival time at a point is later than t,_, , the stress a,_, is equal to the melt hydrostatic stress state at arrival time and if the arrival time is earlier than t, _1, it is the stress calculated at t = t,_, . Combining (40) and (41) gives

I

fl,(L)

B’b,, At* d0 =

N’r,, dT -

I Jfl,(t,)

I fi,(t,)

B’q_,

dR

.

(42)

Finally, using (33), (39) and (42), the following rate form of the principle of virtual work is obtained:

I

%0,)

B’DB At* d&?ri =

I Jf-$O,,) +

I

fqr,)

N%,, dT -

I fgf”)

B’q_,

~$*+i~)Af*d& BtD(

da (43)

The solution scheme used to solve (43) for the displacement rate ti is similar to that introduced by Zienkiewicz and Cormeau [31]. After the nodal displacement rates are found, they are then used to calculate stress rates at element Gauss points. For the non-moving region, the stress rates are integrated at the Gauss points. For the moving region, the stress rates at the Gauss points are extrapolated to obtain nodal stress rates using the local least squares technique [6]. Finally, these node1 stress rates are integrated in a way similar to that indicated in (41). The difference between the time integration of stress and temperature rates should be emphasized. In the stress problem, the primary unknowns are nodal displacement rates and the integration of stress rates is performed at material points. But for the temperature problem where the primary unknowns are nodal temperatures, the integration of temperature rates is performed at nodal points rather than at material points.

346

N. Zabaras et al., Front tracking thermomechanical

model

3.3. Constitutive model In the simulations reported in this paper, a hyperbolic-sine constitutive law [lo, 111 is used to prescribe the inelastic deformation. This constitutive law has the following form: [sinh Be]” c ‘ij 7 i, j = 1,2,3 7 > T + 273 a

w

where A, B and C are material constants given in Table 2 for pure aluminium, 6 is the effective stress defined as 6 = d= and sii are deviatoric stresses defined as sij = aij f ukkSij. The temperature T is in degrees Celsius and C in degrees Kelvin. The above constitutive model was previously employed by Heinlein et al. [ll] for unidirectional solidification problems. 3.4. Summary of the algorithm for the thermal stress problem For each time step t, = t, _, + At, n = 1,2, . . . (t, = 0) A. Calculate the temperature field and front position as discussed in Section 2.4. B. For the moving elements only: 1. On the old mesh (at time t = tn_l), transfer the stress rate &,_i from the Gauss points to the corresponding element nodal points using a local least squares method and then integrate the nodal stress rates to obtain nodal stresses at time t = t,_, as discussed in Section 3.2. 2. Interpolate the nodal stresses found at step B.l to obtain the stresses at the integration points of the new mesh found in step A. C. Evaluate the thermal strain rates using (45), tir= 11

P-P1 At*

a(T)‘,

>

i, j = 1,2,3 ,

(45)

where ?’ and ?‘- ’ are the temperatures at a Gauss integration point obtained by interpolating nodal temperatures at time t = t, and t = t,*_ 1, respectively. D. To start iterations, if it is the first time step or if new elements were introduced, the thermal strain rates .6f is used as an initial guess of “$ + if, otherwise, previous estimates of z$: + 6: are used as an initial guess. E. Calculate the stiffness matrix and the load vector due to boundary conditions as B’DB At* dL! , F. For each iteration i = 1,2, . . . , maximum number of iterations, 1. Evaluate the load vector due to the thermal and inelastic strain rates, B’D(B= + dN) At* da {Fs) =In,,, 2. Solve (43) for ti and compute &,, = DBzi, - D(8= + kN).

.

the stress rate at Gauss points as (46)

N. Zabaras et al., Front tracking thermomechanical

model

347

3. Integrate the stress rate using (41). 4. Calculate the inelastic strain rate using (44). 5. Check convergence using the effective inelastic strain rate (47) where EN denotes the sum of effective inelastic strain rate over all the Gauss integration points and is defined as

with the effective strain rate EE =G, i. ti.. where Nc is the number of Gauss points in the solid phase. If error < tolerance, stop iteration. If error > tolerance, continue iteration (return to F. 1). End of iteration G. Check the number of iterations Z, at t = t, with respect to the number Z,_i at t = t,_, . If Z, > I,_, , decrease time step as At + 0.95 At. If z, 6 In-,, increase time step as At + 1.05 Ar, If the time limit (complete solidification) has not been reached, continue time stepping. End of time stepping. 4. Numerical results In the simulations reported here, four-noded bilinear quadrilateral elements have been employed and a temperature cooling boundary condition of the following form is used: T,(t) = T, + (rr, - T,) eWRf,

(48)

where T, is the final steady state temperature, T,,, is the melting temperature and R is a cooling rate parameter. The integration parameter y in (17)-(X) is 0.85. For the stress problem, a plane stress condition is assumed. The material used in the following examples is pure aluminium with thermal and mechanical properties listed in Tables 1 and 2, respectively. Table 1 Thermal properties

of aluminium

Heat conductivity in solid Heat conductivity in liquid Heat capacity in solid Heat capacity in liquid Latent heat Density Initial temperature Melting temperature

(After [ll])

KS 4 cs CL L P T, T,

0.0548 0.0548 0.2526 0.2526 94.44 2650 660 660

kcal/ms-“C kcallms-“C kcallkg-“C kcal/kg*“C kcal/ kg kg/m3 “C “C

348

N. Zabaras et al., Front tracking thermomechanical model

Table 2 Mechanical properties

of aluminium

(After (111)

Poisson’s ratio v = 0.37. Young’s modulus E(T)“. Temperature

“C

E(T) MPa C.

0

660

6.93 x lo4

4.0458 x lo4

Coefficients of constitutive Coefficients Values

d. Thermal expansion Temperature a(T) m/m*“C

“C

law (eq. (44)). A

B

c

n

0.382 x 10” set-’

0.037 lIMPa

18849 “K

3.84

coefficient a(T)‘. 25

300

400

660

23.19 x lo+

27.86 x 1O-6

30.23 x 1O-6

38.355 x 1o-6

a The variation of E(T) is linear within the temperature interval O-660 “C. b The variation of a(T) is assumed to be piecewise linear within the temperature intervals 25-300 “C, 3~-4~ “C and 400-660 “C.

4.1.

A unidirectional solidification example

To demonstrate the finite element model described previously, a ore-dimensional solidification problem is undertaken. The geometry of the problem is depicted in Fig. 3. The region of 0.01 m x H m is initially occupied with high purity liquid aluminium at the melting temperature. The boundary side x = 0 is assumed to be fixed in x direction and cooled with temperature T,,(t) given by (48). The boundaries at y = 0 and y = 0.01 m are insulated and restrained from motion in the y direction. The hydrostatic pressure applied on the freezing front at x = h(t) is given as p(t) = pg(H - h(t)), where H denotes the height of the free surface and g the gravity constant. The first case considered here is that H = 0.25 m, Ta = 300 “C and R = 0.023 set-’ and the second case is that H = 0.25 m, T, = 51.7 “C and R = 0.023 set-‘. In this unidirectional solidification example, only two finite elements were placed in the y direction, while in the x direction the simulation started with two elements and finally ended with 38 elements for the case T, = 300 “C and 46 elements for the case T, = 51.7 “C, respectively. Both cases have been studied by Heinlein et al. [ll] where the boundary element method is used for the solution of the heat transfer problem and a direct method is employed to calculate the stress field. The front position versus time is plotted in Fig. 4 for both cases. As expected, the front moves faster for the smaller value of 1-,. Figures 5 and 6 show the temperature histories for the same cases at four different locations. All these results are in perfect agreement with those reported by Heinlein et al. [ll]. The calculated stress histories at various locations for the two cases are shown in Figs. 7 and 8. The stress in y direction, u , is tensile, in general. As shown in these figures, aY builds up any location with the a&al of the front due to thermal contractions. However, the stress rate due to thermal contractions will eventually decrease

N. Zabaras et al., Front tracking thermomechanical

model

349

q=o

u,=(

I

Y

Fig. 3. Geometry

for unidirectional

solidification.

0.3

100

200

300

Time (set) Fig. 4. Front positions for the unidirectional

solidification

problem.

350

N. Zabaras et al., Front tracking thermomechanical

model

Q

0

200

100

300

Time (set) Fig. 5. Temperature

history at various locations for the unidirectional

solidification

problem;

case T, = 300°C.

800 -

x=O.lOm

600

0

300

200

100

Time (set) Fig. 6. Temperature

history at various locations for the unidirectional

solidification

problem;

case T. = 51.7 “C.

N. Zabaras et al., Front tracking thermomechanical

model

351

5 -

I

x=0.03287

m

I

200

100

300

Time (set) Fig. 7. Lateral stress history at various locations for the unidirectional

solidification

200

100

case T, = 300°C.

300

Time (set) Fig. 8. Lateral stress history at various locations for the unidirectional

problem;

solidification

problem;

case T, = 51.7 “C.

352

N. Zabaras et al., Front tracking ther~o~echanical

model

and inelastic relaxation will become dominant. For the case of Ta = 300 “C (Fig. 7), it is shown that the stress near x = 0 eventually relaxes. However, for the case of T, = 51.7 “C (Fig. 8), a near saturation of the stress occurs in the region close to x = 0. Further away from the boundary, x = 0, the stresses are smaller since the temperature drop there is much smaller than that at x = 0. Also with T, = 51.7 “C, the overall cooling rate is much higher than that in the case of T, = 300 “C. As a result, the residual stresses left in the material when solidification is complete are small for the case T, = 300 “C and much higher for the case T, = 51.7 “C. Note that for both cases, a, is calculated to be approximately equal to -p(t). The above stress results are similar in shape to those reported in [ll]. However, their magnitude is about 15% less. This is expected since the displacement boundary restrictions imposed by Heinlein et al. are close to plane strain rather than to plane stress. To demonstrate the effect of melt pressure, a higher pressure p(t) = pg(H - h(t)) with H = 25 m is assumed with T, = 300°C and the rest of the data remains the same. The distribution of the lateral stress is shown in Fig. 9. It can be concluded that high pressure reduces the overall stress distribution. The effect of high pressure can be very significant in the early stage of solidification of constrained castings when the fluid pressure forces the solid shell against the mold (not considered here), thus delaying the formation of air-gaps in the shell/mold interface. 4.2. Solidi~c~tio~ of a circular cylinder Solidi~~ation of an axially symmetric cylinder filled with liquid pure aIuminium at melting temperature is considered here. Circumferential variations are ignored and all variables are

200

100

Time (set) Fig. 9. Lateral stress history at various locations 0.649 MPa, T, = 300 “C.

for the unidirectional

solidifcation

problem;

case p,,,

=

N. Zabaras et al., Front tracking t~erm~mec~u~ica~ mod.4

Fig. 10. Geometry

353

for a solidifying circular section.

assumed to depend on the radius Y and time t. A quarter of the circular cylinder is studied as depicted in Fig. 10, where the boundary conditions for mechanical and thermal problems are also shown. For ali the cases discussed in this example, an initial mesh of six elements was used with one element in the Ydirection, while during simulations more elements were added in the r direction and by the time the whole circular section was solidified the finite element mesh was consisting of 96 elements. The present problem can be considered the same as the problem of analyzing a section of a continuously cast metal strand as it proceeds from a completely molten state in the mold through the water spray cooled region until eventually the section has completely solidified. Cylindrical casting experiments were performed by Nishida et al. [23] and commercial finite element software was employed by Smelser and Richmond [Zl] to numerically simulate these experiments. The above references cannot be compared directly with the methodology used in this paper since they allow for mold/shell interactions. However, since the material involved was the same as that used here, it is expected that for approximately similar cooling patterns and geometry, quatitatively, the same residual stresses must be obtained. The external boundary r = 0.018 m is cooled as given by (48) with T, = 500 “C and R = 0.10 set-‘. For the case of continuous casting, the melt pressure is given as p(t) = pg(N - z(t)), where z(t) is the z-coordinate of the interface at time t. However, since here dp/dt = -pg dz/dt can be considered very small, one can assume that p = constant. In this work, p = 0.001 MPa. Figures 11 and 12 show the front position and temperature histories at various locations, respectively. Figure 13 shows the stress distribution when the solidification is complete and Fig. 14 shows the stress distribution at t = 7.347 sec. The stress history at location r = 0.01739 m is given in Fig. 15. As seen from these figures, a compressive hoop stress is almost always present at the surface of the cylinder while near the center of the cylinder, tensile stresses are developed. The trend of the stresses found in this paper is similar to that reported by Smelser and Richmond [21].

N. Zabaras et al., Front tracking thermomechanical

354

0

2

4

6

m&e/

8

10

Time (set) Fig. 11. Radial front pasition for a solidifying circular section, T, = 500 “C, R = 0.10 SW-‘.

680

660

3 g

620

i? g

600

580

-

r=O.OlSOm

-

r=O,O135m r=0.0090m r=O.O045m

-

560 0

I

I

I

I

2

4

6

8

10

Time (SW) Fig. 12. Temperature

history at various locations for a solidifying circular section, T, = 500 “C, R = 0.10 set-‘.

N. Zabaras et al., Front tracking thermomechanical

v

-

-2 f

I 0.01

0.00

355

model

Radial stress Hoop stress

0.02

r Cm> Fig. 13. Stress distribution p = 0.001 MPa.

at the end of solidification,

t = 8.18sec,

for a solidifying

circular

section

0.01 r @-0 Fig. 14. Stress distribution

at time t = 7.347 set for a solidifying circular section with p = 0.001 MPa.

with

356

N. Zabaras et al., Front trucking t}lermome~h~n~cal model

0

2

4

6

8

10

Time (set) Fig. IS. Stress history at r = 0.01739 m for a soEdifying circular section with P = 0.001 MPa.

Finally, to demonstrate the effect of melt pressure, a case with the same conditions as before except p = 0.01 MPa, which is ten times larger, is considered. As shown in Fig. 16, the final stresses obtained in this case are similar to those in the case with smaller melt pressure, but the stresses at early time are quite different. Therefore, one can conclude that the melt pressure has a signi~cant effect on the stresses at early time. The surface radial displacement is also plotted in Fig. 17 as a function of time for the different values of the melt pressure. As expected, the cylinder initially expands due to the melt pressure, but finally it contracts.

Consider a square section of 0.04 m x 0.04 m occupied by pure aluminium liquid at melting temperature. Due to symmetry, one quarter of the square is considered. The temperature and mechanical boundary conditions are shown in Fig. 18. The boundary temperature T,(t) is defined from (48) with 7’a= 500 “C and R = 0.10 set-l. The melt pressure is assumed to be constant at all time with p(t) = 0,001 MPa. The finite element discretization started with 16 and ended with 144 finite elements. The propagation of the solid/liquid interface front and the temperature history at various locations are plotted in Figs. 19 and 20, respectively. The horizontal and vertical displacements of the boundary y = 0 are plotted in Fig. 21 at various times. The principle stresses at the center of each finite element are plotted in Figs. 22-24 at different times. As seen from these figures, compressive lateral stresses develop in the region close to the surface of the prism and tensile stresses occur inside. As solidification proceeds, the lateral stresses becomes

N. Zabaras et al., Front tracking thermomechanical

v * 0

2

4

model

357

Radial stress Hoop stress

6

8

10

Time (set) Fig. 16. Stress history at r = 0.01739 m for a solidifying circular section with P = 0.010 MPa.

0.02

p

2 .%

p=O.OlOh4Pa MPa

o.oa I ‘f

3 !A -0.02 E B ‘5 ?I .*

-0.04

2 p! -0.06

-0.08 0

2

4

6

8

10

Time (set) Fig. 17. Radial surface displacement

history for a solidifying circular section.

N. Zabaras et al., Front tracking thermomechanical

358

model

Y

uy=o q=o

Solid

q=o u,=o X

Fig. 18. Geometry

for the problem with a square cross section.

0.02

t=10.2487 set t= 9.8451 set

t= 8.8688 set

E

0.01

A

$%, ; ; ;

t=2.?96 seci t=1.6713 set t=0.6407 set

II II

0.00 If 0.01

0.00

Fig. 19. Front propagation

for solidification

0.02

x Cm> in a square region, T, = 500 “C, R = 0.10 set-‘.

N. Zabaras et al., Front tracking thermomechanical

359

model

620 600

580

* I

x=0.020 m, y=O.oOSm x=0.020m, y=O.OlOm Boundarytemperature I I I

0

2

560

540

x=0.005m, y=O.OOS m

4

\ -\ I

6

’

8

I 10

12

Time (see) Fig. 20. Temperature

history at various locations for solidification

in a square region, T, = 500 “C, R = 0.10 set-‘.

0.08

--o-

----fr-

-0.02 -

0.00 Fig. 21. Displacements of the boundary R = 0.10 set-‘, p = 0.001 MPa.

t= 0.0459sec

t=

5.0405set t=10.2488sec

I

- - - -

-

x-displacement y-displacement

0.01 x @I y = 0 at various time for solidification

0.02 in a square region,

T, = SOO”C,

N. Zabaras et al., Front tracking thermomechanical

360

1.4 MPa -

0.020

model compression tension

= -

,

- II _

tl

0.015- II _ II g

A

O.OlO- II c

VI

120

Fig. 22. Stress distributions p = 0.001 MPa.

at time t = 5.0405 set for solidification

in a square

1.4 MPa -

~~~:,,~~~:_:,.,

0.000

;,

0.005

0.010

compression tension

-

,;,

T, = 500 “C, R = 0.10 set-‘,

region,

;,

,;

0.015

0

x (m> Fig. 23. Stress distributions p = 0.001 MPa.

at time t = 8.2422 set for solidification

in a square

region,

T, = 500 “C, R = 0.10 set-‘,

361

320 Fig. 24. Stress distributions p = 0.001 MPa.

at time t = 10.2488 SW for solidification

in a square region, T, = 500 “C, R = 0.10 set-‘,

larger. The maximum stress appears at the center of the surface. At the end of solidification, large tensile stresses are also developed at the center of the prism. In the region close to the corner of the prism, stresses are generally small.

5. Cundusion In this report, a two-dimensional deforming finite element model is developed to analyze the temperature and stress fields in a solidifying body. For one-dimensional problems, the lateral stress was found to be always tensile and larger in the region close to the cooling boundary where solidification starts. For the cylindrical section, the residual hoop stresses were compressive in the region close to the outer surface and tensile in the region close to the center. The radial stress was tensile everywhere. Finally, the melt pressure was shown to significantly alter the deformation at the early stage of solidification. A similar residual stress pattern was also obtained for the square prism. Some of the advantages of the present model are easy implementation and its ability to accurately account for the solid/liquid interface conditions. It was shown that by using the freezing interface velocity and location as primary variables of the analysis and by involving a rate form of the principle of virtual work, one can easily account for the initial hydrostatic stress state of a material particle immediately upon solidification as well as for the Stefan condition on the freezing front. In the so called fixed domain techniques, where the problem is

N. Zabaras et al., Front tracking thermomechanical

362

model

formulated on a fixed grid without an explicit reference to the freezing interface, one cannot account directly for thermal or mechanical interface conditions. For example, when the enthalpy method is used to solve the heat transfer part of a solidification problem, the freezing interface energy balance is incorporated into the governing equations via the definition of appropriate source terms. However, no attempt has yet been made to extend these ideas to the deformation analysis of phase change problems and to account for the hypoelastic stress condition on the freezing front. Further research has been undertaken to extend the present work to general axially symmetric geometries. Problems where the solid is initially in contact with a mold will be considered as well. It is important to examine the effect of melt pressure on the formation of air-gaps in the solid shell/mold boundary. Inelastic constitutive models with state variables need to be implemented. Such models will not only predict deformation, but also the evolution of microstructure during solidification.

Acknowledgment The heat transfer part of this work has been supported by NSF grant CBT-8802069 to the University of Minnesota. The computing has been supported with a grant from the Minnesota Supercomputer Institute.

Appendix A. Matrix [A] and vector {H} in eq. (26) In this appendix,

the components

of matrix [A] and the vector {H} in (26) are given, that

is,

(A4 Am = fikSih

Hi2 = a, ,

(no sum on i) ,

(A.2)

(A.41

where (N,b, Nb,) are the directional cosines of the unit normal Nb to the segment b, (tix, tiy) are the directional cosines of the unit tangent vector t to 130, at the interface node i. 8, is the Kronecker delta, i.e. ai, = 1 when i = h and 6;, = 0 when i # h. The shape functions !J’t and Vi are the linear shape functions defined over the segment b and a, and a2 are zero in this work. As an example, the matrix [A] and the vector {H} for one linear interface segment b with nodes 1 and 2 (Fig. 25) are given as

N. Zabaraf et al., Front tracking thermomechanical

model

363

i2-L?J!+

1

Fig. 25. The configuration

of a solid/liquid

freezing interface segment

-I a

If

anh [&VT,

- K,VT,] a1

I

anb [I&VT, - K,VT,] 1

. NbTiqdr'

64.6)

.

Nb!P; dT

%

Note that all variables involved in (A.l)-(A.6)

are calculated

at time L

References 111H.S. Carslaw and J.C Jaeger, Conduction of Heat in Solids, 2nd Ed. (Oxford University Press, Oxford, 1959). PI J. Crank, Free and Moving Boundary Value Problems (Clarendon Press, Oxford, 1984). [31 D.R. Lynch and K. O’Neill, Continuously deforming finite elements for the solution of parabolic problems with and without phase change, Internat. J. Numer. Methods Engrg. 17 (1981) 81-96. [41 R.D. Lynch, Unified approach to simulation on deforming elements with application to phase change probems, J. Comput. Phys. 47 (1982) 387-411. PI N. Zabaras and Y. Ruan, Deforming finite element method analysis of inverse Stefan problems, Internat. J. Numer. Methods Engrg. 28 (1989) 295-313. N. Zabaras and Y. Ruan, Moving finite element simulation of two dimensional Stefan problems, University of PI Minnesota Supercomputer Institute Research Report, UMSI 89155, 1989. [71 K. O’Neill, Boundary integral equation for moving boundary phase change problems, Internat. J. Numer. Methods Engrg. 19 (1983) 1825-1850. 181 N. Zabaras and S. Mukherjee, An analysis of solidification problems by the boundary element method, Internat. J. Numer. Methods Engrg. 24 (1987) 1879-1900. [91 0. Richmod and R. Tien, Theory of thermal stress and air-gap formation during the early stages of solidification in a rectangular mold, J. Mech. Phys. Solids 19 (1971) 273-284. [lOI R. Tien and 0. Richmond, Theory of maximum tensile stresses in solidifying shell of a constrained rectangular casting, ASME J. Appl. Mech. 49 (1982) 481-486.

364

[ll]

N. Zabaras et al., Front tracking thermomechanical

model

M. Heinlein, S. Mukherjee and 0. Richmond, A boundary element method analysis of temperature fields and stresses during solidification, Acta Mech. 59 (1986) 59-81. [12] L. Anand, Constitutive equations for the rate-dependent deformation of metals at elevated temperatures, ASME J. Engrg. Mater. Techn. 104 (1982) 12-17. (131 E.W. Hart, Constitutive relations for the non-elastic deformation of metals, ASME J. Engrg. Mater. Techn. 98 (1976) 193-202. [14] J.O. Kristiansson, Thermal stresses in the early stage of solidification, J. Thermal Stresses 5 (1982) 315-330. [15] J.H. Weiner and B.A. Boley, Elasto-plastic thermal stresses in a solidifying body, J. Mech. Phys. Solids 11 (1963) 145-154. (16) T. Inoue and Z.G. Wang, High temperature behavior of steels with phase transformation and the simulation of quenching and welding processes, in: J. Carlsson and N.G. Ohlson, eds., Mechanical Behavior of Materials 14 (Pergamon Press, Oxford, 1984) 1005-1010. [17] P. Perzyna, Fundamental problems in viscoplasticity, Adv. Appl. Mech. 9 (1966) 243-377. [18] B.G. Thomas, I.V. Samarasekera and J.K. Brimacombe, Mathematical model of the thermal processing of steel ingots: Part I. Heat flow model, Metallurgical Transactions B 18B (1987) 119-130. 1191 B.G. Thomas, I.V. Samarasekera and J.K. B~ma~ombe, Mathematical model of the thermal processing of steel ingots: Part II. Stress model, Metallurgical Transactions B 183 (1987) 131-147. [20] P.J. Wray, Modelling of casting and welding processes, AIME Conference Proceedings (1980) 245-257. [21] R.E. Smelser and 0. Richmond, Constitutive model effects on stresses and deformations in a solidifying circular cylinder, Proceedings of the Fourth Engineering Foundation Conference On Modeling of Casting and Welding Processes, Palm Coast, Fl (1988) 17-22. f22] L.G. Hector Jr. and 0. Richmond, A thermome~hanical model of air gap nucleation during nonuniform directional solidification, Alcoa Laboratories Technical Report, No. 52-88-11, 1988. [23] Y. Nishida, W. Droste and S. Engler, The air-gap formation process at the casting-mold interface and the heat transfer mechanism through the gap, Metallurgical Transactions B 17B (1986) 833-844. [24] T.J.R. Hughes, Unconditionally stable algorithm for nonlinear heat conduction, Comput. Methods Appl. Mech. Engrg. 10 (1977) 135-139. [25] B.A. Boley and J.H. Weiner, Theory of Thermal Stresses (Wiley, New York, 1960). [26] VM. Sample and L.A. Lalli, Effects of thermomechanical history of hardness of aluminium, Mater. Sci. Techn. 3 (1987) 28-35. [27] R. Bodner and Y. Partom, Constitutive equations for elastic viscoplastic strain hardening materials, ASME J. Appl. Mech. 42 (1975) 385-389. [28] K. Morgan, R.W. Lewis and K.N. Seetharamu, Modelling heat flow and thermal stress in ingot casting, Simulation 36 (1981) 55-63. f29] 0. Richmond, Models of Stresses and Deformations in Solidifying Bodies, Modeling of Casting and Welding Processes (ASME, New York, 1982) 215-222. [30] J.R. Williams, R.W. Lewis and K. Morgan, An elasto-viscoplastic thermal stress model with applications to the continuous casting of metals, Internat. J. Numer. Methods Engrg. 14 (1979) l-9. [31] O.C. Zienkiewicz and I.C. Cormeau, Viscoplasticity-plasticity and creep in elastic solids-A unified numerical solution approach, Internat. J. Numer. Methods Engrg. 8 (1974) 821-846.