Michel Bellet

Brian G Thomas

Ch. 27 Solidification Macroprocesses (Thermal-mechanical Modeling of Stress, Distortion and Hot Tearing)

Keywords:

INTRODUCTION

The application of numerical methods to the mechanical modeling in solidification analysis has received a continuously growing interest over the last twenty years. After having concentrated their efforts on the thermal and microstructural predictions, research teams have been more and more interested in the coupled thermomechanical analysis. These developments have been motivated by the efforts done by the casting industry to increase the quality of final products while lowering the costs. Manufacturers are then interested in the development of new numerical tools able to model the thermomechanical response of castings during the processes. Accurate calculation of stress and distortions during casting is just the first step, however, as engineers are more interested in their practical consequences. These include residual stress and distortion, and defects such as segregation and the formation of cracks such as hot tears. As computing power and software tools advance, it is becoming increasingly possible to perform useful mechanical analysis of castings and these important related behaviors.

Despite the considerable progress in computational structural mechanics, and in the capabilities of intensive computing, the thermomechanical analysis of castings is still a challenge nowadays, for the following reasons:

• Many interacting physical phenomena are involved in stress-strain formation. Stress arises primarily from the mismatch of strains caused by the steep temperature gradients of solidification, and depends on the time-and microstructure-dependent inelastic flow of the material.

• Predicting distortions and residual stresses in cast products means being able to describe and calculate the history of the cast product and its environment on huge temperature intervals. This makes the mechanical problem highly non-linear, involving liquid-solid interaction and quite complex constitutive equations. Also the identification of reliable values of the numerous parameters involved in those relations is a very difficult task.

• The coupling between the thermal and the mechanical problems is an additional difficulty. This coupling comes from the mechanical interaction between the casting and the mold components, through gap formation or the build-up of contact pressure, modifying locally the heat exchange. This adds some complexity to the non-linear heat transfer resolution.

• As it does not make sense to perform such analyses without accounting for the presence of molds and their interaction with the castings, the problem to be solved is multidomain, often involving numerous deformable and interacting components.

• Shapes of cast parts are by essence complex, which first brings out frequent difficulties regarding the interface between CAD design and the mechanical solvers, and second demands great computational resources.

• In the case of solidification of semi-finished products, by continuous casting, the computational demand is also outstanding, because of the characteristic dimensions to be considered: meters, tens of meters.

This chapter summarizes some of the issues and approaches in performing computational analyses of mechanical behavior, distortion and hot-tearing during solidification. The governing equations are presented first, followed by a brief description of the methods used to solve them.

Finally, a few examples of recent applications in shape castings and continuous casting are introduced.

CONSTITUTIVE MODELS FOR METALLIC ALLOYS IN LIQUID, MUSHY AND SOLID STATE

The modeling of mechanical behavior requires solution of the equilibrium equations (relating force and stress), constitutive equations (relating stress and strain) and compatibility equations (relating strain and displacement). In casting analysis, the cast material may be in the liquid, mushy or solid state. Therefore, mechanical modeling of casting processes has to consider constitutive models for each of these states.

Liquid state: Newtonian model

Metallic alloys are generally considered as Newtonian fluids. Including thermal dilatation effects, the constitutive equation can be expressed as follows. in which the strain rate tensor ε is split into a mechanical part, showing the linear relation between the strain rate tensor and the stress deviator s, and a thermal part. In this equation, l µ is the dynamic viscosity of the liquid alloy, ρ is the density, and I is the identity tensor. Taking the trace of this expression,

v ε ⋅ ∇ = tr
, the mass conservation equation is recovered:

0 ) ( d d = ⋅ ∇ + ∂ ∂ = ⋅ ∇ + v v ρ ρ ρ ρ t t (27 .2)
In casting processes, the liquid flow may be turbulent, even after mold filling. This may occur because of buoyancy forces or forced convection like in jets coming out of the nozzle outlets in continuous casting processes. The most accurate approach, direct numerical simulation, is generally not feasible for industrial processes, owing to their complex shaped domains and high turbulence. To compute just the large-scale flow features, turbulence models are used, that increase the liquid viscosity according to different models of the small-scale phenomena. These models include the simple "mixing length" models, the two-equation models such as k-ε, and large eddy simulation (LES) models, which have been compared with each other and with measurements of continuous casting. [START_REF] Thomas | Comparison of Four Methods to Evaluate Fluid Velocities in a Continuous Casting Mold[END_REF][START_REF] Yuan | Study of Computational Issues in Simulation of Transient Flow in Continuous Casting[END_REF][START_REF] Yuan | Computational and Experimental Study of Turbulent Flow in a 0.4-Scale Water Model of a Continuous Steel Caster[END_REF]

Mushy state: Non-Newtonian model

Metallic alloys in the mushy state are very complex two-phase liquid-solid media. Their mechanical response is highly dependent on the local microstructural evolution, which involves several complex physical phenomena. To overcome this difficulty, the mushy state may be considered in a first approach as a single continuum. The mushy material is then modeled as a non-Newtonian fluid, according to the following equations:

         - = = + = - I ε s ε ε ε ε t K th m eq

Solid State: Elastic-Visco-Plastic Models

In the solid state, metallic alloys can be modeled either as elastic-plastic or elastic-viscoplastic materials. In this latest class of models, one of the simpler is expressed as follows, but it should be mentioned that a lot of models of different complexity can be found in the literature. [START_REF] Estrin | A versatile unified constitutive model based on dislocation density evolution[END_REF], 5

             - = - =       ∂ ∂ -       + ∂ ∂ + - + = + + = I ε s ε I σ σ I σ σ ε ε ε ε ε t K E T T E T T E E th m eq eq in el th in el d d 3 1 2 3 ) tr( 1 ) tr( 1 / 1 0 ρ ρ σ σ σ ν ν ν ν (27 .5)
The strain rate tensor ε is split into an elastic component, an inelastic (non reversible) component, and a thermal component. Equation (27.5b) yields the hypoelastic Hooke's law, where E is Young's modulus, ν the Poisson's coefficient, and σ a time derivative of the stress tensor σ . Equation (27.5c) gives the relation between the inelastic strain rate tensor in ε and the stress deviator, s, in which 0 σ denotes the scalar static yield stress, below which no inelastic deformation occurs (the expression between brackets is reduced to zero when negative). In these equations, the temperature dependency of all the involved variables should be considered. The effect of strain hardening may appear in such a model by the increase of the static yield stress 0 σ and the plastic consistency K with the accumulated inelastic strain eq ε , or with another state variable that is representative of the material structure. The corresponding scalar equation relating stress and inelastic strain rate von Mises invariants is:

0 ( ) m eq eq K σ σ ε = + (27 .6)
Inserting this into Eq. ( 27.5c) simplifies it to:

s ε s ε eq eq in eq eq in σ ε σ ε 2 d 3 d form, l incrementa in or, , 2 3 = = (27 .7)
Although metallic alloys show a significant strain rate sensitivity at high temperature, they are often modelled in the literature using elastic-plastic models, neglecting this important effect.

In this case, Eq. (27.7) still holds, but the flow stress is then independent of the strain rate. It may depend on the accumulated plastic strain because of strain hardening.

Implementation Issues

As stresses and distortions are generated mainly in the solid phase, the mechanical modeling of a cast part may be restricted to its solidified region at any instant. However, this approach, often used in the literature, has several drawbacks. First, doing so, the volumetric shrinkage that affects the mushy zone cannot be taken into account. Second, the liquid and mushy regions may alter the distortion and stress of the solidified regions. Finally, the mechanical behavior of the mushy zone is of crucial importance to quality problems. Thus, it is preferable to consider the entire casting, including the mushy and liquid regions.

To avoid managing different constitutive equations for the different states, numerous authors have modeled the liquid and mushy state behavior simply by lowering the value of the Young modulus and taking Poisson coefficient close to ½, using then a single elastic-viscoplastic (or elastic-plastic) constitutive model for the entire casting [START_REF] Weiner | Elasto-plastic thermal stresses in a solidifying body[END_REF][START_REF] Smelser | Constitutive model effects on stresses and deformations in a solidifying circular cylinder, Modeling of Casting and Welding Processes IV[END_REF][START_REF] Bellet | Thermomechanics of the cooling stage in casting processes: 3D finite element analysis and experimental validation[END_REF][START_REF] Fjaer | ALSPEN -A mathematical model for thermal stresses in DC-cast Al billets[END_REF][START_REF] Wiese | Modelling stress development during the solidification of grey iron castings[END_REF][START_REF] Inoue | Simulation of solidification and viscoplastic stresses during vertical semi-continuous direct chill casting of aluminium alloy[END_REF][START_REF] Moitra | Thermo-mechanical model of Steel Shell Behavior in the Continuous Casting Mold[END_REF][START_REF] Celentano | A finite element model for thermomechanical analysis in casting processes[END_REF] However, this approach may suffer from numerical difficulty, and furthermore fails to account accurately for the significant thermal dilatation and shrinkage that affect the mushy and liquid regions.

To overcome these difficulties, other authors use a different constitutive equation for the different physical states of the alloy. The whole casting is modeled, and the constitutive equation is chosen according to the local state. One implementation [START_REF] Li | Thermo-Mechanical Finite-Element Model of Shell Behavior in Continuous Casting of Steel[END_REF] simply changes the constants in Eq. In each of these previous mechanical models, the liquid, mushy, and solid zones are considered as a single continuum. The velocity of the liquid phase is not distinguished from the velocity of the solid phase, and the individual dendrites and grain boundaries are not resolved, so fluid feeding, porosity formation, and hot tearing are clearly oversimplified. In the context of stress-strain prediction, this approximation seems valid. To predict defect formation from first principles, however, might require a multiphase approach, involving the different phases, including liquid, solid, and even gas. Recent approaches, like those developed by Nicolli et al [START_REF] Nicolli | Modeling of macrosegregation caused by volumetric deformation in a coherent mushy zone[END_REF] and Fachinotti et al [START_REF] Fachinotti | Two-phase thermo-mechanical and macrosegregation modelling of binary alloys solidification with emphasis on the secondary cooling stage of steel slab continuous casting[END_REF] take into account the mechanical interaction between a deformable solid skeleton and the liquid phase, the momentum transfer between the two phases being expressed by a Darcy law.

Example of Constitutive Equations

Material property data are needed for the specific alloy being modeled and in a form suitable for the constitutive equations just discussed. This presents a significant challenge for quantitative mechanical analysis, because measurements are not presented in this form, and only rarely supply enough information on the conditions to allow transformation to an alternate form. As an example, the following elastic-visco-plastic constitutive equation was developed for the austenite phase of steel by Kozlowski et al [START_REF] Kozlowski | Simple Constitutive Equations for Steel at High Temperature[END_REF] by fitting constant strain-rate tensile tests from Wray

m eq C eq C f eq eq f T f C C T f T m T T ε σ σ σ ε σ σ - - -   × = - -     = × + × + × = - × = - + × = - × (27 .8)
This equation, and a similar one for delta-ferrite, have been implemented into the finite-element codes CON2D [START_REF] Li | Thermo-Mechanical Finite-Element Model of Shell Behavior in Continuous Casting of Steel[END_REF] and THERCAST [START_REF] Costes | 3D thermomechanical simulation of the secondary cooling zone of steel continuous casting[END_REF] and applied to investigate several problems involving mechanical behavior during continuous casting.

Elastic modulus is a crucial property that decreases with increasing temperature. It is difficult to measure at the high temperatures important to casting, owing to the susceptibility of the material to creep and thermal strain during a standard tensile test, which results in excessively low values. Higher values are obtained from high-strain-rate tests, such as ultrasonic measurements. [START_REF] Donsbach | Ultrasonic Measurement of Elastic Constants at Temperatures from 20 to 1100 oC[END_REF] Elastic modulus measurements in steels near the solidus temperature range from ~1 GPa 25 to 44 GPa. [START_REF] Hub | Measurement of Velocity and Attenuation of Sound in Iron up to the Melting Point[END_REF] Typical modulus data by Mizukami et al [START_REF] Mizukami | Mechanical Properties of Continuously Cast Steels at High Temperatures[END_REF] include values ~10 GPa near the solidus and have been used in previous analyses. [START_REF] Li | Maximum Casting Speed for Continuous Cast Steel Billets Based on Sub-Mold Bulging Computation[END_REF][START_REF] Li | Thermo-Mechanical Finite Element Model of Bulging and Hot Tearing During Continuous Casting of Steel Billets[END_REF] The density needed to compute thermal strain in Eqs. (27.1), (27.4) or (27.5d) can be found from a weighted average of the values of the different solid and liquid phases, based on the local phase fractions. For the example of plain low carbon steel, the following equations were compiled [START_REF] Li | Thermo-Mechanical Finite-Element Model of Shell Behavior in Continuous Casting of Steel[END_REF] based on the solid data for ferrite (α), austenite (γ), and delta (δ) from Harste et al [START_REF] Harste | Shrinkage and Formation of Mechanical Stresses during Solidification of Round Steel Strands[END_REF][START_REF] Harste | Investigation of the shinkage and the origin of mechanical tension during the solidification and successive cooling of cylindrical bars of Fe-C alloys[END_REF] and the liquid (l) measurements from Jimbo and Cramb 32 :

l l o o o o o l kg m f f f f T C T C T C C C T C C C C C T C α α γ γ δ δ α γ δ ρ ρ ρ ρ ρ ρ ρ ρ ρ - = + + + = - -× - = - + - = - + = - - - - (27 .9)
Specialized experiments to measure mechanical properties for use in computational models will be an important trend for future research in this field.

MODEL DEVELOPMENT

Thermo-mechanical Coupling

Coupling between the thermal and mechanical analyses arises from several sources. First, regarding the mechanical problem, besides the strain rate due to thermal expansion and solidification shrinkage, the material parameters of the preceding constitutive equations strongly depend on temperature and phase fractions, as shown in the previous section. Second, in the heat transfer problem, the thermal exchange between the casting and the mold strongly depends on local conditions such as the contact pressure or the presence of a gap between them (as a result of thermal expansion and solidification shrinkage). This is explained in the next two paragraphs.

Air gap formation: conductive-radiative modeling

In the presence of a gap between the casting and the mold, resulting from their relative deformation, the heat transfer results from concurrent conduction through the gas within the gap and from radiation. The exchanged thermal flux, q gap , can then be written: Therefore, the conductive heat exchange coefficient g k h gas cond = should be limited by a finite value 0 h , corresponding to the "no-gap" situation, and depends on the roughness of the casting surface. A recent benchmark exercise has demonstrated the significant impact of the consideration of gap formation on temperature prediction in shape casting. [START_REF] Kron | Comparison of numerical simulation models for predicting temperature in solidification analysis with reference to air gap formation[END_REF] Specific examples of these gap heat transfer laws are provided elsewhere for continuous casting with oil lubrication, [START_REF] Kelly | Initial Development of Thermal and Stress Fields in Continuously Cast Steel Billets[END_REF] and continuous casting with mold flux. [START_REF] Meng | Heat Transfer and Solidification Model of Continuous Slab Casting: CON1D[END_REF]

1 1 1 ) ( ) (

Effective contact : heat transfer as a function of contact pressure

In the case of an effective contact, the conductive heat flux increases with the contact pressure according to a power law. [START_REF] Madhusudana | Contact Heat Transfer -The last decade[END_REF] Still denoting 0 h the heat exchange coefficient corresponding to no gap and no contact pressure, we can write, in order to ensure the continuity with the gap situation:

) )( ( 0 m c B c contact T T Ap h q - + = (27. 1 1)
with c p the contact pressure, A and B two parameters which depend on the materials, the presence of coating or lubricating agent, the surface roughness, and the temperature. The parameters and possibly the laws governing their evolution need to be determined experimentally.

Numerical solution

We focus here on the resolution of the momentum conservation, from which distortions and stresses can be calculated. The energy conservation is not discussed in this paper.

Momentum conservation

At any time, and in any location of the solidifying material, whatever its current state (liquid, k gas mushy or solid) the momentum conservation is expressed by:

0 = - + ∇ - ⋅ ∇ = - + ⋅ ∇ γ g s γ g σ ρ ρ ρ ρ p (27.1 2)
where g denotes the gravity and γ γ γ γ the acceleration. The acceleration is actually noticeable only in liquid pools, when they are affected by convection.

The weak form of the preceding equation can be obtained. Keeping the velocity and pressure as primitive unknown variables, it is written as [START_REF] Rappaz | Numerical Modeling in Materials Science and Engineering[END_REF] :

       = ∀ = ⋅ + ⋅ - ⋅ - ⋅ ∇ - ∀ ∫ ∫ ∫ ∫ ∫ ∫ Ω Ω Ω Ω ∂ Ω Ω 0 tr 0 d d : ) ( * * * * * * * * dV p p dV t dV dS dV p dV in ε v v v g v T v ε v s v ρ ρ (27.1 3)
where T is the external stress vector, g the gravity. In such a velocity-pressure formulation, the second equation is a weak form of the incompressibility of inelastic deformations.

Finite element formulation and numerical implementation

The thermal-mechanical modeling equations just presented must be solved numerically, owing to the complex shape of the casting process domain. Although finite-difference approaches are popular for fluid flow and solidification analysis to compute the temperature field, the finite element formulation is preferred for the mechanical analysis, owing to its historical advantages with unstructured meshes and accurate implicit solution of the resulting simultaneous algebraic equations.

Spatial discretization using finite elements, combined with time discretization using finite differences, yield a set of non-linear equations

0 ) , ( = P V R mech
, in which the unknowns are the velocity components and the pressure value at each node of the finite element mesh. In some formulations, the incompressibility condition is directly included in the momentum equation using a penalty technique to enforce it. This results in a velocity formulation 0 ) ( = V R mech in which the velocity field is the only unknown, but which is known to give rise to pressure locking problems when the incompressible inelastic strains become too large. [START_REF] Zienkiewicz | The Finite Element Method[END_REF] In the preceding formulations, nodal velocities can be replaced by nodal displacements, as they are linked by the time integration scheme.

Solving the previous global finite element non-linear system 0 = mech R can be a daunting task that is subject to significant convergence difficulties. Generally, it is solved using a full or modified Newton-Raphson method, [START_REF]ABAQUS Theory Manual[END_REF] which iterates to minimize the residual error in Eq. (27.13). These methods require local consistent tangent operators

ε s ∂ ∂ and ε ∂ ∂p in order to form a global consistent stiffness matrix V R ∂ ∂ mech
. [START_REF] Zienkiewicz | The Finite Element Method[END_REF] At the local level, when the constitutive equations involve strain rate dependency, then an algorithm is also required to integrate the constitutive equations to provide an estimate of ) (ε s and ) (ε p and the previously mentioned tangent operators. When the constitutive equations are highly non-linear, it is very useful to employ an implicit algorithm to provide better estimates of inelastic strain at the local level. Many methods have been developed [START_REF] Nemat-Nasser | An explicit Algorithm for large-strain, large-strain rate elastic-viscoplasticity[END_REF][START_REF] Lush | An Implicit Time-integration Procedure for A Set of Internal Variable Constitutive Equations for Isotropic Elastic-Viscoplasticity[END_REF][START_REF] Zhu | Coupled thermal-mechanical finite-element model with application to initial solidification[END_REF] which require solving two or more ordinary differential equations at each local integration point. The bounded Newton-Raphson method developed by Lush et al [START_REF] Lush | An Implicit Time-integration Procedure for A Set of Internal Variable Constitutive Equations for Isotropic Elastic-Viscoplasticity[END_REF] and later improved upon by Zhu [START_REF] Zhu | Coupled thermal-mechanical finite-element model with application to initial solidification[END_REF] was implemented into a user-subroutine in ABAQUS and found to improve greatly accelerate the solution. [START_REF] Koric | Efficient Thermo-Mechanical Model for Solidification Processes[END_REF] Alternatively, an operator splitting method can be used to march through time by alternating between the global and local levels without iteration at either level. [START_REF] Li | Thermo-Mechanical Finite-Element Model of Shell Behavior in Continuous Casting of Steel[END_REF][START_REF] Koric | Efficient Thermo-Mechanical Model for Solidification Processes[END_REF]

Boundary conditions: modeling of contact conditions. Multidomain approaches

At the interface between the solidifying material and the mold, a unilateral contact condition (i.e. including contact release) generally applies:

     = ⋅ ≥ ≤ ⋅ 0 ) ( 0 0 g g n σn n σn (27.1 4)
where g is the local interface gap width (positive when air gap exists effectively) and n is the local outward unit normal to the part. The fulfillment of Eq. ( 27.14) can be obtained by means of a penalty condition, which consists in applying a normal stress vector T proportional to the normal distance or the normal velocity difference (or a combination of both) via a penalty constant χ p :

n n v v σn T ⋅ - - = = ) ( mold p χ (27.1 5)
Different methods of local adaptation of the penalty coefficient χ p have been developed, among which the augmented Lagrangian method. [START_REF] Glowinski | Augmented Lagrangian and Operator-Splitting Methods in Non-linear Mechanics[END_REF] More complex and computationally expensive methods, such as the use of Lagrange multipliers may also be used. [START_REF] Wriggers | On contact between three-dimensional beams unergoing large deflections[END_REF] The possible tangential friction effects between part and mold can be taken into account by a friction law, such as a Coulomb model for instance. In this case, the previous stress vector has a tangential component, T τ , given by:

) ( 1 mold mold c f p v v v v T - - - = µ τ (27.1 6)
where

n n ⋅ = - = σ σ n c p
is the contact pressure, and µ f the friction coefficient.

The previous approach can be extended to the multidomain context in order to account for the deformation of mold components. The local stress vectors calculated by Eq. ( 27.15) can be applied onto the surface of the mold, contributing then to its deformation. For most casting processes, the mechanical interaction between the cast product and the mold is sufficiently slow (i.e. its characteristic time remains significant with respect to the process time) to permit a staggered scheme within each time increment: the mechanical problem is successively solved in the cast product and in the different mold components. A global updating of the different configurations is then performed at the end of the time increment. This simple approach gives access to a prediction of the local air gap size g, or alternatively of the local contact pressure p c , that are used in the expressions of the heat transfer coefficient, according to Eqs. (27.10) and (27.11). [START_REF] Jaouen | A numerical mechanical coupling algorithm for deformable bodies: application to part/mold interaction in casting process[END_REF]

3 Treatment of Regions in the Solid, Mushy and Liquid States

Solidified regions: Lagrangian formulation

In casting processes, the solidified regions generally encounter small deformations. It is thus natural to embed the finite element domain into the material, with each node of the computational grid corresponding with the same solid particle during its displacement. The boundary of the mesh corresponds then to the surface of the casting. This method, called Lagrangian formulation, provides the best accuracy when computing the gap forming between the solidified material and the mold. It is also the more reliable and convenient method for time integration of highly non-linear constitutive equations, such as elastic-(visco)-plastic laws presented in section 27.2.3.

Mushy and liquid regions: ALE modeling

When the mushy and liquid regions are modeled in the same domain as the solid (cf. discussion in section 27.2.4), they are often subjected to large displacements and strains arising from solidification shrinkage, buoyancy, or forced convection. Similar difficulties are generated in casting processes such as squeeze casting, where the entire domain is highly deformed. In these cases, a Lagrangian formulation would demand frequent remeshings in order to avoid mesh degeneracy, which is both computationally costly, and detrimental to the accuracy of the modeling. It is then preferable to use a so-called arbitrary Lagrangian Eulerian formulation (ALE). In an Eulerian formulation, material moves through the computational grid, which remains stationary in the "laboratory" frame of reference. In the ALE formulation, the updating of the mesh is partially independent of the velocity of the material particles in order to maintain the quality of the computational grid. Several methods can be used, including the popular "barycentering" technique which keeps each node at the geometrical centroid of a set of its neighbors. This method involves significant extra complexity to account for the advection of material through the domain, and the state variables such as temperature and inelastic strain must be updated according to the relative velocity between the mesh and the particles. In doing this, some surface constraints must be enforced in order to ensure mass conservation, expressing that the fluxes of mesh velocity and of fluid particle velocity through the surface of the mesh should remain identical. A review on the ALE method in solidification modeling is available, together with some details on its application. [START_REF] Bellet | ALE method for solidification modelling[END_REF]

Thermomechanical coupling

Because of the interdependency between the thermal and mechanical analyses, as presented in section above, their coupling should be taken into account all during the cooling process. In practice, the cooling time is decomposed into time increments, each increment requiring the solution of two problems: the energy conservation and the momentum conservation. With the highly nonlinear elastic-visco-plastic constitutive equations typical of solidifying metals, the incremental steps required for the mechanical analysis to converge are generally much smaller than those for the thermal analysis. Thus, these two analyses are generally performed in succession and only once per time increment. However, in the case of very rapid cooling, these solutions might be preferably performed together (including thermal and mechanical unknowns in a single set of non-linear equations), or else separate but iteratively until convergence at each time increment, otherwise the time step has to be dramatically reduced.

Hot tearing analysis

Hot tearing is one of the most important consequences of stress during solidification. Hot tearing is caused by a combination of tensile stress and metallurgical embrittlement. It occurs at temperatures near the solidus when strain concentrates within the interdendritic liquid films, causing separation of the dendrites and intergranular cracks at very small strains (on the order of 1 percent). This complex phenomenon depends on the ability of liquid to flow through the dendritic structure to feed the volumetric shrinkage, the strength of the surrounding dendritic skeleton, the grain size and shape, the nucleation of supersaturated gas into pores or crack surfaces, the segregation of solute impurities, and the formation of interfering solid precipitates.

The subsequent refilling of hot tears with segregated liquid alloy can cause internal defects that are just as serious as exposed surface cracks. The hot tearing of aluminum alloys is reviewed elsewhere. [START_REF] Eskin | Mechanical properties in the semi-solid state and hot tearing of aluminum alloys[END_REF] Hot tearing phenomena are too complex and insufficiently understood to model in detail, so several different criteria have been developed to predict hot tears from the results of a thermal-mechanical analysis.

Casting conditions that produce faster solidification and alloys with wider freezing ranges are more prone to hot tears. Thus, many criteria are solely based on thermal analysis. That of Clyne and Davies 50 simply compares the local time spent between two critical solid fractions 1 s g and 2 s g (typically 0.9 and 0.99, respectively), with the total local solidification time (or a reference solidification time). The "hot cracking susceptibility" is defined as: Criteria based on classical mechanics often assume cracks will form when a critical stress is exceeded, and they are popular for predicting cracks at lower temperatures [START_REF] Kinoshita | Thermal Elasto-plastic Stress Analysis of Solidifying Shell in Continuous Casting Mold[END_REF][START_REF] Kristiansson | Thermal stresses in the early stage of solidification of steel[END_REF][START_REF] Thomas | Mathematical Model of the Thermal Processing of Steel Ingots, Part II: Stress Model[END_REF][START_REF] Okamura | Calculation of Bulging Strain and its Application to Prediction of Internal Cracks in Continuously Cast Slabs[END_REF] . This critical stress depends greatly on the local temperature and strain rate. Its accuracy relies on measurements, such as the submerged split-chill tensile test for hot tearing. [START_REF] Ackermann | In Situ Tensile Testing of Solidifying Aluminum and Al-Mg Shells[END_REF][START_REF] Bernhard | Simulation of shell strength properties by the SSCT test[END_REF][START_REF] Suzuki | In-Situ Measurement of Tensile Strength of Solidifying Steel Shells to Predict Upper Limit of Casting Speed in Continuous Caster with Oscillating Mold[END_REF] Measurements often correlate hot tear formation with the accumulation of a critical level of mechanical strain while applying tensile loading within a critical solid fraction where liquid feeding is difficult. This has formed the basis for many hot-tearing criteria. That of Yamanaka et al. [START_REF] Yamanaka | Measurement of Critical Strain for Solidification Cracking, Modelling of Casting, Welding, and Advanced Solidification Processes[END_REF] accumulates inelastic deformation over a brittleness temperature range, which is defined, for in which the critical strain cr ε is 1.6% at a typical strain rate of 4 10 3 -× s -1 . Careful measurements during bending of solidifying steel ingots have revealed critical strains ranging from 1 to 3.8%. [START_REF] Yamanaka | Measurement of Critical Strain for Solidification Cracking, Modelling of Casting, Welding, and Advanced Solidification Processes[END_REF][START_REF] Yamanaka | Critical Strain for Internal Crack Formation in Continuous Casting[END_REF] The lowest values were found at high strain rate and in crack-sensitive grades (e.g. high-sulfur peritectic steel). [START_REF] Yamanaka | Measurement of Critical Strain for Solidification Cracking, Modelling of Casting, Welding, and Advanced Solidification Processes[END_REF] In aluminum rich Al-Cu alloys, critical strains were reported from 0.09 to 1.6% and were relatively independent of strain rate. [START_REF] Wisniewski | Tensile Behavior of Solidifying Aluminum Alloys, Modelling of Casting, Welding, and Advanced Solidification Processes[END_REF] Tensile stress is also a requirement for hot tear formation. [START_REF] Yamanaka | Measurement of Critical Strain for Solidification Cracking, Modelling of Casting, Welding, and Advanced Solidification Processes[END_REF] The maximum tensile stress occurs just before formation of a critical flaw. [START_REF] Wisniewski | Tensile Behavior of Solidifying Aluminum Alloys, Modelling of Casting, Welding, and Advanced Solidification Processes[END_REF] The critical strain decreases with increasing strain rate, presumably because less time is available for liquid feeding, and also decreases for alloys with wider freezing ranges. Won et al [START_REF] Won | A New Criterion for Internal Crack Formation in Continuously Cast Steels[END_REF] suggested the following empirical equation for the critical strain in steel, based on fitting measurements from many bend tests: 0.3131 0.8638 0.02821

cr B T ε ε = ∆ ( 19
)
where ε is the strain rate and ∆T B is the brittle temperature range, defined between the temperatures corresponding to solid fractions of 0.9 and 0.99.

More mechanistically-based hot-tearing criteria include more of the local physical phenomena that give rise to hot tears. Feurer, [START_REF] Feurer | Mathematisches modell der Warmrissneigung von binären aluminium legierungen[END_REF] and more recently Rappaz et al. [START_REF] Rappaz | A New Hot-Tearing Criterion[END_REF] have proposed that hot tears form when the local interdendritic liquid feeding rate is not sufficient to balance the rate of tensile strain increase across the mushy zone. The criterion of Rappaz et al. predicts fracture when the strain rate exceeds a limit value that allows pore cavitation to separate the residual liquid film between the dendrites:

( )

        - - - ∇ ≥ H v p p T R S L S T C m S L l ρ ρ ρ ρ ρ µ λ ε 180 1 2 2 (27.2 0)
in which 2 λ is the secondary dendrite arm spacing, p m is the local pressure in the liquid ahead of the mushy zone, p C is the cavitation pressure, v T is the velocity of the solidification front. The quantities R and H depend on the solidification path of the alloy:

∫ ∆ =  ⌡ ⌠ =  ⌡ ⌠ ∇ = T T s T T l s T T l s dT g T T F dT g g H dT g T F g T R 2 1 2 1 2 1 ) ( ) ( 1 2 2 3 2 (27.2 1)
where the integration limits are calibration parameters which also have physical meaning. [START_REF] Drezet | Prediction of Hot Tears in DC-Cast Aluminum Billets[END_REF] The upper limit T 1 may be the liquidus or the coherency temperature, while the lower limit T 2 typically is within the solid fraction range of 0.95-0.99. [START_REF] Hamdi | The importance of viscoplastic strain rate in the formation of center cracks during the start-up phase of direct-chill cast aluminium extrusion ingots[END_REF]

Model Validation

Model validation with both analytical solutions and experiments is a crucial step in any computational analysis and thermo-mechanical modeling is no exception. Weiner and Boley 6 derived an analytical solution for unidirectional solidification of an unconstrained plate with a unique solidification temperature, an elastic-perfectly-plastic constitutive law and constant properties. The plate is subjected to sudden surface quench from a uniform initial temperature to a constant mold temperature. This benchmark problem is ideal for the validation of computational thermal-stress models, as it can be solved with a one-dimensional mesh, as shown in Figure 27.2. Numerical predictions should match with arbitrary precision according to the mesh refinement. For example, the solidification stress analysis code, CON2D [START_REF] Li | Thermo-Mechanical Finite-Element Model of Shell Behavior in Continuous Casting of Steel[END_REF] and the commercial code ABAQUS were applied for the conditions in Table 27.1. [START_REF] Koric | Efficient Thermo-Mechanical Model for Solidification Processes[END_REF] The solidification temperature interval is spread over a small (0.1 o C) artificial temperature range about the fusion temperature of 1494.4°C. The instantaneous surface quench is modeled with high convection, ( ) MPa-s. [START_REF] Koric | Efficient Thermo-Mechanical Model for Solidification Processes[END_REF] This represents a limit case for this elastic-visco-plastic expression. To model the unconstrained plate with a single row of elements, a generalized plane strain condition was imposed in the y and z directions (parallel to the surface) by coupling the displacements of all nodes along the edges of the slice domain as shown in Figure 27.2. [START_REF] Koric | Efficient Thermo-Mechanical Model for Solidification Processes[END_REF] The constant axial strain assumed under this condition is computed by satisfying a single extra scalar equation:

gap q h T T ∞ = - ,
z zz F dS = ∫ σ (27.

Application to Mold Casting

Braking disks

The finite element software THERCAST for thermomechanical analysis of solidification 66 has been used in automotive industry to predict distortions of braking discs made of grey iron and cast in sand molds. [START_REF] Bellet | Finite elements for a thermomechanical analysis of solidification processes[END_REF] A particular attention has been paid to the interaction between the deformation of internal sand cores and the cast parts. This demands a global coupled thermomechanical simulation, as presented above. This deformation causes a difference in thickness between the two braking tracks of the disc. Such a defect needs heavy and costly machining operations to get qualified parts. Instead, process simulation allows the manufacturer to test alternative geometries and process conditions in order to minimize the defect. <INSERT FIGURES 27.5, 27.6 and 27.7 HERE> Similar thermomechanical calculations have been made in the case of plain discs, leading to comparisons with residual stress measurements by means of neutrons and X-ray diffraction. [START_REF] David | Numerical simulation of casting processes. Benefits of thermomechanical simulation in automotive industry[END_REF] As shown in Figure 27.8, calculations are consistent with measurements, the difference being less than 10 MPa.

<INSERT FIGURE 27.8 HERE>

Application to Continuous Casting of Steel, Primary Cooling in Mold

The two-dimensional transient finite-element thermal-mechanical model, CON2D [START_REF] Li | Thermo-Mechanical Finite-Element Model of Shell Behavior in Continuous Casting of Steel[END_REF][START_REF] Koric | Efficient Thermo-Mechanical Model for Solidification Processes[END_REF] has been applied to predict temperature, displacement, strain, stress and hot tearing during primary cooling in continuous casting of steel. This Lagrangian model tracks a transverse slice through the strand as it moves downwards at the casting speed to reveal the entire 3-D stress state. The 2-D assumption produces reasonable temperature predictions because axial (z-direction) conduction is negligible relative to axial advection. [START_REF] Meng | Heat Transfer and Solidification Model of Continuous Slab Casting: CON1D[END_REF] In-plane mechanical predictions should be reasonable because bulging effects are small and the undiscretized casting direction is modeled with the appropriate condition of generalized plain strain.

The mechanical properties of steel used with this model were discussed in Sections 27.2.5 and 27.3.2 and in Table 27.1 and feature temperature-dependent elastic modulus and composition-dependent elastic-viscoplastic constitutive behavior and density. Two specific example applications with this model are presented next: 1) predicting the ideal taper of the mold to minimize gap formation during slab casting and 2) finding the critical casting speeds to avoid quality problems related to bulging below the mold.

Ideal Taper of a Slab Casting Mold

Thermal and mechanical behavior were simulated in a slice through the centerline of the wide face of the solidifying steel shell (Figure 27.2), as it moves down through the continuous casting mold. The total shrinkage strain predicted for this slice has been shown to provide an accurate and economical estimate of the ideal taper that should be applied to the narrow faces of the mold, in order to avoid formation of an air gap. The corner effects are reasonably small and were ignored in this analysis. The total heat flux (integrated from the heat flux profile) was forced to match an empirical equation that was obtained from a curve fit of many measurements under different conditions at a typical slab caster. [START_REF] Cicutti | Mould Thermal Evaluation in a Slab Continuous Casting Machine[END_REF] Computing ideal taper is useful application of computational models, in order to minimize the problems associated with excessive air gap formation, or pressure on defect formation and productivity.

Sample results are given in Figure 27.9 for a typical 200 mm thick slab mold, based on standard conditions of 1000 mm width, 800 mm working mold length, and 1.5 m/min casting speed. These figures show the effect of steel grade on heat flux, shell thickness, surface temperature, and ideal narrow-face taper, all as a function of distance down the mold wall. Taper is presented as %/mold. A mold flux with 1215 o C solidification temperature was assumed for the 0.13%C peritectic steel and a 1120 o C mold flux for the low and high carbon steels, which is typical of industrial practice. The practice of adopting mold powders with high solidification temperature and low viscosity was proposed by Wolf to produce lower, but more uniform heat transfer rates to help avoid cracks in depression-sensitive grades, such as peritectic steels. [START_REF] Wolf | Continuous Casting: Initial Solidification and Strand Surface Quality of Peritectic Steels[END_REF] Slags with opposite properties are used for low and high carbon steels, to help avoid sticker problems.

<INSERT FIGURE 27.9 HERE>

The results show that the higher solidification-temperature mold flux, combined with the rougher shell surface produce a lower heat flux for the peritectic steel (Figure 27.9a). The lower heat flux produces a thinner shell (Figure 27.9b). It also produces a hotter shell surface temperature (Figure 27.9c). This effect appears to outweigh the importance of the extra shrinkage of the peritectic steels. Thus, peritectic steels experience less shrinkage and require less taper than either low or high carbon steels (Figure 27.9d).

The low carbon steel (0.07%C) experiences more inelastic strain than other grades, owing to its microstructure being mainly composed of soft, delta phase. The delta-ferrite phase creeps faster than the stronger austenite phase found in other steels. The extra creep generated in the compressive solid surface layer tends to allow the interior shrinkage to have more effect in this grade. The net effect is that low carbon steel experiences a thicker shell with more shell shrinkage and taper than for the other grades.

As for every case studied in this work, significantly more mold taper is needed just below the meniscus than near mold exit. Thermal shrinkage strain dominates the need for taper.

However, before implementing new taper designs into an operating casting mold, these results must be modified to account for several phenomena which affect ideal taper. These include the thermal distortion of the narrow face, relative to its distortion at the meniscus, the expansion and thermal distortion of the wide face, the decrease in wide face perimeter due to the change in mold cavity dimensions down a funnel mold, and finally, the variation in thickness of the resolidified mold flux layers down the mold. More details are provided by Thomas and Ojeda, [START_REF] Thomas | Ideal Taper Prediction for Slab Casting, ISSTech Steelmaking Conference[END_REF] but clearly, more work is needed before the model results can be safely used in practice.

Maximum Casting Speed to Avoid Hot Tears in Billet Casting

The model was next applied to predict the maximum casting speed allowable, while avoiding excessive bulging and hot tearing. Simulations start at the meniscus, 100 mm below the top of the mold, and extend through the 800-mm long mold and below, for a caster with no sub-mould support.

The model domain is an L-shaped region of a 2-D transverse section, shown in Figure 27.10. Removing the center portion of the section, which is always liquid, saves computational cost and allows a pressure boundary condition that avoids stability problems related to element "locking" that may occur with fixed-displacement boundary conditions.

<INSERT FIGURE 27.10 HERE>

The instantaneous heat flux, given in Eq. (27.23), was based on fitting many plant measurements of total mold heat flux and differentiating. [START_REF] Li | Maximum Casting Speed for Continuous Cast Steel Billets Based on Sub-Mold Bulging Computation[END_REF] It was assumed to be uniform around the perimeter of the billet surface in order to simulate ideal taper and perfect contact between the shell and mold. Below the mold, the billet surface temperature was kept constant at its circumferential profile at mold exit. This eliminates the effect of spray cooling practice imperfections on sub-mold reheating or cooling and the associated complication for the stress/strain development. A typical plain carbon steel was studied (0.27%c, 1.52%Mn, 0.34%Si) with 1500.7°C liquidus temperature, and 1411.8 °C solidus temperature. <INSERT FIGURE 27.12 HERE> Figure 27.12(b) shows contours of "hoop" stress constructed by taking components in the x direction across the dendrites in the horizontal portion of the domain and the y direction in the vertical portion. High values appear at the off-corner sub-surface region, due to a hinging effect that the ferrostatic pressure over the entire face exerts around the corner. This bends the shell around the corner and generates high subsurface tensile stress at the weak solidification front in the off-corner subsurface location. This tensile stress peak increases slightly and moves towards the surface at higher casting speed. Stress concentration is less and the surface hoop stress is compressive at the lower casting speed. This indicates no possibility of surface cracking. However, tensile surface hoop stress is generated below the mold at high speed in 27.12(b) at the face center due to excessive bulging. This tensile stress, and the accompanying hot-tear strain, might contribute to longitudinal cracks which penetrate the surface.

Hot tearing was predicted using the criterion in Eq. 18 with the critical strain given in Eq. 19, and a 90% temperature of 1459.9 °C. Inelastic strain was accumulated for the component oriented normal to the dendrite growth direction, because that is the weakest direction and corresponds to the measurements used to obtain Eq. 27.19. .Figure 27.12(c) shows contours of hot-tear strain in the hoop direction. The highest values appear at the off-corner sub-surface region in the hoop direction. Moreover, significantly higher values are found at higher casting speeds. For this particular example, hot-tear strain exceeds the threshold at 12 nodes, all located near the off-corner subsurface region. This is caused by the hinging mechanism around the corner. No nodes fail at the center surface, in spite of the high tensile stress there. The predicted hot-tearing region matches the location of off-corner longitudinal cracks observed in sections through real solidifying shells, such as the one pictured in Figure 27.13. The bulged shape is also similar.

<INSERT FIGURE 27.13 HERE>

Results from many computations were used to find the critical speed to avoid hot tear cracks as a function of section size and working mold length, presented in Figure 27.14. [START_REF] Li | Thermo-mechanical Finite Element Model of Shell Behavior in Continuous Casting of Steel, Modeling of Casting, Welding and Advanced Solidification Process X[END_REF] These predictions slightly exceed plant practice, which is generally chosen by empirical trial and error. This suggests that plant conditions such as mold taper are less than ideal, that other factors limit casting speed, or those speeds in practice could be increased. The qualitative trends are the same. This quantitative model of hot tearing has enabled many useful insights into the continuous casting process. Larger section sizes are more susceptible to bending around the corner, so have a lower critical speed, resulting in less productivity increase than expected. The trend towards longer molds over the past three decades enables a higher casting speed without cracks by producing a thicker, stronger shell at mold exit.

<INSERT FIGURE 27.14 HERE>

4. 3 Application to Continuous Casting of Steel, Secondary Cooling

Thermomechanical simulations are used by steelmakers to analyze stresses and strains all along the secondary cooling zone. One of their thrust is the prediction the bulging of the solidified crust between the supporting rolls which is responsible for the tensile stress state in the mushy core, which in turn induces central macrosegregations [START_REF] Miyazawa | Macrosegregation in continuously cast steel slabs: preliminary theoretical investigation on the effect of steady state bulging[END_REF][START_REF] Lesoult | Analysis and prevention of centreline segregation during continuous casting of steel related to deformation of the solid phase[END_REF] Two and three-dimensional finite element models have been recently developed, based on an original "global non steady-state" approach which provides results on the whole length of the caster. The approach implemented in THERCAST software is described in detail elsewhere. [START_REF] Costes | 3D thermomechanical simulation of the secondary cooling zone of steel continuous casting[END_REF][START_REF] Bellet | A two-dimensional finite element thermomechanical approach to a global stress-strain analysis of steel continuous casting[END_REF] The constitutive models are those presented in section 27.2.5 The contact with supporting rolls is controlled by a penalty formulation derived from what has been presented in section 27.3.2 using penalty coefficients attached to the different rolls and continuously adapted in order to control the numerical penetration of the strand. Figure 27.15 shows the results obtained on a vertical-curved machine (strand thickness 0.22 m, casting speed 0.9 m/min, material Fe-0.06wt%C) in a region located around 11 m below the meniscus. The pressure distribution reveals a double alternation of compressive and depressive zones. First, along the strand surface, the material is in a compressive state under rolls where the pressure reaches its maximum, 36 MPa. Conversely, it is in a depressive (tensile) state between rolls, where the pressure is minimum (-9 MPa).

<INSERT FIGURE 27.15 HERE>

Examination of the pressure state within the solid shell close to the solidification front (i.e. close to the solidus isotherm), reveals that the stress alternates. The steel is in a tensile state (negative pressure of about -2 MPa) when passing in front of rolls, while it is in a compressive state in between, the value of pressure being around 2 to 3 MPa. These results agree with previous structural analyses of the deformation of the solidified shell between rolls, such as those carried out in static conditions by Wünnenberg, [START_REF] Wunnenberg | Strand bulging between supporting rollers during continuous slab casting[END_REF] Miyazawa and Schwerdtfeger [START_REF] Miyazawa | Macrosegregation in continuously cast steel slabs: preliminary theoretical investigation on the effect of steady state bulging[END_REF] or by Kajitani et al. [START_REF] Kajitani | Numerical simulation of deformation-induced segregation in continuous casting of steel[END_REF] on limited slab sections moving downstream between rolls and submitted to the metallurgical pressure onto the solidification front. The influence of process parameters on the thermomechanical state of the strand can then be studied using such numerical models. An example is given by Figure 27.16, presenting the sensitivity of bulging to the casting speed. It can also be seen that bulging predictions are sensitive to the roll pitch, a larger pitch between two sets of rolls inducing an increased bulging. These numerical simulations can then be used to study possible modifications in the design of continuous casters, such as the replacement of large rolls by smaller ones in order to reduce the pitch and the associated bulging. [START_REF] Triolet | Mastering steel slab internal soundness and surface quality issues through thermomechanical modelling of continuous casting[END_REF]

CONCLUSIONS

Mechanical analysis of casting processes is growing in sophistication, accuracy, and phenomena incorporated. Quantitative predictions of temperature, deformation, strain, stress, and hot tearing in real casting processes are becoming possible. Computations are still hampered by the computational speed and limits of mesh resolution, especially for realistic three-dimensional geometries and defect analysis. Further developments are needed in fundamental models of defect formation over multiple length scales, and in efficient convergence algorithms to solve the equations. Fundamental measurements, including the quantification of interfacial conditions, material constitutive properties and conditions for defect formation are also needed. However, solidification processes are growing in maturity and improvements gained by trial and error in the plant are becoming more expensive. As computing power and software tools continue to advance, advanced computational models will become even more important in the years to come. Future advances to casting processes will increasingly rely on advanced computational models such as the thermo-mechanical models discussed here.

.

  which K is the viscoplastic consistency and m the strain rate sensitivity. Denoting the von Mises equivalent strain rate scalar, Eq. (27.3) yields the well known power law: It can be noticed that the preceding Newtonian model is actually a particular case of the non-Newtonian one: Eq. (27.1) can be derived from Eq. (27.3shrinkage is included in the third equation, as we can write in the solidification interval ρ the densities at the solidus and liquidus temperatures, respectively. Hence, we have:

( 27 . 5 )

 275 to model the different states. Using a physically-reasonable (high) value for the elastic modulus.(27.5) to approximate Eq. (27.3) for a Newtonian fluid, as the generated elastic strains are very small. Another approach[START_REF] Cruchaga | Modeling fluid-solid thermomechanical interactions in casting processes[END_REF][START_REF] Bellet | An ALE-FEM approach to the thermomechanics of solidification processes with application to the prediction of pipe shrinkage[END_REF] uses a viscoplastic equation for the mushy and liquid states and an elastic-visco-plastic equation for the solid state.

  % )][1 0.008(% )] 100[8011 0.47 ( )] [100 (% )][1 0.013(% )] 7100 73(% ) [0.8 0.09(% )][ ( ) 1550]

  the thermal conductivity of the gas, g the gap thickness, c T and m T the local surface temperature of the casting and mold, respectively, c ε and m ε their gray-body emissivities, σ the Stefan- Boltzmann constant. It is to be noted that the conductive part of the flux can be written in more detail to take into account the presence of coating layers on the mold surface: conduction through a medium of thickness coat g , of conductivity ( ) coat k T . It can be seen that the first term tends to infinity as the gap thickness tends to zero: this expresses a perfect contact condition, c T and m T tending towards a unique interface temperature. The reality is somewhat different, showing always non perfect contact conditions.

Figure 27 .

 27 Figure 27.1 is a graphic representation of the heat exchange coefficient in both cases. It should

  -0.15wt%C steel grade. The local condition for fracture initiation is then:

  to lessen the instabilities caused by these extreme conditions. The elastic-perfectly-

2 2 )

 2 Figures 27.3 and 27.4 compare the temperature and stress profiles in the plate at two

Figure 27 .

 27 5 illustrates the discretization of the different domains involved in the calculation. The actual cooling scenario has been simulated: cooling in mold during 45 min, shake out and air cooling during 15 min. Figure 27.6 gives the temperature evolution in different points located in a horizontal cross section at midheight in the disc, showing the influence of different physical phenomena : solidification after 2 min, solid state phase change after 20 min. The deformation of the core has been calculated, bringing out the thermal buckling of core blades. They are exposed to very high temperature, and their dilatation is too much constrained, resulting in their deformation, as shown in Figure 27.7.

  Figure 27.12(a) shows the distorted temperature contours near the strand corner at 200 mm

Figure 27 . 1 .

 271 Figure 27.1. Modeling of the local heat transfer coefficient in the gap and effective contact situations.

Figure 27 . 2 .

 272 Figure 27.2. One-dimensional slice-domain for modeling solidifying plate.

Figure 27 . 3 .

 273 Figure 27.3. Temperatures through solidifying plate at different times comparing analytical solution and numerical predictions.

Figure 27 . 4 .

 274 Figure 27.4. Transverse (Y and Z) stress through solidifying plate at different times comparing analytical solution and numerical predictions.

Figure 27 . 5 .

 275 Figure 27.5. Finite element meshes of the different domains: part, core, and two half molds.

Figure 27 . 6 .

 276 Figure 27.6. Temperature evolution in the part at different points located in the indicated section.

Figure 27 . 9

 279 Figure 27.9 Effect of steel grade on heat flux (a), shell thickness (b), surface temperature (c), and

Figure 27 . 10 .Figure 27 . 10 .

 27102710 Figure 27.10. Model domain.

| | | | | | | | | | | | | | | | | | | | | | || | | | X X X X X X X X X X X X X X X X X X X X X X X X XDISTANCE BELOW M ENISCUS (mm) SHELL THICKNESS (mm)

| | | | | | | | | | | | | | | | | | | | | | | | X X X X X X X X X X X X X X X X X X X X X X X

Acknowledgements

The authors wish to thank the Continuous Casting Consortium and the National Center for Supercomputing Applications at the University of Illinois, the French Ministry of Industry, the French Technical Center of Casting Industries (CTIF) and the companies Arcelor, Ascometal,

Atlantic Industrie, Aubert et Duval, Erasteel, Industeel and PSA Peugeot-Citroën, for support of this work.

Temperature [°C]

Vertical section Figure 27.7. Deformation of core blades in a radial section, after a few seconds of cooling. On the left, displacements have been magnified by a factor 100. The temperature distribution is superimposed. On the right, the difference in thickness between the two braking tracks is shown.

DISTANCE-BELOW -M ENISCU S(mm) SHELL SURFACE TEMPERATURE (C)


Bastien Fernandez

Lev S Tsimring

Co-Repressive Interaction and Clustering of Degrade-and-Fire Oscillators

Keywords: numbers: 05.45Xt, 87.16Yc, 87.18Cf

Strongly nonlinear degrade-and-fire (DF) oscillations may emerge in genetic circuits having a delayed negative feedback loop as their core element. Here we study the synchronization of DF oscillators coupled through a common repressor field. For weak coupling, initially distinct oscillators remain de-synchronized. For stronger coupling, oscillators can be forced to wait in the repressed state until the global repressor field is sufficiently degraded, and then they fire simultaneously forming a synchronized cluster. Our analytical theory provides necessary and sufficient conditions for clustering and specifies the maximum the number of clusters which can be formed in the asymptotic regime. We find that in the thermodynamic limit a phase transition occurs at a certain coupling strength from the weakly-clustered regime with only microscopic clusters to a strongly clustered regime when at least one giant cluster has to be present.

Many gene regulatory networks contain negative feedback loops as their core elements [START_REF] Alon | An introduction to systems biology: design principles of biological circuits[END_REF]. The negative feedback provides robustness and noise-resistance to signaling pathways [START_REF] Stelling | [END_REF]. They also lie at the center of many natural rhythmic circuits (such as circadian clocks [3]) and synthetic gene oscillators [4][5][6]. In our recent work [7] we identified the delayed negative feedback as the mechanism of oscillations in the synthetic two-gene oscillator [5]. It is well known that delayed auto-repression can lead to oscillatory gene expression even with only a single regulatory element [8][9][10][START_REF] Bratsun | Proc. Natl. Acad. Sci. USA[END_REF][START_REF] Coutinho | [END_REF]. Typically, the period of delay-induced oscillations is comparable with the delay time and thus this mechanism seemingly could not explain rather slow oscillations observed in gene circuits. However, we showed that in a strongly nonlinear regime, the period of oscillations is determined by the rate of enzymatic degradation of the repressor protein and therefore can be arbitrarily longer than the transcriptional and translational delay. The essential mechanism of oscillations in this circuit is based on the periodic alternation of two regimes -slow degradation of repressor protein, and the following fast production of the repressor (firing), hence we termed these sawtooth-like oscillations "degrade-and-fire" in analogy with integrate-and-fire oscillations in neuronal circuits. Bacterial cells containing the DF circuit demonstrated robust oscillations, but due to cell-to-cell variability and stochastic effects, in the absence of cell-cell communication the bacterial colony quickly becomes desynchronized. In our subsequent work [13] we succeeded in synchronizing gene oscillators across bacterial colony, however the synchronized regime was achieved with a different circuit design that incorporated a quorum-sensing machinery producing small molecule AHL which served as signaling element necessary for cellcell synchronization.

In this paper we study synchronization of DF oscillators through purely co-repressive interaction. To enable the analytical calculations, we replace the original non-linear delay-differential model of the DF oscillator [7] by a piecewise linear model which assumes that the concentration x degrades from the maximum value (which can be scaled to 1) linearly with unit rate, ẋ = -1 (mimicking enzymatic decay with high enzyme affinity), until it it reaches the (small) threshold value η > 0, after which repressor concentration x is immediately reset to 1 ("fire"), and the process repeats. It is easy to see that the protein concentration oscillates with period 1 -η. Now we consider a population of N degrade-and-fire (DF) oscillators coupled through a common repressor field. We postulate that the dynamics of the concentration of repressor protein x i is given by the following rules • if positive, x i degrades linearly with rate 1, i.e. ẋi = -1, or remains constant if it has reached 0,

• when the locally averaged concentration

χ i (t) = (1 -ǫ)x i (t) + ǫX(t) (where X(t) = 1 N N i=1
x i (t) and 0 ǫ 1 is the coupling strength parameter) reaches threshold η, the ith oscillator fires, and its concentration is reset to 1, i.e. x i (t+) = 1.

With these simple evolution rules, one can readily see that all oscillators must fire indefinitely (i.e. there cannot be "oscillator death" regime for any of the oscillators). Indeed, starting from an arbitrary configuration {x i } N i=1 for which all χ i > η, all concentrations decay towards 0 with time (those that have reached x i = 0, remain at zero) and so do all χ i . Thus the oscillator with the lowest x i (possibly, more than one if several oscillators have identical concentrations x i ) eventually fires when the corresponding χ i reaches η. After that, the oscillator j with new lowest x j has to fire when its χ j reaches η, and so on [14]. It is also clear that if any two oscillators in a population are in sync at certain time t * , i.e. x i = x j , they will remain in sync for all t > t * . What is not obvious however, is under which conditions oscillators that are initially out of sync will synchronize in the course of the dynamics, and what the properties of the resulting clusters are. To answer these questions, we begin with the simple case of two oscillators.

Interaction of two DF oscillators. In this case oscillator i (i = 1, 2) fires when χ i = (1 -ǫ/2)x i + ǫx 3-i /2 = η. Without loss of generality, we can assume that initially, one oscillator has just fired (i.e. (x 1 , x 2 ) = (x, 1) with η < x 1) and denote t f (x) > 0 the firing time associated with this configuration. The dynamics can be described through the return map R between ordered concentration pairs just after firings, (x, 1) → R(x, 1) := (1 -t f (x), 1).

The 2-dimensional dynamics of two DF oscillators is contained in a parameter dependent subset of the unit square [0, 1] 2 (see Fig. 1). For ǫ = 0, firing occurs when either of the individual concentrations x i = χ i touches the threshold η. A simple computation shows that R 2 (x, 1) = (x, 1) for all x, i.e. we have a continuum of periodic orbits equivalent to the rigid rotation on a 2-dimensional torus (see Fig. 1a).

(t (1- )) = 1 η = χ 1 x 2 x 1 1 0 η 0 η η χ = 2 a x 1 x (0) = 2 x x f +η 1 f x (t ( ))
x For ǫ > 0, the dynamics changes drastically, as expected. Instead of a continuum of neutral periodic orbits, a single stable periodic orbit emerges which attracts all trajectories, except for the unstable periodic orbit lying on the diagonal x 1 = x 2 (see Fig. 1b). These are direct consequences of Lemmas 2 and 3 below. Furthermore, it can be shown that a unique and globally stable periodic orbit exists in arbitrary systems of two coupled DF oscillators with any monotonous degradation of both concentrations. Note that this result implies a somewhat unexpected corollary that two initially distinct oscillators never synchronize, and always remain distinct. Figure 2 shows typical trajectories of two coupled oscillators for η = 0.01 and different values of ǫ. As seen from the Figure, the two oscillators in the asymptotic regime are always in anti-phase, and, for large ǫ, one oscillator remains repressed (x i = 0) almost the entire time when the other repressor is present (more precisely, it fires when the other concentration reaches the small value 2η/ǫ), and vice versa.

2 η = χ 2 x 1 η = χ η η 1 0 1 0 b x * 1 ( ,1)
Many coupled DF oscillators. Now we turn to the analysis of the clustering dynamics of a population of N corepressively coupled DF oscillators. Before we proceed with analytical results, we illustrate the typical behavior of the system numerically. Figure 3 shows the time series of 50 oscillators with initially distinct values of x for η = 0.01 and different values of the coupling parameter ǫ. For small ǫ = 0.01 (Fig. 3a), all oscillators fire before they ever reach zero, and so they remain distinct (no clustering). For larger ǫ = 0.03 (Fig. 3b) some oscillators reach zero, and some of them (not all) synchronize, and thus asymptotically, only 8 clusters remain. Note that this number is not universal, as it depends on specific initial conditions. For ǫ = 0.05 (Fig. 3c), only two clusters remain. Their "weights" however are not equal, in this example one cluster contains 24 oscillators, and another 26 oscillators. Again, this weight distribution is nonuniversal. And for very large ǫ = 0.5 (Fig. 3d) oscillators produce repressor protein in almost exact anti-phase, one one repressor is present, the other one is absent. By grouping oscillators with identical value of x i into one cluster, the population dynamics can be described via {(n k , x k )} K k=1 where n k ∈ {1, ..., N } denotes the size of the cluster k and x k the corresponding repressor concentration (K N is the total number of clusters). In this viewpoint, the cluster size distribution {n k } obviously remains unaffected in time unless two clusters k and k ′ fire together.

As before, we consider the return map between consecutive firings. Any ordering in {(n k , x k )} is irrelevant because of the permutation symmetry in this system. However, it is more convenient to deal with ordered values of x k . Thus we assume that 0

< x 1 < x 2 < • • • < x K-1 < x K = 1.
Given {n k }, this defines the state configuration set where the return map R ≡ R {n k } is effectively defined, viz.

T {n k } = {(x 1 , ..., x K ) : 0 < x 1 < • • • < x K = 1 and (1 -ǫ)x 1 + ǫ N K k=1 n k x k > η .
In order to maintain the ordering of protein concentrations in time, we must include cyclic permutations of indices in the return map dynamics. The absence of clustering requires that, when starting in T {n k } , the after-firing configuration has to lie in T {n k+1 } (where all indices are understood mod K) and this should hold for each k = 1, ..., K per where K per is the minimal size distribution period. This property globally holds in state configuration sets provided that the coupling strength is small enough.

Lemma 1 R {n k } T {n k } ⊂ T {n k+1 } for all k = 1, ...K per iff ǫ η min k=1,...,Kper N N -n k -n k+1
where ǫ η := ǫ/η (1)

In particular, for ǫ η < N/(N -2) which is the minimal RHS in (1), no clustering can ever occur, independently of the initial configuration. Notice that N/(N -2) > 1 and → 1 as N → ∞. Lemma 1 is rigorously proven in [15], but it can be intuitively understood as follows.

The only way two clusters can merge is when one cluster reaches zero (collapse) and remains at zero at least until the following cluster also reaches zero. It is easy to see that in the worst case scenario two clusters 1 and 2 should be very close to zero when the other K -2 clusters are very close to 1. In the limit x 1,2 → 0, x 3 , ..., x K → 1, after x 1 has collapsed, χ 1 reaches η when the value of x 2 is equal to [N -ǫ η (N -n 1 -n 2 )]/(N ǫ η ), and in order to avoid merging of n 1 and n 2 , this value has to be positive. Of course, the same condition has to be satisfied for all consecutive pairs n k , n k+1 , hence the minimum in inequality [START_REF] Alon | An introduction to systems biology: design principles of biological circuits[END_REF]. Independently of (1), we may ask about the fate of the trajectories for which clustering occurs. Since the sequence of merging events is always finite, this actually would determine every possible asymptotic regime. In absence of clustering, the return map in T {n k } becomes R * {n k } (x 1 , ..., x K ) := (x 2 -t f , ..., x K -t f , 1) where t f still denotes the firing time. (This map is actually

K -1- dimensional because (R * {n k } x) K = x K = 1.) Since the image R * {n k }
x belongs to T {n k+1 } , one needs to iterate further until the point returns back to T {n k } ; hence the dynamics to study is the composed map

R * {n k+Kper-1 } • • • • • R * {n k+1 } • R * {n k } .
It can be shown [15] that this map is always a pure contraction.

Lemma 2 For every {n

k } K k=1 and 0 < ǫ 1, there is a norm in R K-1 for which R * {n k+Kper -1 } • • • • • R * {n k+1 } • R * {n k } is a global contraction.
It follows that, in absence of clustering, every trajectory initially in T {n k } must approach a unique periodic orbit (whose single element in T {n k } is the fixed point of the composed map). More generally, every trajectory in a population of N oscillators must converge to the periodic orbit associated with some {n k }. Hence, the asymptotic cluster size distributions and periodic orbits are in one-toone correspondence. Accordingly, to get asymptotically attainable cluster distributions, it suffices to compute the existence domains of the corresponding periodic orbits. These domains are given by the following statement.

Lemma 3 The periodic orbit in T {n k } exists iff

ǫ η < ζ({n k }) := 2N 2 N 2 -k n 2 k 1 1 -min k n k N -n k+1 (2)
This claim is proved in [15] where we also give the explicit expressions for the periodic orbit. Note that when n K = N -n 1 (i.e. K = 2), we have ζ({n k }) = +∞, which means that every two-cluster periodic orbit exists for any ǫ. In particular, for n 1 = n 2 = 1 this observation implies the result for two coupled oscillators presented above.

For K > 2, the critical value ζ({n k }) is minimal for the (equi-)distribution where all n k = 1 (in this case K = N ) and the corresponding minimum is ζ c = 2N/(N -2). So all possible periodic orbits exist for ǫ η up to ζ c . For larger ǫ η , some periodic orbits disappear -certainly the one associated with equi-distribution -and the number of clusters K must eventually be less than N .

A closer look at the formula (2) reveals that when {n k } is only composed of microscopic clusters, i.e. when k n 2 k = O(N ), the corresponding critical value ζ({n k }) approaches 2 in the thermodynamic limit N → ∞. Therefore, such distributions do not perdure asymptotically beyond ǫ η ≃ 2 when N is large, and at least some of the clusters size n k reach O(N ) in the course of time. Moreover, distributions that contain macroscopic cluster(s) (i.e. n k ≃ ρ k N with ρ k > 0 for some k) all have critical value ζ({n k }) larger than 2 in the thermodynamic limit. Thus, there is a sharp transition at a certain ǫ c = ζ c η → 2η in the system for large N , from a dynamical regime in which all kinds of cluster size distributions may exist asymptotically, to a regime where every asymptotic distribution contains at least one giant O(N ) cluster.

Maximal number of clusters. The transition at ǫ η = ζ c can be quantitatively described by estimating the number of asymptotic clusters as a function of ǫ η . Since the precise number actually depends on the initial condition and can be as little as 2, the appropriate quantity to compute is the maximal possible number K max . It can be proven [15] that for any given K 3 and N K (and N = 6 if K = 3), the critical value ζ({n k } K k=1 ) in (2) reaches its maximum for the distribution consisting of K -1 unitary clusters of weight n k = 1 and one big cluster n K = N -K + 1. In this case, we have and therefore all distributions with K clusters exist for 0 ǫ η < ζ max . This immediately yields the maximum number clusters N and :

ζ({n k }) = ζ max (K) = 2N 2 (N -1) (K -1)(2N -K)(N -2)
(ǫ η ) = N if ǫ η < ζ c = N - 1 2 + 1 4 + N (N -1) 1 - ζ c ǫ η (3) if ζ c ǫ η < ζ max (3) = 2 if ǫ η ζ max (3)
This result indicates that a kind of second order phase transition takes place at ǫ η = ζ c with K max that behaves like N (1-1 -ζ c /ǫ η ) slightly above the threshold. Formula (3) implies that for each 1 K K max (ǫ η ), exists non-empty set which produces a stable periodic trajectory with K clusters. However, the "typical" number of clusters emerging from an arbitrary set of initial conditions can be much less than this upper bound. We performed numerical simulations of the full model with initial conditions randomly selected from the [0, 1] N hypercube. Figure 4 shows that the typical number of clusters for a given ǫ can be well approximated by the number K e corresponding to the uniformly distributed initial condition x i = i/N . However, the number of asymptotic clusters depends strongly on the width of the distribution of initial states. We performed simulations for initial conditions x i equi-distributed between 1 -δ and 1, and found that the number of asymptotic clusters strongly increases as the width of the initial distribution δ gets smaller, and approaches the values close to the upper bound K max (Fig. 4, inset).

In summary, we presented an analytically solvable model of co-repressive coupling of degrade-and-fire oscillators in a strongly nonlinear regime. At a certain coupling strength this model exhibits a phase transition from a regime when all cluster distributions including nonclustered state can be attained to the regime of strong clustering characterized by the appearance of at least one giant O(N ) cluster. Our model neglected the variability in parameters of individual oscillators and possible delays in the coupling. We plan to address these effects in our future work. The experimental realization of coupled gene oscillators [13] incorporated a different (co-excitatory) mechanism of coupling through a quorum-sensing mechanism. However, we believe that co-repressive coupling may also be realized experimentally, if for example oscillators produce an enzyme that degrade a freely diffusing inducer or the quorum sensing molecule activates a repressor protein. Our results demonstrate that a strong clustering can be expected in this case, and large out-of-phase clusters are likely to emerge.

xFIG. 1 :FIG. 2 :

 12 FIG. 1: Typical phase trajectories of two coupled DF oscillators. Simple arrows indicate motions between firings when both concentrations decay with rate 1. Triple arrows denote firings when the trajectory jumps from critical line χi = η to xi = 1. (a) ǫ = 0, every trajectory is periodic. (b) ǫ > 0, for every initial condition except x1 = x2 the system asymptotically converges to a unique periodic trajectory passing through (x * , 1) where x * = 1 -2(1 -η)/(4 -ǫ) (thick red line).

FIG. 3 :

 3 FIG. 3: Typical time series of 50 coupled DF oscillators and of their mean value X(t) (thick/blue line) for η = 0.01 and different values of ǫ: (a) 0.01; (b) 0.03; (c) 0.05; (d) 0.5.

FIG. 4 :

 4 FIG. 4: Number of clusters the 1000 coupled with η = 0.01 and 1000 different random initial conditions for each ǫ. The solid blue line indicates the upper bound (3), and the solid red line shows the number of clusters Ke(ǫ) for the uniform initial distribution. Inset: the number of clusters as a function of the initial distribution width δ for ǫ = 0.03

Acknowledgements. LT acknowledges useful discussions with Jeff Hasty and William Mather and is grateful to CPT for hospitality during his stay at the University of Provence, Marseille. BF thanks the Courant Institute (NYU) for hospitality. This work was supported by EU Marie Curie fellowship PIOF-GA-2009-235741 and NIH grants RO1-GM69811-01 and P50-GM085764.