Brian G Thomas
email: bgthomas@uiuc.edu

Michel Bellet
email: michel.bellet@ensmp.fr

Modeling of Stress, Distortion and Hot Tearing

Introduction

Computational modeling of mechanical behavior during solidification is growing in importance. This is because thermal and microstructural simulations alone are insufficient to predict the quality of the final product that is desired by the casting industry. Accurate calculation of displacements, strains, and stresses during the casting process is needed to predict residual stress and distortion, and defects such as the formation of cracks such as hot tears. It also helps in predicting porosity and segregation. As computing power and software tools for computational mechanics advance, it is becoming increasingly possible to perform useful mechanical analysis of castings and these important related behaviors.

The thermomechanical analysis of castings presents a formidable challenge for many reasons:

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

• Predicting distortions and residual stresses in cast products requires calculation of the history of the cast product and its environment over huge temperature intervals. This makes the mechanical problem highly non-linear, involving liquid-solid interaction and complex constitutive equations. Even identifying the numerous metallurgical parameters involved in those relations is a daunting 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.

• Accounting for the mold and its interaction with the casting makes the problem multidomain, usually involving numerous deformable components with coupled interactions.

• Cast parts usually have very complex three-dimensional shapes, which puts great demands on the interface between CAD design and the mechanical solvers, and on computational resources.

• The important length scales range from microns (dendrite arm shapes) to tens of meters (metallurgical length of a continuous caster), with similar huge order-of-magnitude range in time scales. 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, and a few examples of recent applications in shape castings and continuous casting.

Governing Equations

The modeling of mechanical behavior requires solution of 1) the equilibrium or momentum equations, relating force and stress; 2) the constitutive equations, relating stress and strain and 3) compatibility equations, relating strain and displacement. This is because the boundary conditions specify either force or displacement at different boundary regions of the domain Ω :

Equilibrium and compatibility equations

At any time and location in the solidifying material, the conservation of force (steady-state equilibrium) or momentum (transient conditions) can be expressed by:

d 0 d t ∇⋅ + ρ -ρ = v σ g (2)
where σ is the stress tensor, ρ is the density, g denotes gravity, v is the velocity field and t d d/ denotes the total (particular) time-derivation.. Stress can be further split into the deviatoric stress tensor and the pressure field. The different approaches for simplifying and solving these equations are discussed in Section 2.5.

Once the material has solidified, the internal and gravity forces dominate, so the inertia terms in Eq.

(2) can be neglected. Furthermore, the strains which dominate thermo-mechanical behavior during solidification are on the order of only a few percent, or cracks will form. With small gradients of spatial displacement,

x u u ∂ ∂ = ∇
/ , and the compatibility equations simplify to the following [START_REF] Mase | Continuum mechanics for engineers[END_REF] :

where ε ε ε ε is the strain tensor and u is the displacement vector. This small-strain assumption simplifies the analysis considerably. The compatibility equations can also be expressed as a rate formulation, where strains become strain rates, and displacements become velocities. This formulation is more convenient for a transient computation with time integration involving fluid flow and / or large deformation.

In casting analysis, the cast material may be in the liquid, mushy or solid state. Each of these states has different constitutive behavior, as discussed in the next 3 sections.

Liquid state constitutive models

( ) T ) ( 2 1 u u ε ∇ + ∇ = (3)
Metallic alloys generally behave as Newtonian fluids. Including thermal dilatation effects, the constitutive equation can be expressed as follows.

I s ε t l d d 3 1 2 1 ρ ρ µ - = & (4)
The strain rate tensor ε & is split into two components: a mechanical part, which varies linearly with the deviatoric stress tensor s, and a thermal part. In this equation, l µ is the dynamic viscosity of the liquid, ρ 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 (5)
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] 2.3 Mushy-state constitutive models Metallic alloys in the mushy state are two-phase liquid-solid media. Their mechanical response depends greatly on the local microstructural evolution, which involves several complex physical phenomena. An accurate description of these phenomena is useful for studying hot tearing or macrosegregation. Knowledge of the liquid flow in the mushy zone is necessary to calculate the transport of chemical species (alloying elements) [START_REF] Beckermann | Macrosegregation[END_REF] . Knowledge of the deformation of the solid phase is important when it affects liquid flow in the mushy zone by "sponge-effects". [START_REF] Lesoult | Segregation during solidification with spongy deformation of the mushy zone[END_REF] In such cases, two-phase models must be used. Starting from microscopic models describing the intrinsic behavior of the liquid phase and the solid phase, spatial averaging procedures must be developed to express the behavior of the compressible solid continuum and of the liquid phase that flows through it [START_REF] Nicolli | Modeling of macrosegregation caused by volumetric deformation in a coherent mushy zone[END_REF][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][START_REF] Ludwig | Rheological behaviour of Al-Cu alloys during solidification: constitutive modelling, experimental identification and numerical study[END_REF] .

If a detailed description is not really needed, such as in the analysis of residual stresses and distortions, the mushy state can be approximated as a single continuum that behaves as a non-Newtonian (i.e. viscoplastic) fluid, according to the following equations ( 6) to (8). Thus, the liquid phase is not distinguished from the solid phase, and the individual dendrites and grain boundaries are not resolved.

th vp ε ε ε & & & + = (6) s ε m eq vp K - = 1 ) ( 2 3 ε& & (7) I ε t th d d 3 1 ρ ρ - = & (8)
K is the viscoplastic consistency and m the strain rate sensitivity. Denoting ij ij eq s s 2 3 = σ the von Mises equivalent stress scalar, and

vp ij vp ij eq ε ε ε & & & 3 2 =
the von Mises equivalent strain rate scalar, Eq. ( 7) yields the well known power law:

m eq eq K ) (ε σ & =
. Note that the preceding Newtonian liquid model is actually a particular case of this non-Newtonian one: Eq. ( 4) can be derived from Eq. ( 6), ( 7) and ( 8) taking 1 = m and

l K µ 3 =
. The solidification shrinkage is included in Eq. ( 8): writing

L l S s g g ρ ρ ρ + =
in the solidification interval, ( S ρ , L ρ densities at the solidus and liquidus temperatures respectively, g s , g l volume fractions of solid and liquid respectively), the thermal strain rate is defined as follows:

t g t g t s L S L s L S th d d d d ) ( 1 d d 1 tr ρ ρ ρ ρ ρ ρ ρ ρ - ≈ - - = - = ε & (9)

Solid-state constitutive models

In the solid state, metallic alloys can be modeled either as elastic-plastic or elastic-viscoplastic materials. In the latter 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][START_REF] Agelet De Saracibar | On the constitutive modeling of coupled thermomechanical phase-change problems[END_REF] th in el

ε ε ε ε + + = (10)
I σ σ I σ σ ε ) tr( 1 ) tr( 1       ∂ ∂ -       + ∂ ∂ + - + = E T T E T T E E el ν ν ν ν (11) s ε m eq eq in K / 1 0 2 3 σ σ σ - = & (12)
I ε t th d d 3 1 ρ ρ - = & (13)
The strain rate tensor ε & is split into an elastic component, an inelastic (non reversible) component, and a thermal component. Equation (11) is the hypoelastic Hooke's law, where E is Young's modulus, ν the Poisson's coefficient, and σ & a time derivative of the stress tensor σ . Equation (12) 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 both 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 σ σ ε = + & ( 14
)
Inserting this into Eq. ( 12) simplifies it to:

s ε s ε eq eq in eq eq in σ ε σ ε 2 d 3 d form, l incrementa in or, , 2 3 = = & & (15)
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. ( 15) still holds, but the flow stress is independent of the strain rate. It may depend on the accumulated plastic strain because of strain hardening.

Implementation Issues

One of the major difficulties in the thermomechanical analysis of casting processes is the concurrent presence of liquid, mushy and solid regions which move with time as solidification progresses. Several different strategies have been developed, according to the process and model objectives.

• A first strategy consists in extracting the solidified regions of the casting domain based on the thermal analysis results. Then, a small-strain thermo-mechanical analysis is carried out on just this solid subdomain, using a standard solid-state constitutive model. Beside difficulties with the extraction process, especially when the solidified regions have complex unconnected shapes, this method may have numerical problems with the application of the liquid hydrostatic pressure onto the new internal boundary of the solidified region. However, this simple strategy is very convenient for many practical problems, especially when the solidification front is stationary, such as the primary cooling of continuous casting of aluminum, [START_REF] Sengupta | Quantification of temperature, stress, and strain fields during the start=up phase of direct chill casting process by using a 3D fully coupled thermal and stress model for AA5182 ingots[END_REF] and steel [START_REF] Kelly | Initial Development of Thermal and Stress Fields in Continuously Cast Steel Billets[END_REF][START_REF] Huespe | Thermomechanical model of a continuous casting process[END_REF] . For transient problems, such as the prediction of residual stress and shape (butt-curl) during startup of the aluminum DC continuous casting process, the domain can be extended in time by adding layers. [START_REF] Sengupta | Quantification of temperature, stress, and strain fields during the start=up phase of direct chill casting process by using a 3D fully coupled thermal and stress model for AA5182 ingots[END_REF] • A second strategy considers the entire casting, including the mushy and liquid regions. The liquid, mushy and solid regions are modeled as a continuum by adopting the constitutive equations for the solid phase (Section 2.4) for all regions by adjusting material parameters such as K, m, E, ν, σ 0 , and ρ according to temperature. For example, liquid can be treated by setting the strains to zero when the temperature is above the solidus temperature. This ensures that stress development in the liquid phase is suppressed. In the equilibrium equation, Eq. 2, acceleration terms are neglected, and a small-strain analysis can be performed. The primary unknowns are the displacements, or displacement increments. This popular approach can be used with structural finite element codes, such as MARC [15] or ABAQUS [START_REF]ABAQUS Theory Manual v6.0[END_REF] and with commercial solidification codes or special-purpose software, such as ALSIM [START_REF] Hakonsen | Modeling of Casting, Welding and Advanced Solidification Processes VII[END_REF] / ALSPEN, [START_REF] Fjaer | ALSPEN -A mathematical model for thermal stresses in DC-cast Al billets[END_REF] CASTS, [START_REF] Laschet | Thermo-mechanical analysis of cast/mould interaction in casting processes[END_REF] 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] Magmasoft, [START_REF] Magmasoft | [END_REF] and Procast [START_REF] Procast | [END_REF][START_REF] Samonds | Coupled thermal-fluid-stress analysis of castings[END_REF] It has been applied successfully to simulate deformation and residual stress in shape castings , [START_REF] Drezet | Thermomechanical Effects in D.C. Casting of Aluminum Alloy: A Numerical Benchmark Study[END_REF][START_REF] Fjaer | The Mechanism of Pull-in during DC-Casting of Aluminum Sheet Ingots[END_REF] DC casting of aluminum, [START_REF] Sengupta | Quantification of temperature, stress, and strain fields during the start=up phase of direct chill casting process by using a 3D fully coupled thermal and stress model for AA5182 ingots[END_REF][START_REF] Hakonsen | Modeling of Casting, Welding and Advanced Solidification Processes VII[END_REF][START_REF] Fjaer | ALSPEN -A mathematical model for thermal stresses in DC-cast Al billets[END_REF][START_REF] Mo | [END_REF][28][29] and continuous casting of steel [START_REF] Li | Thermo-Mechanical Finite-Element Model of Shell Behavior in Continuous Casting of Steel[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] Despite its efficiency, this approach may suffer from several drawbacks. First, it cannot properly account for fluid flow and the volumetric shrinkage that affects flow in the liquid pool, fluid feeding into the mushy zone, and primary shrinkage depressions that affect casting shape. In addition, incompressibility of the metal in the liquid state is accounted for by increasing Poisson's ratio close to 0.5, which sometimes makes the solution prone to numerical instability [START_REF] Zienkiewicz | The Finite Element Method[END_REF][START_REF] Drezet | Stress-strain computations of the DC casting process of aluminum alloy: a sensitivity study on material properties[END_REF] • A third strategy has been recently developed, which addresses the above issues. It still simulates the entire casting, but treats the mass and momentum equations of the liquid and mushy regions with a mixed velocity-pressure formulation. The primary unknowns are the velocity (time derivative of displacement) and pressure fields, which makes it easier to impose the incompressibility constraint (see section 4). Indeed, the velocity-pressure formulation is also applied to the equilibrium of the solid regions, in order to provide a single continuum framework for the global numerical solution. This strategy has been implemented into codes dedicated to casting analysis such as THERCAST, [START_REF] Bellet | An ALE-FEM approach to the thermomechanics of solidification processes with application to the prediction of pipe shrinkage[END_REF][START_REF] Bellet | ALE method for solidification modelling[END_REF]34] and VULCAN [START_REF] Chiumenti | Coupled thermomechanical simulation of solidification and cooling phases in casting processes[END_REF] . If stress prediction is not important so that elastic strains can be ignored, then this formulation simplifies to a standard fluid flow analysis, which is useful in the prediction of bulging and shape in large-strain processes.

Example of Solid-State 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-viscoplastic constitutive equation was developed for the austenite phase of steel [START_REF] Kozlowski | Simple Constitutive Equations for Steel at High Temperature[END_REF] by fitting constant strain-rate tensile tests [START_REF] Wray | Plastic Deformation of Delta-Ferritic Iron at Intermediate Strain Rates[END_REF][START_REF] Wray | Effect of Carbon Content on the Plastic Flow of Plain Carbon Steels at Elevated Temperatures[END_REF] and constant-load creep tests [START_REF] Suzuki | Creep Properties of Steel at Continuous Casting Temperatures[END_REF] to the form required in Eqs. ( 10)-( 13).

m eq C eq C f eq eq f T f C C T f T m T T ε σ σ σ ε σ σ - - -   × = - -     = × + × + × = - × = - + × = - × & ( 16
)
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 [START_REF] Puhringer | Strand Mechanics for Continuous Slab Casting Plants[END_REF] to 44 GPa. [START_REF] Hub | Measurement of Velocity and Attenuation of Sound in Iron up to the Melting Point[END_REF] with typical modulus values ~10 GPa near the solidus. [START_REF] Mizukami | Mechanical Properties of Continuously Cast Steels at High Temperatures[END_REF][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. ( 4), ( 8) or ( 13) 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 solid data for ferrite (α), austenite (γ), and delta (δ) [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 liquid (l) measurements. [START_REF] Cramb | The density of liquid iron-carbon alloys[END_REF]

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 α α γ γ δ δ α γ δ ρ ρ ρ ρ ρ ρ ρ ρ ρ - = + + + = - -× - = - + - = - + = - - - - (17)
Specialized experiments to measure mechanical properties for use in computational models will be an important trend for future research in this field.

Thermomechanical 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 illustrated in Figure 1 and discussed hereunder.

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: T and m T tending towards a unique interface temperature. The reality is somewhat different, showing always non perfect contact conditions. Therefore, the conductive heat exchange coefficient g k h gas cond = should be limited by a finite value 0 h , corresponding to the "nogap" situation, and depends on the roughness of the casting surface. 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 mold flux. [START_REF] Meng | Heat Transfer and Solidification Model of Continuous Slab Casting: CON1D[END_REF] 3.1.2 Effective contact : heat transfer as a function of contact pressure With 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 as the heat exchange coefficient corresponding to no gap and no contact pressure, the interfacial heat flux is:

) )( ( 0 m c B c contact T T Ap h q - + = ( 19
)
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

The thermal-mechanical modeling equations just presented must be solved numerically, owing to the complex shape of the casting process domain, and the highly nonlinear material properties. The calculation depends greatly on the numerical resolution of time and space. Although finite-difference approaches are popular for heat transfer, solidification, and fluid flow analyses, the finite element formulation is usually preferred for the mechanical analysis, owing to its historical advantages with unstructured meshes and accurate implicit solution of the resulting simultaneous algebraic equations. The latter are discussed below.

Finite element formulation and numerical implementation

In the framework of the small strain approach presented above (section 2.5), having displacements for primitive unknowns, the weak form of the equilibrium equation, Eq. ( 2), neglecting inertia terms, is written as:

0 ) ( d : * * * * = ⋅ - ⋅ - ∀ ∫ ∫ ∫ Ω Ω ∂ Ω dV dS dV u g u T u ε σ u ρ ( 20
)
where T is the external stress vector. The vector test-functions * u can be seen as virtual displacements in a statement of virtual work.

If the third strategy described in section 2.5 is adopted, with velocity and pressure as primary unknown variables, the weak form of the momentum equation (Eq. ( 2)) 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 & & ρ ρ (21)
The first equation, contains vector test-functions * v , which can be seen as virtual velocities in a statement of virtual power. Unlike Eq. ( 20), the pressure p is a primary variable, and only the deviatoric part of the constitutive equations is involved (to determine the stress deviator s). This is why the second equation is needed, which consists of a weak form of the incompressibility of inelastic deformations. Equations ( 20) and ( 21) are spatially discretized using the standard finite-element method, as explained in detail in many references. [START_REF] Rappaz | Numerical Modeling in Materials Science and Engineering[END_REF] Combined with time discretization using finite differences, this leads to a set of non-linear equations to be solved at each time increment. In the context of the displacement strategy, Eq. ( 20), this leads to:

0 ) ( = U R (22)
where R is the vector of the nodal equilibrium residues (number of components: 3 x number of nodes, in dimension 3), and U is the vector of nodal incremental displacements (same size).

Adopting the velocity-pressure strategy, Eq. ( 21), leads to a set of non-linear equations:

0 ) , ( ' = P V R (23)
where R' is the vector of the nodal residues (number of components: 4 x number of nodes, in dimension 3), V is the vector of nodal velocities (size: 3 x number of nodes) and P is the vector of nodal pressures (size: number of nodes).

The global finite-element systems (22) or ( 23) are usually solved using a full or modified Newton-Raphson method, [31, , 2004 #3715] which iterates to minimize the norm of the residue vectors R or R'. Alternatively, explicit methods may be employed at this global level.

At the local (finite element) level, an algorithm is also required to integrate the constitutive equations, when they depend on strain-rate or strain. When the constitutive equations are highly non-linear, an implicit algorithm is useful to perform time integration at each Gauss point in order to provide better estimates of inelastic strain at the local level. [START_REF] Nemat-Nasser | An explicit Algorithm for large-strain, large-strain rate elasticviscoplasticity[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] 4.2 Boundary conditions: modeling of contact conditions. Multidomain approaches At the interface between the solidifying material and the mold, a contact condition is required to prevent penetration of the shell into the mold, while allowing shrinkage of the shell away from the mold to create an interfacial gap:

:      = ⋅ ≥ ≤ ⋅ 0 ) ( 0 0 g g n σn
n σn (24) where g is the local interface gap width (positive when air gap exists effectively, as in section 3.1.1) and n is the local outward unit normal to the part. Eq. ( 24) can be satisfied with a penalty condition, which consists in applying a normal stress vector T proportional to the penetration depth (if any) via a penalty constant χ p :

n σn T g p - - = = χ (25)
Here again, the brackets denote the positive part: a repulsive stress is applied only if g is negative (penetration). Different methods of local adaptation of the penalty coefficient χ p have been developed, including 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 - - - = µ τ (26)
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. ( 25) 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. ( 18) and ( 19). [START_REF] Jaouen | A numerical mechanical coupling algorithm for deformable bodies: application to part/mold interaction in casting process[END_REF] 4.3 Treatment of the 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 2.4.

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 2.5), 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] 4.4 Thermomechanical coupling Because of the interdependency of the thermal and mechanical analyses, as presented in section 3, 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.

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 [START_REF] Weiner | Elasto-plastic thermal stresses in a solidifying body[END_REF] 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 estimating the discretization errors of computational thermalstress models, as it can be solved with a simple mesh consisting of one row elements, as shown in Figure 2. Numerical predictions should match with acceptable precision using the same element type and mesh refinement planned for the real problem. 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 typical conditions of steel casting. [START_REF] Koric | Efficient Thermo-Mechanical Model for Solidification Processes[END_REF] Figure 2. One-dimensional slice-domain for modeling solidifying plate.

Figure 3 and Figure 4 compare the temperature and stress profiles in the plate at two times. The temperature profile through the solidifying shell is almost linear. Because the interior cools relative to the fixed surface temperature, its shrinkage generates internal tensile stress, which induces compressive stress at the surface. With no applied external pressure, the average stress through the thickness must naturally equal zero, and stress must decrease to zero in the liquid. Stresses and strains in both transverse directions are equal for this symmetrical problem. The close agreement demonstrates that both computational models are numerically consistent and have an acceptable mesh resolution. Comparison with experimental measurements is also required, to validate that the modeling assumptions and input data are reasonable. 6 Example Applications

Sand Casting of Braking disks

The finite element software THERCAST for thermomechanical analysis of solidification [34] has been used in the automotive industry to predict distortion of grey-iron braking discs cast in sand molds. [START_REF] Bellet | Finite elements for a thermomechanical analysis of solidification processes[END_REF] 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. Figure 5 illustrates the discretization of the different domains involved in the calculation. The actual cooling scenario has been simulated: cooling in mold for 45 min, shake out and air cooling for 15 min. Figure 6 shows the temperature evolution at different points in a horizontal cross section at mid-height in the disc, revealing: solidification after 2 min, and solid state phase change after 20 min. The calculated deformation of the core ASM Handbook, Volume 15, Casting, Division 4: Modeling and Analysis of Casting Processes, American Society of Metals (2008) 449-461 13 blades shows thermal buckling due to the very high temperature, and constraint of their dilatation, as shown in Figure 7. 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 produce quality parts. Instead, process simulation allows the manufacturer to test alternative geometries and process conditions in order to minimize the defect. Similar thermomechanical calculations have been made for 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 8, calculations are consistent with measurements to within 10 MPa.

Continuous Casting of Steel Slabs

Thermomechanical simulations are used by steelmakers to analyze stresses and strains both in the mold and in the secondary cooling zone below. One goal is to quantify 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 internal cracks and macrosegregation [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, for the entire length of the caster using THERCAST, as described 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 were presented in section 2. Contact with supporting rolls is simulated with the penalty formulation discussed in section 4.2, adapting penalty coefficients for the different rolls continuously to control numerical penetration of the strand. shows results for a vertical-curved machine (strand thickness 0.22 m, casting speed 0.9 m/min, material Fe-0.06wt%C) at around 11 m below the meniscus. The pressure distribution reveals a double alternation of compressive and depressive zones. First, the strand surface is in a compressive state under the 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). Near the solidification front close to the solidus isotherm), the stress alternates between tension (negative pressure of about -2 MPa) beneath the rolls, and compression in between, (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 small 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 in Figure 10, 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] . Slab bulging calculated at two different casting speeds: 0.9 m/min and 1.2 m/min. The slab bulging increases with the casting speed. [START_REF] Triolet | Mastering steel slab internal soundness and surface quality issues through thermomechanical modelling of continuous casting[END_REF] Pressure Hot tearing analysis Hot tear crack formation 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, too small-scale, 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.

Thermal-analysis Criteria

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. One [START_REF] Clyne | Comparison between experimental data and theorical predictions relating to dependence of solidification cracking on composition[END_REF] 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:

Classical mechanics criteria

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[END_REF] accumulates inelastic deformation over a brittleness temperature range, which is defined, for example as (28) 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[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[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[END_REF] Tensile stress is also a requirement for hot tear formation. [START_REF] Yamanaka | Measurement of Critical Strain for Solidification Cracking[END_REF] The maximum tensile stress occurs just before formation of a critical flaw. [START_REF] Wisniewski | Tensile Behavior of Solidifying Aluminum Alloys[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 ε ε = ∆ & ( 29
)
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.

Mechanistically-based criteria

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 & ( 30
)
in which l µ is the dynamic viscosity of the liquid phase, 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: 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] 7.4 Case Study: Billet Casting Speed Optimization A Lagrangian model of temperature, distortion, strain, stress, hot tearing has been applied to predict the maximum speed for continuous-casting of steel billets without forming off-corner internal cracks. 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 used to track a transverse slice through the solidifying steel 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 (zdirection) 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 are also reasonable because bulging effects are small and the undiscretized casting direction is modeled with the appropriate condition of generalized plain strain. Other applications with this model include the prediction of ideal taper of the mold walls, [START_REF] Thomas | Ideal Taper Prediction for Slab Casting[END_REF] and quantifying the effect of steel grade on oscillation mark severity during level fluctuations [START_REF] Sengupta | Effect of a Sudden Level Fluctuation on Hook Formation During Continuous Casting of Ultra-Low Carbon Steel Slabs[END_REF] .

∫ ∇ =  ⌡ ⌠ =  ⌡ ⌠ = 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 R 2 1 2 1 2 1 ) ( ) (
The model domain is an L-shaped region of a 2-D transverse section, shown in Figure 11. Removing the central liquid region saves computation and lessens stability problems related to element "locking". Physically, this "trick" is important in two-dimensional domains because it allows the liquid volume to change without generating stress, which mimics the effect of fluid flow into and out of the domain that occurs in the actual open-topped casting process. Simulations start at the meniscus, 100 mm below the mold top, and extend through the 800-mm long mold and below, for a caster with no sub-mould support. The instantaneous heat flux, given in Eq. ( 32), was based on plant measurements. [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. Sample results are presented here for one-quarter of a 120 mm square billet cast at speeds of 2.0 and 5.0 m/min. The latter is the critical speed at which hot-tear crack failure of the shell is just predicted to occur. The temperature and axial (z) stress distributions in a typical section through the wideface of the steel shell cast at 2.0 m/min are shown in Figure 12 and Figure 13 at four different times during cooling in the mold. Unlike the analytical solution in Figure 3, the surface temperature drops as time progresses. The corresponding stress distributions are qualitatively similar to the analytical solution (Figure 4). The stresses increase with time, however, as solidification progresses. The realistic constitutive equations produce a large region of tension near the solidification front. The magnitude of these stresses (and the corresponding strains) are not predicted to be enough to cause hot tearing in the mold, however. The results from two different codes reasonably match, demonstrating that the formulations are accurately implemented, convergence has been achieved, and that the mesh and time-step refinement are sufficient. The shell becomes hotter and thinner with increasing casting speed, owing to less time in the mold. This weakens the shell, allowing it to bulge more under the ferrostatic pressure below the mold. Figure 14(b) shows contours of "hoop" stress constructed by taking the stress component acting perpendicular to the dendrite growth direction, which simplifies to σ x in the lower right portion of the domain and σ y in the upper left 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 offcorner 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 Figure 14(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. ( 28) with the critical strain given in Eq. ( 29). 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. ( 29). Figure 14(c) shows contours of hot-tear strain in the hoop direction. The highest values appear at the off-corner subsurface 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 offcorner 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 15. The bulged shape is also similar. 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 16. [START_REF] Li | Thermo-Mechanical Finite Element Model of Bulging and Hot Tearing During Continuous Casting of Steel Billets[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 provides 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. Comparison of critical casting speeds, based on hot-tearing criterion, and typical plant practice. [START_REF] Howard | Development of High Speed Casting[END_REF] 8 Conclusions

 - ≤  =  >   (32)
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.

  u ˆ are prescribed displacements on boundary surface portion u Ω ∂ , and T ) are boundary surface forces or "tractions" on portion T Ω ∂ . The next sections first present the equilibrium and compatibility equations, and then introduce constitutive equations for the different material states during solidification.

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

  l

  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

Figure 1 .

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

4

 4

  Gap width g (m) k gas ASM Handbook, Volume 15, Casting, Division 4: Modeling and Analysis of Casting Processes, American Society of Metals (2008) 449-461

Figure 3 .

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

Figure 4 .

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

Figure 5 .

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

Figure 6 .Figure 7 .

 67 Figure 6. Temperature evolution in the part at different points located in the indicated section.

Figure 8 .

 8 Figure 8. Residual hoop stresses (left) and radial stresses (right) in a radial section on as-cast plain discs made of grey iron. Top line: calculated values; bottom line: measured values.

Figure 9

 9 Figure9shows results for a vertical-curved machine (strand thickness 0.22 m, casting speed 0.9 m/min, material Fe-0.06wt%C) at around 11 m below the meniscus. The pressure distribution reveals a double alternation of compressive and depressive zones. First, the strand surface is in a compressive state under the 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). Near the solidification front close to the solidus isotherm), the stress alternates between tension (negative pressure of about -2 MPa) beneath the rolls, and compression in between, (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 small slab sections moving downstream between rolls and submitted to the metallurgical pressure onto the solidification front.

Figure 9 .

 9 Figure 9. Predictions in the middle of the secondary cooling zone, about 11 m below the meniscus. The finite element mesh, (top left) features a fine band of 20 mm. The pressure distribution (right) reveals alternating stress, including tension near the solidification front (the mushy zone is materialized by 20 lines separated by an interval l g ∆ = 0.05).

  Figure 10. Slab bulging calculated at two different casting speeds: 0.9 m/min and 1.2 m/min. The slab bulging increases with the casting speed.[START_REF] Triolet | Mastering steel slab internal soundness and surface quality issues through thermomechanical modelling of continuous casting[END_REF]

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

Figure 11 .

 11 Figure 11. Model domain.

Figure 12 .

 12 Figure 12. Temperature distribution along the solidifying slice in continuous casting mold.

Figure 13 .

 13 Figure 13. Lateral (y and z) stress distribution along the solidifying slice in continuous casting mold.

Figure 14 (

 14 Figure 14(a) shows the distorted temperature contours near the strand corner at 200 mm below the mold exit, for a casting speed of 5.0 m/min. The corner region is coldest, owing to two-dimensional cooling. The shell becomes hotter and thinner with increasing casting speed, owing to less time in the mold. This weakens the shell, allowing it to bulge more under the ferrostatic pressure below the mold.

Figure 14 .

 14 Figure 14. Distorted contours at 200 mm below mold exit.

Figure 15 .

 15 Figure 15. Off-corner internal crack in break-out shell from a 175 mm square bloom.

  Figure16. Comparison of critical casting speeds, based on hot-tearing criterion, and typical plant practice.[START_REF] Howard | Development of High Speed Casting[END_REF]

  7 °C liquidus temperature, and 1411.8 °C solidus temperature. Constitutive equation and properties are given in Sections 2.4 and 2.6.

	( q MW m /	2	)	( ) ( ) 0.504 5 0.2444 sec. 4.7556 sec. t t -	t t	1.0sec. 1.0sec.

& & & &

& & & & &

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 ArcelorMittal, Ascometal, Aubert et Duval, Industeel and PSA Peugeot-Citroën, for support of this work.


Clément Dombry
email: clement.dombry@math.univ-poitiers.fr

Ingemar Kaj
email: ikaj@math.uu.se

The on-off network traffic model under intermediate scaling

Keywords: on-off process, workload process, renewal process, intermediate scaling, fractional Poisson motion, fractional Brownian motion, Lévy motion, heavy tails, long range dependence AMS Subject classification. Primary: 60G22; Secondary: 60F05 60F17 60K05

The result provided in this paper helps complete a unified picture of the scaling behavior in heavy-tailed stochastic models for transmission of packet traffic on high-speed communication links. Popular models include infinite source Poisson models, models based on aggregated renewal sequences, and models built from aggregated on-off sources. The versions of these models with finite variance transmission rate share the following pattern: if the sources connect at a fast rate over time the cumulative statistical fluctuations are fractional Brownian motion, if the connection rate is slow the traffic fluctuations are described by a stable Lévy process, while the limiting fluctuations for the intermediate scaling regime are given by fractional Poisson motion.

Introduction

It is well-known that packet traffic on high-speed links exhibit data characteristics consistent with long-range dependence and self-similarity. To explain the possible mechanisms behind this behavior, various network traffic models have been developed where these features arise as heavy-tailed phenomena; see [START_REF] Resnick | Heavy-tail phenomena, Probabilistic and statistical modeling[END_REF] [START_REF] Resnick | Heavy-tail phenomena, Probabilistic and statistical modeling[END_REF]. A natural basis for modeling such systems, applied early on during these developments, is the view of packet traffic composed of a large number of aggregated streams where each source alternates between an active onstate transmitting data and an inactive off-state. The traffic streams generate on average a given mean-rate traffic, they have stationary increments and they are considered statistically independent. In particular, the transmission channel is able to accommodate peak-rate traffic corresponding to all sources being in the on-state. To capture in this model the strong positive dependence manifest in empirical trace data measurements, it is assumed that the duration of on-periods and/or off-periods are subject to heavy-tailed probability distributions. It is then interesting to analyze the workload of total traffic over time and understand the random fluctuations around its cumulative average. Our continued interest in these questions comes from the finding that several scaling regimes exist with disparate asymptotic limits.

The first result of the type we have in mind is Taqqu, Willinger, Sherman (1997) [START_REF] Taqqu | Proof of a fundamental result in self-similar traffic modeling[END_REF], which introduces a double limit technique. In this sequential scheme, if the on-off model is averaged first over the level of aggregation and then over time the resulting limit process is fractional Brownian motion. As the fundamental example of a Gaussian self-similar process with long-range dependence, this limit preserves the inherent long-range dependence of the original workload fluctuations. On the other hand, averaging first over time and then over the number of traffic sources the limit process is a stable Lévy process. This alternative scaling limit is again self-similar but lacks long memory since the increments are independent. Moreover, having infinite variance the limiting workload is itself heavy-tailed. In Mikosh, Resnick, Rootzén, Stegeman (2002), [START_REF] Mikosch | Is network traffic approximated by stable Lévy motion of fractional Brownian motion[END_REF], the double limits are replaced by a single scheme where instead the number of sources grows at a rate which is relative to time. Two limit regimes of fast growth and slow growth are identified and two limit results corresponding to these are established, where again fractional Brownian motion and stable Lévy motion appear as scaled limit processes of the centered on-off workload. The purpose of this paper is to show that an additional limit process, fractional Poisson motion, arises under an intermediate scheme which can be viewed as a balanced scaling between slow and fast growth. In this case the scale of time grows essentially as a power function of the number of traffic sources. As will be recalled, fractional Poisson motion does indeed provide a bridge between fractional Brownian motion and stable Lévy motion.

The intermediate limit scheme discussed here is indicated in [START_REF] Kaj | Stochastic Modeling in Broadband Communications Systems[END_REF] [START_REF] Kaj | Stochastic Modeling in Broadband Communications Systems[END_REF], and introduced in Gaigalas and Kaj (2003) [START_REF] Gaigalas | Convergence of scaled renewal processes and a packet arrival model[END_REF], where limit results are given for a different but related class of traffic models under three scaling regimes referred to as slow, intermediate and fast connection rate. The workload process is again the superposition of independent traffic streams with stationary increments but now each source generates packets according to a finite mean renewal counting process with heavy-tailed interrenewal cycle lengths. The link to the class of on-off models is that each pair of an on-period and a successive off-period forms a renewal cycle and the number of such on-off cycles generate a heavy-tailed renewal counting process. Moreover, if we associate with each renewal cycle a reward given by the length of its on-period and apply a suitable interpretation of partial rewards, then the corresponding renewal-reward process coincides with the on-off workload process.

To explain briefly the limit result in [START_REF] Gaigalas | Convergence of scaled renewal processes and a packet arrival model[END_REF] under intermediate connection rate, let (N i (t)) i be i.i.d. copies of a stationary renewal counting process associated with a sequence of interrenewal times of finite mean µ and a regularly varying tail function F (t) ∼ L(t)t -γ , characterized by an index γ, 1 < γ < 2, and a slowly varying function L. Let m → ∞ and a → ∞ in such a way that mL(a)/a γ-1 → µc γ-1 for some constant c > 0. Then the weak convergence holds, 1 a m i=1

(N i (at) - at µ ) =⇒ - 1 µ c Y γ (t/c),
where Y γ (t) is an almost surely continuous, positively skewed, non-Gaussian and non-stable random process, which is defined by a particular representation of the characteristic function of its finite-dimensional distributions. Additional properties of the limit process are obtained in Kaj (2005) [9] and Gaigalas (2006) [START_REF] Gaigalas | A Poisson bridge between fractional Brownian motion and stable Lévy motion[END_REF], where it is shown with two different methods that Y γ can be represented as a stochastic integral with respect to a Poisson measure N (dx, du)

on R × R + with intensity measure n(dx, du) = γdxu -γ-1 du. Indeed,

Y γ (t) = R×R + t 0 1 [x,x+u] (y) dy (N (dx, du) -γ dxu -γ-1 du), t ≥ 0, ([7] uses F (t) ∼ L(t)γ -1 t -γ , consequently n(dx, du) = dxu -γ-1 du).
We call this process fractional Poisson motion with Hurst index H = (3γ)/2 ∈ (1/2, 1). With

σ 2 γ = 2 (γ -1)(2 -γ)(3 -γ) = 1 2H(1 -H)(2H -1) = σ 2 H , (1)
B H (t) = 1 σ H R×R + t 0 1 [x,x+u] (y) dy M (dx, du),
where M (dx, du) is a Gaussian random measure on R × R + which is characterized by the control measure (3 -2H) dxu -2(2-H) du. The covariance functions of B H and P H coincide. The fast connection rate limit for the model of aggregated renewal processes applies if mL(a)/a γ-1 → ∞ and the slow connection rate limit if mL(a)/a γ-1 → 0. For suitable normalizing sequences, the limit processes under these assumptions are fractional Brownian motion with Hurst index H = (3γ)/2 in the case of fast growth and a stable Lévy process with self-similarity index 1/γ in the slow growth situation, see [START_REF] Gaigalas | Convergence of scaled renewal processes and a packet arrival model[END_REF].

A number of other models have been suggested for the flow of traffic in communication networks. The superposition of independent renewal-reward processes applies more generally to sources which attain random transmission rates at random times, and not merely switch between on and off. For a model where the length of a transmission cycle as well as the transmission rate during the cycle are allowed to be heavy-tailed, Levy and Taqqu (2000) [START_REF] Lévy | Renewal reward processes with heavy-tailed interrenewal times and heavy-tailed rewards[END_REF], [START_REF] Pipiras | The limit of a renewal-reward process with heavy-tailed rewards is not a linear fractional stable motion[END_REF] [START_REF] Pipiras | The limit of a renewal-reward process with heavy-tailed rewards is not a linear fractional stable motion[END_REF], and Pipiras, Taqqu and Levy (2004) [START_REF] Pipiras | Slow, fast and arbitrary growth conditions for renewal reward processes when the renewals and the rewards are heavy-tailed[END_REF], established results for slow and fast growth scaling analogous to those for the on-off model. In addition, they obtained as a fast growth scaling limit a stable, self-similar process with stationary but not independent increments, coined the telecom process. A further category of models for network traffic with long-range dependence over time starts from the assumption that long-lived traffic sessions arrive according to a Poisson process. The sessions carry workload which is transmitted either at fixed rate, at a random rate throughout the session, or at a randomly varying rate over the session length. Such models, called infinite source Poisson models, are widely accepted as realistic workload processes for internet traffic. Indeed, it is natural to assume that web flows on a non-congested backbone link are initiated according to a Poisson process while the duration of sessions and transmission rates are highly variable. The conditions under which slow, intermediate and fast scaling results exist and fractional Brownian motion, fractional Poisson motion, stable Lévy processes and telecom processes arise in the asymptotic limits are known in great detail for variants of the infinite source Poisson model, see [START_REF] Kaj | Convergence to fractional Brownian motion and to the Telecom process: the integral representation approach[END_REF] [START_REF] Kaj | Convergence to fractional Brownian motion and to the Telecom process: the integral representation approach[END_REF]. In [START_REF] Kaj | Convergence to fractional Brownian motion and to the Telecom process: the integral representation approach[END_REF], Y γ is called the intermediate telecom process. Mikosh and Samorodnitsky (2007) [START_REF] Mikosch | Scaling Limits for Cumulative Input Processes[END_REF] consider scaling limits for a general class of input processes, which includes as special cases the models already mentioned as well as other cumulative cluster-type processes. It is shown that fractional Brownian motion is a robust limit for a variety of models under fast growth conditions, whereas the slow growth behavior is more variable with a number of different stable processes arising in the limit.

Our current result completes the picture for the intermediate scaling regime, where neither of the mechanisms of fast or slow growth are predominant. In this case, where the system workload is under simultaneous influence of Gaussian and stable domains of attraction, we show that the fluctuations which build up in the on-off model are robust and again described by the fractional Poisson motion, parallel to what is known to be valid for infinite source Poisson and renewal-based traffic models. In the next section 2 we introduce properly both the on-off model and the renewal-based model to be used as an approximation and we state the relevant background results for these models. In section 3 we state the main result and give the structure of the proof. Section 4 is devoted to remaining and technical aspects of the proof.

The on-off model and background results

We begin by introducing the on-off model using similar notations as in [START_REF] Mikosch | Is network traffic approximated by stable Lévy motion of fractional Brownian motion[END_REF]. Let X on , X 1 , X 2 , . . . be i.i.d. non-negative random variables with distribution F on representing the lengths of onperiods. Similarly let Y off , Y 1 , Y 2 , . . . be i.i.d. non-negative random variables with distribution F off representing the lengths of off-periods. The X-and Y -sequences are supposed to be independent. For any distribution function F we write F = 1 -F for the right tail. We fix two parameters, α on and α off , such that

1 < α on < α off < 2, (2)
and assume that Fon (x) = x -αon L on (x) and Foff (x

) = x -α off L off (x), x → ∞, (3)
with L on , L off arbitrary functions slowly varying at infinity. Hence both distributions F on and F off have finite mean values µ on and µ off but their variances are infinite. Assumption (2) agrees with that of [START_REF] Mikosch | Is network traffic approximated by stable Lévy motion of fractional Brownian motion[END_REF]. However, thanks to a simple symmetry argument, we can also cover the case α on > α off . The case α on = α off , for which the on-off process is an alternating renewal process, falls outside of the class of processes we are able to study within the methodology developed here.

We consider the renewal sequence generated by alternating on-and off-periods. For the purpose of stationarity we introduce random variables (X 0 , Y 0 ) representing the initial onand off-periods as follows: let B, X eq on , Y eq off be independent random variables, independent of {X on , (X n ), Y off , (Y n )}, and such that B is Bernoulli with

P(B = 1) = 1 -P(B = 0) = µ on /µ,
and X eq on and Y eq off have distribution functions

F eq on (x) = 1 µ on x 0 Fon (s) ds and F eq off (x) = 1 µ off x 0
Foff (s) ds, respectively. Now, let X 0 = BX eq on and Y 0 = BY off + (1 -B)Y eq off .

Note that X 0 and Y 0 are conditionally independent given B but not independent. At time t = 0 the system starts in the on-state if B = 1 and in the off-state if B = 0. With this initial distribution, the alternating renewal sequence is stationary and the probability that the system is in the on-state at any time t is µ on /µ. Renewal events occur at the start of each on-period. Inter-renewal times are given by the independent sequence

Z i = X i + Y i , i ≥ 0,
where Z i has distribution F = F on * F off and mean µ = µ on + µ off for i ≥ 1, and Z 0 has distribution function

F eq (x) = 1 µ x 0 F (s) ds.
The renewal sequence (T n ) n≥1 with delay T 0 is defined by

T n = n-1 i=0 Z i ,
and we denote by N (t) the associated counting process

N (t) = n≥0 1 (0,t] (T n ).
Note that N (t) has stationary increments and expectation E[N (t)] = t/µ. Moreover, because of (2), the tail behavior of the inter-renewal times is given by

F (x) ∼ L on (x)x -αon , x → ∞, (4)
see Asmussen [START_REF] Asmussen | Ruin Probabilities[END_REF], Chapter IX, Corollary 1.11. The on-off input process is the indicator process for the on-state defined by

I(t) = 1 [0,X 0 ) (t) + n≥0 1 [Tn,Tn+X n+1 ) (t), t ≥ 0.
The source is in the on-state if I(t) = 1 and in the off-state if I(t) = 0. The input process I(t) is strictly stationary with mean

E[I(t)] = P(I(t) = 1) = µ on /µ.
The associated cumulative workload defined by

W (t) = t 0 I(s) ds, t ≥ 0
is a stationary increment process with mean E[W t ] = tµ on /µ.

Let (I j , W j , N j ) j≥1 denote i.i.d. copies of the input process I, the accumulative workload process W , and the renewal counting process N for the stationary on-off model. For m ≥ 1, consider a server fed by m independent on-off sources. We define the cumulative workload of the m-server system as the superposition process

W m (t) = m j=1 W j (t), t ≥ 0, m ≥ 1,
and the renewal-cycle counting process for m aggregated traffic sources by

N m (t) = m j=1 N j (t), t ≥ 0, m ≥ 1.
In this paper, we are mainly concerned with the asymptotic properties of the cumulative workload when the number of sources, m, increases and time t is rescaled by a factor a > 0. Thus, we consider the centered and rescaled process

W m (at) -matµ on /µ b(a, m) = 1 b(a, m) m j=1 at 0 (I j (s) - µ on µ ) ds, t ≥ 0,
where the renormalization b(a, m) will be precised in the sequel. The asymptotic is considered when both m → ∞ and a → ∞. The relative growth of m and a have a major impact on the limit. Let a = a m be the sequence governing the scaling of time and suppose a m → ∞ as m → ∞ (we will often omit the subscript m). Following the notation in [START_REF] Gaigalas | Convergence of scaled renewal processes and a packet arrival model[END_REF], we consider the following three scaling regimes:

• fast connection rate mL on (a)/a αon-1 → ∞;

(FCR)

• slow connection rate mL on (a)/a αon-1 → 0; (SCR)

• intermediate connection rate

mL on (a)/a αon-1 → µ c αon-1 , 0 < c < ∞. (ICR)
In [START_REF] Mikosch | Is network traffic approximated by stable Lévy motion of fractional Brownian motion[END_REF], the asymptotic behavior of the cumulative total workload is investigated under conditions (FCR) and (SCR).

Theorem 2.1 (Mikosch et al.) Recall assumptions (2) and (3).

• Under condition (FCR) and with the normalization b(a, m) = (a 3-αon L on (a)m) 1/2 , the following weak convergence of processes holds in the space of continuous functions on

R + : W m (at) -matµ on /µ b(a, m) =⇒ σ αon µ on µ 3/2 B H (t), t ≥ 0
where B H (t) is a standard fractional Brownian motion with index H = (3α on )/2.

• Under condition (SCR) and with the normalization b(a, m) = (1/ Fon ) ← (am) := inf{x ≥ 0 : Fon (x) ≤ 1/am}, then in the sense of convergence of finite dimensional distributions,

W m (at) -matµ on /µ b(a, m) f dd -→ σ 0 µ off µ 1+1/αon X αon,1,1 (t), t ≥ 0,
where X αon,1,1 (t) is a standard α on -stable Lévy motion totally skewed to the right, i.e. such that X αon,1,1 (1) ∼ S αon (1, 1, 0), and

σ αon 0 = Γ(2 -α on ) cos(πα on /2) 1 -α on .
The intermediate regime for renewal processes was investigated in [START_REF] Gaigalas | Convergence of scaled renewal processes and a packet arrival model[END_REF]. The formulation adopted here is given in [START_REF] Kaj | Limiting fractal random processes in heavy-tailed systems[END_REF], and is an immediate consequence of (4).

N m (at) -mat/µ a =⇒ - 1 µ σ αon c P H (t/c),
where σ αon is given in [START_REF] Asmussen | Ruin Probabilities[END_REF] and P H (t) is the standard fractional Poisson motion

P H (t) = 1 σ αon R×R + t 0 1 [x,x+u] (y) dy (N (dx, du) -dx α on u -αon-1 du) (5)
with Hurst index H = (3α on )/2.

3 Intermediate limit for the on-off model

In this section, we investigate the intermediate scaling limit for the on-off model. The following is our main result.

Theorem 3.1 Under condition (ICR) and with the normalization b(a, m) = a, the following convergence of processes holds in the space of continuous functions on R + :

W m (at) -matµ on /µ a =⇒ σ αon µ off µ c P H (t/c),
with σ αon in [START_REF] Asmussen | Ruin Probabilities[END_REF] and P H (t) the standard fractional Poisson motion in [START_REF] Mikosch | Scaling Limits for Cumulative Input Processes[END_REF].

Remarks 3.2

The fractional Poisson motion is not self-similar but does have a property of aggregate-similarity, introduced in [START_REF] Kaj | Limiting fractal random processes in heavy-tailed systems[END_REF], which allows for an interpretation of the scaling parameter c. Consider for each integer m ≥ 1 the sequence c m = m 1/(αon-1) . Then

c m P H (t/c m ) f dd = m i=1 P i H (t),
where P 1 H , P 2 H , . . . are i.i.d. copies of P H . Consider also the sequence c ′ m = m -1/(αon-1) . For any m,

m i=1 c ′ m P i H (t/c ′ m ) f dd = P H (t).
Hence, by tracing the limit process in Theorem 3.1 as c m → ∞, we recover in distribution the succession of all aggregates 1≤i≤m P i H , m ≥ 1. Also, by letting c ′ m → 0 we find that the limit process represents successively smaller fractions which sum up to recover fractional Poisson motion.

These relations explain the fact that fractional Poisson motion acts as a bridge between the stable Levy process and fractional Brownian motion. First, {c H P H (t/c)} converges weakly

to {B H (t)}, as c → ∞. Indeed, c H m P H (t/c m ) f dd = 1
√ m 1≤i≤m P i H (t) and the central limit Theorem yields the Gaussian limit as m → ∞. The required tightness property is shown in [START_REF] Gaigalas | A Poisson bridge between fractional Brownian motion and stable Lévy motion[END_REF]. Moreover, it is shown in [START_REF] Gaigalas | A Poisson bridge between fractional Brownian motion and stable Lévy motion[END_REF] that c 1/αon P H (t/c) converges in distribution as c → 0 to the α on -stable Levy process. To see that the limit must be α on -stable, take d = c • c ′ m for any c > 0. Then

c 1/αon P H (t/c) f dd = 1 m 1/αon m i=1 d 1/αon P i H (t/d), m ≥ 1,
and, assuming that the rescaled process (c 1/αon P H (t/c)) t≥0 converge to some non-trivial limit process L, we must have as c → 0 (and hence d → 0)

L(t) f dd = 1 m 1/αon m i=1 L i (t), m ≥ 1.
This indicates that the limit L must be α on -stable.

Heuristics of the proof of Theorem 3.1 To motivate that the limit process in the intermediate connection rate limit appears naturally, we discuss a decomposition of the centered on-off process based on its representation as a renewal-reward model. We first note that the single source cumulative workload has the form

W (t) = X 0 ∧ t + N (t) i=1 X i -(T N (t)-1 + X N (t) -t) + .
Similarly, focusing on off-periods rather than on-periods, we have

t -W (t) = Y 0 ∧ t + N (t) i=1 Y i -(T N (t) -t) ∧ Y N (t) .
The centered single source workload is therefore

W (t) - µ on µ t = -(t -W (t)) + µ off µ t = -µ off (N (t) -t/µ) - N (t) i=1 (Y i -µ off ) + R(t) with R(t) = (T N (t) -t) ∧ Y N (t) -Y 0 ∧ t.
Thus, for the workload of m sources,

W m (t) - µ on µ mt = -µ off (N m (t) -mt/µ) - m j=1 N j (t) i=1 (Y j i -µ off ) + m j=1 R j (t) (6)
using obvious notations. The balancing of terms under the scaling relation (ICR), makes it plausible that both terms

1 a m j=1 N j (at) i=1 (Y j i -µ off ), 1 a m j=1
R j (at) vanish in the scaling limit. This suggests asymptotically,

W m (at) -µ on mat/µ a ∼ -µ off N m (at) -mat/µ a , (7)
and so Theorem 2.2 would imply Theorem 3.1. In the next final section, we will compare rigorously the two processes in [START_REF] Gaigalas | Convergence of scaled renewal processes and a packet arrival model[END_REF].

4 Proof of Theorem 3.1

The proof of Theorem 3.1 relies on the following three lemmas:

W m (at) -matµ on /µ a f dd =⇒ σ αon µ off µ c P H (t/c),
is a consequence of the decomposition given in [START_REF] Gaigalas | A Poisson bridge between fractional Brownian motion and stable Lévy motion[END_REF]. By Lemma 4.3 the sequence is tight in the space of continuous functions on R + . Hence weak convergence holds in the space of continuous functions and Theorem 3.1 is proved.

Proof of Lemma 4.1. We construct an alternative representation of the random variable in the left hand side of (8). Define

N 1 (at) = inf{k ≥ 0; X 1 0 + k i=1 Z 1 i ≥ at}
and for j ≥ 2

N j (at) = inf{k ≥ 0; X j 0 + k i=1 Z 1 N j-1 (at)+i ≥ at}.
For m ≥ 1, let N m (at) = m j=1 N j (at). The random variables N j (at), j ≥ 1 are i.i.d and for each fixed t ≥ 0

1 a m j=1 N j (at) i=1 (Y j i -µ off ) and 1 a Nm(at) i=1 (Y 1 i -µ off ),
have the same distribution (note that the uni-dimensional marginal distributions are equal but not the multidimensional distributions). This representation will enable us to prove that under assumption (ICR), in the space of cád-lág functions on R + endowed with the Skorokhod topology, we have the convergence

1 a amu i=1 (Y 1 i -µ off ) u≥0 =⇒ 0. (9)
Moreover,

1 am Ñm (at) =⇒ t µ . (10)
Equations ( 9) and ( 10) together imply

1 a Ñm(at) i=1 (Y 1 i -µ off ) =⇒ 0
and this proves the lemma. Thus, it remains to prove ( 9) and [START_REF] Kaj | Convergence to fractional Brownian motion and to the Telecom process: the integral representation approach[END_REF].

To this aim, recall that the random variables Y 1 i , i ≥ 1, are i.i.d. with distribution such that the tail function Foff is regularly varying with index -α off . Hence there exists a regularly varying function L such that the centered and rescaled sum

1 (am) 1/α off L(am) amu i=1 (Y 1 i -µ off ) u≥0
converges in the space of cád-lág functions to some α off -stable Lévy process (see [START_REF] Petrov | Sums of Independent Random Variables[END_REF], the exact form of L or of the limit process are not needed here). This implies the convergence property (9), since a >> (am) 1/α off L(am) under the scaling assumption (ICR) with α on < α off .

We now prove equation [START_REF] Kaj | Convergence to fractional Brownian motion and to the Telecom process: the integral representation approach[END_REF]. The stationary renewal process N (t) has mean t/µ and variance given asymptotically by Var(N (t)) ∼ σ 2 αon 1 µ 3 t 3-αon L on (t), t → ∞, see [START_REF] Gaigalas | Convergence of scaled renewal processes and a packet arrival model[END_REF], Equation (30), and references therein. Hence, 1 am N m (at) has mean t/µ and variance under scaling (ICR), such that

Var 1 am N m (at) = 1 a 2 m Var[N (at)] ∼ a 1-αon L on (at) m σ 2 αon 1 µ 3 t 3-αon → 0.
This shows that 1 am N m (at) converges in distribution to t/µ, which is [START_REF] Kaj | Convergence to fractional Brownian motion and to the Telecom process: the integral representation approach[END_REF]. This ends the proof of 4.1

Proof of Lemma 4.2. Since |R j (t)| ≤ Y j 0 + (T j N j (t) -t) ∧ Y j N j (t) , it is enough to prove that 1 a m j=1 Y j 0 =⇒ 0 and 1 a m j=1 (T j N j t -t) ∧ Y j N j (t) =⇒ 0.
By stationarity, the random variables Y j 0 and (T j N j t t) ∧ Y j N j (t) have the same distribution; they represent the remaining time after 0 and t, respectively, of the first off-period after 0, and after t. Since both sums have the same distribution, we only consider the first one. Using Karamata's Theorem (see [START_REF] Bingham | Regular Variation[END_REF]), the tail function F eq off satisfies

F eq off (x) = 1 µ off ∞ x Foff (s)ds ∼ 1 µ off x -(α off -1) α off -1 L off (x)
as x → ∞. This implies that the random variable Y 0 has a regularly varying tail with index -(α off -1) and hence belongs to the domain of attraction of an (α off -1)-stable distribution. Therefore there exists a slowly varying function L, such that Proof of Lemma 4.3. The proof given in [START_REF] Mikosch | Is network traffic approximated by stable Lévy motion of fractional Brownian motion[END_REF] for fast scaling (FCR) can be adapted to our settings. We recall only the main lines. According to Billingsley [START_REF] Billingsley | Convergence of probability measures[END_REF], Theorem 12.3, it is enough to prove that for any t 1 , t 2 with |t 1t 2 | ≤ 1 and for some ε > 0, there exists a constant C > 0 and an a 0 > 0, such that for all a ≥ a 0 Using the definition of W m , centering and stationarity, it is enough to prove that for all t ∈ [0, 1] and a ≥ a 0 , m a 2 Var[W at ] ≤ Ct 1+ε

1 m 1/(α off -1) L(m)
(the constant C may change from one appearance to another). However, according to [START_REF] Mikosch | Is network traffic approximated by stable Lévy motion of fractional Brownian motion[END_REF], Equation By [START_REF] Petrov | Sums of Independent Random Variables[END_REF], the function a → Var[W a ] is regularly varying with index 3α on . Then, using Potter bounds (see [START_REF] Bingham | Regular Variation[END_REF]), we conclude that there exist a 0 > 0 and ε < 1α on /2, such that for all t ∈ (0, 1) and a ≥ a 0 /t, Var[W at ] Var[W a ] ≤ 1 1ε t 3-αon-ε .

(see the proof of Lemma 13 in [START_REF] Mikosch | Is network traffic approximated by stable Lévy motion of fractional Brownian motion[END_REF] for details). This implies that for all t ∈ (0, 1) and all a such that at ≥ a 0 , m a 2 Var[W at ] ≤ C 1ε t 3-αon-ε ≤ Ct 1+ε .

On the other hand, if t ≤ a 0 /a, then, for a large enough, In the last inequality, we use the fact that 2α onε > 0 and so a 2-αon-ε L on (a) → ∞ as a → ∞; taking a 0 large enough, we can suppose that for a ≥ a 0 , a 2-αon-ε L on (a) remains bounded away from zero.

By combining the estimates for the cases at ≥ a 0 and at ≤ a 0 we obtain [START_REF] Lévy | Renewal reward processes with heavy-tailed interrenewal times and heavy-tailed rewards[END_REF], which completes the proof.

Theorem 2 . 2 (

 22 Gaigalas-Kaj) Under condition (ICR) and with the normalization b(a, m) = a, the following convergence of processes holds:

Lemma 4 . 1 8 ) 4 . 2 Lemma 4 . 3

 4184243 In the scaling (ICR), for all t ≥ 0, iµ off ) =⇒ 0. (Lemma In the scaling (ICR), for all t ≥ 0, In the scaling (ICR), the sequence of processesW m (at)matµ on /µ a , t ≥ 0, m ≥ 1is tight in the space of continuous functions on R + .Proof of Theorem 3.1. By Theorem 2.2, Lemma 4.1 and Lemma 4.2, the convergence of finite-dimensional distributions,

  to a stable law of index α off -1. Under scaling (ICR) with α on < α off , we have a >> m 1/(α off -1) L(m)

  m (at 2 )mat 2 µ on /µ) -(W m (at 1 )mat 1 µ on /µ)| 2 ≤ C|t 2t 1 | 1+ε .

) ∼ σ 2 αon µ 2 on µ 3 t 3 - 2 αon µ 2 µ 2 on 2

 23222 αon L on (t), t → ∞.[START_REF] Petrov | Sums of Independent Random Variables[END_REF] This relation and the scaling (ICR) together imply, as a → ∞,a 2 m ∼ c 1-αon µ a 3-αon L on (a) ∼ c 1-αon σ Var[W a ],and so there is C > 0, such that for a large enoughm a Var[W at ] ≤ C Var[W at ] Var[W a ] .

m a 2 a 3 -so m a 2

 232 Var[W at ] ≤ Ca 2 t 2 Var[W a ] ≤ C a 2 t 2 a 3-αon L on (a) ≤ C (at) 1+ε a 1-ε 0 αon L on (a)andVar[W at ] ≤ C a 1-ε 0 t 1+ε a 2-αon-ε L on (a) ≤ C t 1+ε .