E Frénod

Existence result for a model of Proteus mirabilis swarm *

In this paper we present a modification of the usual Proteus mirabilis Swarm model. For the obtained model (which is a two phase model with a non-linear diffusion term containing memory) we set up a collection of a priori estimates. Those estimates allow to get an existence and uniqueness result.

Introduction and results

Proteus mirabilis is a bacterium that can be either a short cell we call "swimmer" or an elongated cell capable of translocation we call "swarmer". A model of behaviour of Proteus mirabilis colony has been proposed by Esipov and Shapiro [START_REF] Esipov | Kinetic model of Proteus mirabilis swarm colony development[END_REF], based on ideas of Gurtin [START_REF] Gurtin | A system of equations for age-dependent population diffusion[END_REF]. In this paper, we prove an existence result to a model which is, in a way, a generalization but also a regularization of the Esipov and Shapiro [START_REF] Esipov | Kinetic model of Proteus mirabilis swarm colony development[END_REF] model.

The model under consideration here is a two phase model with a non-linear diffusion term containing memory for one of the two phases. It involves two functions ρ and Q. The function ρ = ρ(t, a, x) refers, at time t ∈ [0, T ), 0 < T < +∞, to the density of swarmers of age a ∈ [0, A), 0 < A ≤ +∞ at position x ∈ Ω, where Ω is a regular sub-domain of R 2 , with boundary ∂Ω. In each point x of ∂Ω, -→ ν = -→ ν (x) stands for the unit normal vector pointing outside Ω. The function Q = Q(t, x) stands, at time t, for the biomass density of swimmers on Ω. For a constant τ , those two functions are supposed to satisfy the following system:

∂Q ∂t = 1 -ξ τ Q + A 0 ρ(•, a,
•)e a/τ µ(•, a, •) da + χ(A)ρ(•, A, •)e A/τ , on [0, T ) × Ω, (1.1)

∂ρ ∂t + ∂ρ ∂a = -µρ + ∇ • D(M, Q, P ) + d ∇ρ , on [0, T ) × [0, A) × Ω, (1.2)
ρ(•, 0, •) = ξ τ Q, on [0, T ) × Ω, (1.3)
ρ(0, •, •) = ρ 0 , on [0, A) × Ω, (1.4) Q(0, •) = Q 0 , on Ω, (1.5) ∂ρ ∂ -→ ν = 0, on [0, T ) × [0, A) × ∂Ω. (1.6)
Above, χ(A) is an artifice allowing to take into account a possible maximum age A beyond which swarmers cannot exist. It has the following definition:

χ(A) = 1, ∀A ∈ R, χ(+∞) = 0. (1.7)
Denoting by C k b the space of functions having continuous and bounded derivatives up to order k, µ = µ(t, a, x) is a function such that µ ∈ C 2 b ([0, T ) × [0, A) × Ω), µ ≥ 0, lim a→A µ(t, a, x) = µ uniformly in x and t, (

with µ ≥ c(1 -χ(A)), for a constant c > 0. The function ξ = ξ(t, Q) satisfies

ξ ∈ C 2 b ([0, T ) × R), 0 ≤ ξ ≤ 1.
(1.9)

In the second equation, ∇ stands for the gradient with respect to the x-variable, and ∇• for the divergence. The diffusion coefficient is the sum of a constant, a priori small, d > 0, and of D = D(M, P, Q), a non negative C 1 b function of its arguments.

(1.10)

In (1.10), P is defined for 0 ≤ a min < A by P (t, x) = A amin ρ(t, a, x)e a/τ da, (1.11) and Q is given by the first equation of the system. The memory (or hysteresis) term M = M [P ] (t, x) keeps information on the value of P in the past.

For four thresholds, P min < p min < p max < P max , with P min close to p min and P max close to p max , M is defined as the solution to:

∂M ∂t = 1 P max -p max H r P -p max P max -p max H r (1 -M) - 1 p min -P min H r p min -P p min -P min H r (M), M(0, •) = M 0 , (1.12)
with, denoting P 0 = P (0, •), M 0 ∈ C 1 b (Ω), 0 ≤ M 0 ≤ 1, M 0 = 0 where P 0 < P min and M 0 = 1 where P 0 > P max , (1. [START_REF] Huang | An age-dependent population model with nonlinear diffusion in R n[END_REF] and with

H r (p) = 0 if p ≤ 0, H r (p) = p if 0 ≤ p ≤ 1 and H r (p) = 1 if p ≥ 1. (1.14)
We now turn to the statement of the main result of this paper.

Theorem 1.1 Under assumptions (1.8), (1.9), (1.10) and (1.13), if

ρ 0 ≥ 0 ∈ L 1 ∩W 2,2 ∩W 1,4 ([0, A)× Ω, e a/τ dadx) satisfies Ω (ρ 0 ) 2 e 2a/τ dx ≤ b Ω (ρ 0 ) 2 e a/τ dx, ∀a ∈ [0, A), A 0 Ω |∇ρ 0 | 4 e 4a/τ dadx ≤ b A 0 Ω |∇ρ 0 | 4 e a/τ dadx, A 0 Ω |∇ρ 0 | 2 e 2a/τ dx ≤ b A 0 Ω |∇ρ 0 | 2 e a/τ dx, (1.15)
for a constant b, and if

Q 0 ≥ 0 ∈ L 1 ∩ W 2,2 ∩ W 1,4
(Ω); then, there exists a unique solution

(Q, ρ) ∈ L ∞ 0, T ; ((L 1 ∩W 1,2 (Ω))× (L 1 ∩W 1,2 ([0, A)× Ω, e a/τ dadx)) to system (1.1)-(1.6) coupled
with (1.11) and (1.12). Moreover, Q ≥ 0 and ρ ≥ 0.

The precise definitions of the spaces at work in the Theorem are given in the beginning of section 4.

We now give references where modelling and mathematical methods are developed on agestructured population problem: Gurtin and Mac Camy [START_REF] Gurtin | Non-linear age-dependent population dynamics[END_REF], Marcati [START_REF] Marcati | Asymptotic behavior in age-dependent population dynamics with heredity renewal law[END_REF], Andreasen [START_REF] Andreasen | Diesease regulation of age-structured host populations[END_REF][START_REF] Andreasen | The effect of age-dependent host mortality on the dynamics of an endemic disease[END_REF][START_REF] Andreasen | Instability in SIR-model with age-dependent susceptibility[END_REF]; possibly with diffusion: Gurtin [START_REF] Gurtin | A system of equations for age-dependent population diffusion[END_REF], Di Blasio and Lamberti [START_REF] Blasio | An initial-boundary value problem for age-dependent population diffusion[END_REF], Di Blasio [START_REF] Blasio | Non linear age-dependent population diffusion[END_REF], Mac Camy [START_REF] Mac Camy | A population model with nonlinear diffusion[END_REF], Gurtin and Mac Camy [START_REF] Gurtin | Diffusion model for age-structured population[END_REF], Busenberg and Iannelli [START_REF] Busenberg | A class of nonlinear diffusion problems in age-dependent population dynamics[END_REF], Langlais [START_REF] Langlais | A nonlinear problem in age-dependent population diffusion[END_REF][START_REF] Langlais | Large time behavior in a nonlinear age-dependent population dynamics problem with spatial diffusion[END_REF], Kubo and Langlais [START_REF] Kubo | Periodic solutions for population dynamics problem with agedependent and spatial structure[END_REF], Huang [START_REF] Huang | An age-dependent population model with nonlinear diffusion in R n[END_REF] and Esipov and Shapiro [START_REF] Esipov | Kinetic model of Proteus mirabilis swarm colony development[END_REF]. For simulation methods we refer to Lopez and Trigiante [START_REF] Lopez | A finite difference scheme for stiff problem arising in the numerical solution of a population dynamic model with spatial diffusion[END_REF], Milner [START_REF] Milner | A numerical method for a model of population dynamics with spatial diffusion[END_REF] Kim [START_REF] Kim | Galerkin methods for a model of population dynamics with nonlinear diffusion[END_REF], Esipov and Shapiro [START_REF] Esipov | Kinetic model of Proteus mirabilis swarm colony development[END_REF], Medvedev, Kapper and Koppel [START_REF] Medvedev | A reaction-diffusion system with periodic front dynamics[END_REF], Ayati and Dupont [START_REF] Ayati | Galerkin methods in age and space for a population model with nonlinear diffusion[END_REF]. Concerning the biological description of Proteus mirabilis colony behaviour, we refer for instance to Rauprich et al [START_REF] Rauprich | Periodic phenomena in Proteus mirabilis swarm colony development[END_REF], Gué, Dupont, Dufour and Sire [START_REF] Gué | Bacterial swarming: A biological time-resolved FTIR-ATR study of Proteus mirabilis swarm-cell differentiation[END_REF] and theirs references.

The paper is organized as follows: in section 2 we present the way to go from the Esipov and Shapiro model to system (1.1)- (1.6). Then section 3 is devoted to a priori estimates for the solution to (1.1)- (1.6). By a usual procedure consisting in linearizing and passing to limit, we prove the Theorem in section 4.

Acknowledgements: I would like to thank O. Sire for having introduced to me the swarm model of Esipov and Shapiro [START_REF] Esipov | Kinetic model of Proteus mirabilis swarm colony development[END_REF] and for stimulating discussions. I would also like to thank F. Granger who, despite stopping its PhD thesis for personal reasons, made the first steps towards the result.

Model

Proteus mirabilis is a pathogenic bacterium of urinary tract that when standing in liquid medium, consists in a usual short "swimmer cell" or "swimmer". When placed on agar medium, if the bacterium density is large enough, it may undergo a differentiation process producing an elongated cell with several nuclei called "swarmer cell" or "swarmer". Those swarmers are capable of translocation allowing the bacterial colony to colonise the medium.

The macroscopic model built by Esipov and Shapiro [START_REF] Esipov | Kinetic model of Proteus mirabilis swarm colony development[END_REF] describes this swarm phenomenon at the colony scale. We shall explain this model now. The swarmer behaviour depends on their own age. This dependence is taken into account by introducing the age dependent density of swarmers ρ(t, a, x). The link between this density ρ and the biomass density is a consequence of the fact that the mass of each cell is in direct proportion with its length and that the length growth of a swarmer is also in direct proportion with its length. Then the biomass density, at time t, of swarmers of age a at position x is ρ(t, a, x)e a/τ , where τ is the growth rate of the biomass. The first age-depending behaviour of the swarmers is that they actively participate in group migration only after an age a min . Then the definition of the biomass density P of swarmers capable of active translocation is given by (1.11). The second age-depending behaviour is that the swarmers dedifferentiate themselves and give swimmers. On this topic, Esipov and Shapiro [START_REF] Esipov | Kinetic model of Proteus mirabilis swarm colony development[END_REF] consider two situations. In the first one (Model A), the swarmers dedifferentiate themselves at a given age a max . The second situation consider that swarmers may dedifferentiate at each time with a probability 1/a (Model B). Now, we are able to write the swimmer evolution equation for the biomass density Q. Its evolution results from the classical cellular division, with a characteristic time which is the biomass growth rate τ , subtracting the proportion of bacteria undergoing differentiation and adding the dedifferentiation product. In the case of Model A, the evolution equation for Q is then

∂Q ∂t = 1 -ξ τ Q + ρ(•, a max , •)e amax/τ , Q(0, •) = Q 0 . (2.1)
Here, Q 0 stands for the initial swimmer density and ξ/τ is the fraction of swimmer population to produce swarmers. In the case of Model B, the evolution equation is

∂Q ∂t = 1 -ξ τ Q + t 0 ρ(•, a, •) e a/τ a da, Q(0, •) = Q 0 . (2.2)
We turn to the evolution of the swarmer density ρ. Its evolution is linked with ageing and the dedifferentiation process, but also to swarm. This last phenomenon is modelled by a non linear diffusion term with memory. The evolution equation for ρ is then

∂ρ ∂t + ∂ρ ∂a = ∇ • D(M, Q, P )∇ρ , (2.3)
in the case of Model A; and, in the case of Model B, it is

∂ρ ∂t + ∂ρ ∂a = - 1 a ρ + ∇ • D(M, Q, P )∇ρ . (2.4)
Both of those equations are equiped with the following initial and boundary conditions:

(a): ρ(•, 0, •) = ξ τ Q, (b): ρ(0, •, •) = 0, and (c): ∂ρ ∂ -→ ν ∂Ω = 0. (2.5)
The first of those three conditions means the fraction ξ/τ of swimmers undergoing the differentiation process produces swarmers of age 0. The initial condition on ρ means that, at the beginning of the process, there is no swarmer. The boundary condition means that swarmer cannot leave the domain Ω. In the diffusion term ∇ • D(M, Q, P )∇ρ appearing in (2.3) and (2.4), and modelling the swarm, the diffusion factor D(M, Q, P ) depends on the present value of Q and P but also on the history of P .

The term M then keeps in memory informations concerning the history of the swarmer density. Esipov and Shapiro [START_REF] Esipov | Kinetic model of Proteus mirabilis swarm colony development[END_REF] defines M as being set to 1 in a given point x if the value of P in x reaches a threshold P max . Then it remains at the value 1 until the value of P in x reaches another value P min < P max . Then they suggest to take

D(M, Q, P ) = D 0 M γ P P max exp -Q Q sat , (2.6)
for given values of D 0 and Q sat and with, H being such that H(p) = 0 if p < 0 and

H(p) = 1 if p > 0, γ(p) = (p - P min P max )H(p - P min P max ) or γ(p) = p - P min P max or γ(p) = p 2 or γ(p) = 1. (2.7)
We also mention that Medvedev, Kapper and Koppel [START_REF] Medvedev | A reaction-diffusion system with periodic front dynamics[END_REF] studied this model taking, for a given value of k, ) is a generalization and a regularization of the Esipov and Shapiro [START_REF] Esipov | Kinetic model of Proteus mirabilis swarm colony development[END_REF] model. First, because of the initial condition (2.5.b), it is an easy game to see that, at least formally, the solution ρ to (2.3) or (2.4) satisfies ρ(t, a, x) = 0 when a > t.

D(M, Q, P ) = D(Q, P ) = D 0 P P + kQ . ( 2
(2.9)

Hence we can replace (2.2) by

∂Q ∂t = 1 -ξ τ Q + +∞ 0 ρ(•, a, •) e a/τ a da, Q(0, •) = Q 0 . (2.10)
Making this allows to take under consideration, with no loss of consistency, initial data ρ(0, •, •) that are not 0 coming to the more general initial and boundary conditions (1.3)-(1.6). Secondly, equations (2.1) and (2.10) are particular cases of the general equation (1.1) with assumption (1.8). The case of equation (2.1) is recovered setting µ = 0 and A = a max and the case of (2.10) is recovered setting A = +∞ and µ = 1/a. We also see that (1.2) with d = 0 is a general framework inside which (2.3) and (2.4) may enter directly.

Concerning the fraction ξ of swimmers to produce swarmers, it seems to depend on experimental conditions and to be 0 when the swimmer density is high. It seems then to be reasonable to set that ξ satisfies (1.9).

The first regularization effect we consider in our model (1.1)-(1.6) consists in adding d > 0 in (1.2). This may be justified by experimental arguments saying that swarmers always experience a small random motion even before and after swarming.

The second regularization effect which constitutes the most visible modification of the model concerns the memory term M. In order to explain the way to go from the definition of M by Esipov and Shapiro [START_REF] Esipov | Kinetic model of Proteus mirabilis swarm colony development[END_REF] to the definition (1.12), we first notice that the term M of Esipov and Shapiro could be defined formally, in any x, as the solution to the following equation:

∂M ∂t (•, x) = (1 -M(•, x))δ {t/P (t,x)=Pmax} -M(•, x)δ {t/P (t,x)=Pmin} , M(0, x) = 0, (2.11)
where δ {t/P (t,x)=Pmax} stands for the Dirac measure in instant t where P (t, x) = P max . Of course this equation has no real mathematical meaning. But formally, if M = 0 and if the value P max is reached at a given time t, the solution of (2.11) experiences a jump +1. When P = P max , and M = 1, nothing happens and this is what is needed. In the same way, if at a given time t, P = P min , nothing happens if M = 0 and M experiences a jump -1 if M = 1. Now, it is easy to see that the right hang side of (1.12) is nothing but a regularization of the right hand value of (2.11). This is the reason why we make this choice to define M.

A priori estimates

The key point to get the existence result is a collection of a priori estimates satisfied by (Q, ρ). They are mathematical translations of biological properties. In order to set those estimates comfortably, we assume that the assumptions and the conclusions of Theorem 1.1 are satisfied.

L 1 and L 2 estimates

For p > 0, we denote

ρ(t) p = A 0 Ω |ρ(t, a, x)| p e a/τ dxda 1/p , Q(t) p = Ω |Q(t, x)| p dx 1/p , (3.1)
and we have the following estimates saying that the total biomass grows exponentially with a growth rate τ .

Lemma 3.1 If the assumptions of Theorem 1.1 are valid, and if the solution (Q, ρ) given by this same Theorem exists, then it satisfies

ρ(t) 1 + Q(t) 1 = ρ 0 1 + Q 0 1 e t/τ . (3.2)
Proof. Multiplying equation (1.2) by e a/τ , integrating then with respect to x and a yields:

d ρ 1 dt + lim a→A Ω ρ e a/τ dx - Ω ρ(•, 0, •) dx = 1 τ ρ 1 - A 0 Ω µρ e a/τ dxda, (3.3)
which also reads

d ρ 1 dt + χ(A)e A/τ Ω ρ(•, A, •) dx - Ω ξ τ Q dx = 1 τ ρ 1 - A 0 Ω
µρ e a/τ dxda.

(3.4)

Integrating now (1.1) with respect to x, we obtain

d Q 1 dt = Ω 1 -ξ τ Q dx + A 0 Ω µρ e a/τ dxda + χ(A)e A/τ Ω ρ(•, A, •) dx. (3.5)
Summing up (3.4) and (3.5) finally gives

d ρ 1 + Q 1 dt = 1 τ ρ 1 + Q 1 , (3.6)
proving the Lemma.

The second Lemma concerns the L 2 -norms of ρ and Q. It mathematically translates that the biomass cannot be so gathered that a null area set contains a positive biomass quantity. Lemma 3.2 If the assumptions of Theorem 1.1 are realized and if the solution (Q, ρ) exists, then for any t ∈ [0, T ), it satisfies

ρ(t) 2 + Q(t) 2 ≤ c ρ 0 2 + Q 0 2 , (3.7)
for a constant c (c only depends on A, τ , T and μ = sup{µ(t, a, x), t ∈ [0, T ), a ∈ [0, A), x ∈ Ω}.)

Proof. First, integrating (1.2) with respect to x, we get:

∂ ∂t + ∂ ∂a Ω ρ dx = - Ω µρ dx ≤ 0. (3.8)
Making the same, after multiplying (1.2) by ρ, gives:

1 2 ∂ ∂t + ∂ ∂a Ω ρ 2 dx = - Ω µρ 2 dx - Ω (D + d)|∇ρ| 2 dx ≤ 0. (3.9)
Defining, for a fixed α and for p = 1 or 2,

r p : s → Ω ρ p (s, α + s, x) dx, (3.10)
we have

r ′ p (s) = ∂ ∂t + ∂ ∂a Ω ρ p dx (s, α + s) ≤ 0. (3.11)
Now, since for fixed t and a and with α = a -t we have Ω ρ p (t, a, x) dx = r p (t), if t < a the relation

r p (t) ≤ r p (0) = Ω ρ p (0, a -t, x) dx reads Ω ρ(t, a, x) dx ≤ Ω ρ 0 (a -t, x) dx, Ω ρ 2 (t, a, x) dx ≤ Ω ρ 2 0 (a -t, x) dx. (3.12)
In the case when t > a, the relation

r p (t) ≤ r p (t -a) = Ω ρ p (t -a, 0, x) dx gives Ω ρ(t, a, x) dx ≤ Ω ξ τ Q(t -a, x) dxda, Ω ρ 2 (t, a, x) dx ≤ Ω ξ 2 τ 2 Q 2 (t -a, x) dxda, (3.13)
Secondly, we multiply equation (1.2) by e a/τ ρ and we integrate in x and a. Since

and since

Ω ρ 2 (t, 0, x) dx = Ω ξ 2 τ 2 Q 2 dx, (3.15)
we obtain

d ρ 2 2 dt + χ(A)e A/τ Ω ρ 2 (•, A, •) dx + 2 A 0 Ω (D + d)|∇ρ| 2 e a/τ dxda + 2 A 0 Ω µρ 2 e a/τ dxda = Ω ξ 2 τ 2 Q 2 dx + 1 τ ρ 2 2 .
(3.16) As the second, third and fourth terms of the left hand side of equation (3.16) are non negative, we may deduce

d ρ 2 2 dt ≤ 1 τ 2 Q 2 2 + 1 τ ρ 2 2 . (3.17)
In the same way, multiplying now (1.1) by Q and integrating yields

d Q 2 2 dt = 2 Ω 1 -ξ τ Q 2 dx + 2 A 0 Ω µρQ e a/τ dxda + 2χ(A)e A/τ Ω ρ(•, A, •)Q dx. (3.18)
Concerning the second term of the right hand side of the last equality, using Young's inequality and formula (3.12) and (3.13), we get

A 0 Ω µρQ e a/τ dxda ≤ A 0 Ω µ 2 ρ 2 dx 1/2 Ω Q 2 dx 1/2 e a/τ da ≤ μ t 0 Ω ξ 2 τ 2 Q 2 (t -a, •) dx 1/2 e a/τ da + H(A -t)e t/τ A t Ω ρ 2 0 (a -t, •)e 2(a-t)/τ dx 1/2 da Q 2 ≤ μ max(e T /τ , e T /τ τ 2 ) t 0 Ω Q 2 (t -a, •) dx 1/2 da + H(A -t) A t Ω ρ 2 0 (a -t, •)e 2(a-t)/τ dx 1/2 da Q 2 , (3.19)
where H(a) = 0 if a < 0 and H(a) = 1 if a ∈ [0, +∞] and where μ = sup{µ(t, a, x), t ∈ [0, T ), a ∈ [0, A), x ∈ Ω}. The third term of the right hand side of (3.18) may also be estimated:

χ(A)e A/τ Ω ρ(•, A, •)Q dx ≤ χ(A)e A/τ Ω ρ 2 (•, A, •) dx 1/2 Ω Q 2 dx 1/2 . (3.20)
Hence, applying again (3.12) and (3.13),

χ(A)e A/τ Ω ρ(•, A, •)Q dx ≤ χ(A)e A/τ Ω Q 2 (t -A, •) dx 1/2 Q 2 , if t > A, ≤ χ(A)e 2t/τ Ω ρ 2 0 (A -t, •)e 2(A-t)/τ dx 1/2 Q 2 , if t < A.
d Q 2 2 dt ≤ c 1 Q 2 + t 0 Ω Q 2 (t-a, •) dx 1/2 da+c 3 ρ 0 2 +c 4 H(t-A) Ω Q 2 (t-A, •) dx 1/2 Q 2 . (3.22) Setting F (t) = sup s∈[0,T ) Q(s) 2 2 + sup s∈[0,T ) ρ(s) 2 2 , (3.23)
we have

F ′ (t) ≤ max(0, d Q(t) 2 2 dt ) + max(0, d ρ(t) 2 2 dt ), (3.24)
and from (3.17) and (3.22) we deduce

F ′ (t) ≤ (c 1 + c 2 T + c 3 + c 4 ) F (t) F (t) + 1 τ 2 + 1 τ F (t) ≤ c 5 F (t), (3.25)
for a constant c 5 , from which we deduce that

F (t) ≤ e c5T F (0), (3.26)
and, as a consequence, that for a constant c 6 ,

Q(t) 2 2 + ρ(t) 2 2 ≤ c 6 ( Q 0 2 2 + ρ 0 2
2 ).

(3.27)

Finally we get (3.7) as a consequence of (3.27), this ends the proof.

As a consequence of the Lemma 3.2 we have the following Corollary whose biological meaning is: Biomass cannot be created from nothing. In order to establish the previous estimates, we have assumed that Q ≥ 0 and ρ ≥ 0. We can show that this is a consequence of the non negativity of Q 0 and ρ 0 . The proof of this Lemma is close to the one of Lemma 3.2.

Proof. We define ρ -= min(ρ, 0),

Q -= min(Q, 0). (3.29) Of course, ρ -(0, •, •) = 0, Q -(0, •) = 0, ρ -(•, 0, •) = ξ τ Q -. (3.30)
Now, multiplying (3.30) by e a/τ ρ -and integrating, since

∂ρ ∂t ρ -= ∂ρ - ∂t ρ -, (3.31)
proceeding as while establishing (3.17), we get

d ρ -2 2 dt ≤ 1 τ 2 Q -2 2 + 1 τ ρ -2 2 .
(3.32)

Multiplying now (1.1) by Q -and integrating gives

d Q -2 2 dt = 2 Ω 1 -ξ τ Q -2 dx + 2 A 0 Ω µρQ -e a/τ dxda + 2χ(A)e A/τ Ω ρ(•, A, •)Q -dx ≤ 2 Ω 1 -ξ τ Q -2 dx + 2 A 0 Ω µρ -Q -e a/τ dxda + 2χ(A)e A/τ Ω ρ -(•, A, •)Q -dx.
(3.33)

Since, as we had (3.12) and (3.13), we have here:

Ω ρ -2 dx = Ω ρ - 0 2 dx if t < a and Ω ρ -2 dx = Ω ξ 2 τ 2 Q -2 (t -a, •) dx if t > a, (3.34)
we can finish the proof as in the proof of Lemma 3.2 and get

Q -(t) + ρ -(t) ≤ c ′ Q - 0 + ρ - 0 = 0, (3.35)
giving the Lemma.

Estimates on the derivatives

As a by product of the proof of Lemma 3.2, we can deduce from (3.16) the following result insuring a first control on the regularity of ρ.

Corollary 3.5 If the assumptions and the conclusion of Theorem 1.1 are valid then, for any 0 ≤ s ≤ T ,

s 0 ∇ρ 2 2 dt = s 0 A 0 Ω |∇ρ| 2 2 dxdadt ≤ c d , (3.36)
for a constant c (depending only on A, τ , T , μ, ρ 0 2 and Q 0 2 ).

Proof. Integrating (3.16) from 0 to s yields:

χ(A)e A/τ s 0 Ω ρ 2 (•, A, •) dxdt + 2 s 0 A 0 Ω (D + d)|∇ρ| 2 e a/τ dxdadt + 2 s 0 A 0 Ω µρ 2 e a/τ dxdadt = s 0 Ω ξ 2 τ 2 Q 2 dxdt + 1 τ s 0 ρ 2 2 dt + ρ(0) 2 2 -ρ(s) 2 2 . (3.37)
Using the previous estimate concerning ρ and Q, we may deduce

d s 0 A 0 Ω
|∇ρ| 2 e a/τ dxdadt ≤ c, (3.38) and (3.36) follows.

Because of the form of the non linearity in (1.1)-(1.6), we need a supplementary estimate concerning

T 0 ∆ρ 2 2 dt = T 0 A 0 Ω |∆ρ| 2 2 dxdadt. (3.39)
This estimate is a consequence of an estimate on ∇ρ 4 and on

ρ ∞ = sup |ρ(t, a, x)|, t ∈ [0, T ), a ∈ [0, A), x ∈ Ω , Q ∞ = sup |Q(t, x)|, t ∈ [0, T ), x ∈ Ω , (3.40)
that we now set.

Lemma 3.6 The solution (Q, ρ) given by Theorem 1.1 satisfies

ρ ∞ + Q ∞ ≤ k, (3.41)
where k is a constant depending only on T , A, sup t∈[0,T ) ρ(t) 2 and sup t∈[0,T ) Q(t) 2 (which are estimated by Lemma 3.7)

Proof.

As we already see that ρ(t) 2 is bounded, using a method similar to Ladyzenskaja, Solonnikov and Ural'ceva [START_REF] Ladyzenskaja | Linear and Quasi-linear Equation of Parabolic Type[END_REF], (paragraph III-8) we deduce that

χ(A)ρ(•, A, •)e A/τ ≤ k 1 , (3.42)
where the constant k 1 only depends only on A and sup t∈[0,T ) ρ(t) 2 . Then defining

P = A 0 ρe a/τ da, (3.43)
we deduce from (1.2) that P is solution to the following parabolic equation

∂ P ∂t + χ(A)ρ(•, A, •)e A/τ = - A 0 µρe a/τ da + ∇ • (D + d)∇ P + ξ τ Q, (3.44)
from which we get that

P ∞ = sup | P (t, x)|, t ∈ [0, T ), x ∈ Ω ≤ k 2 , (3.45)
where k 2 only depends on sup Proof. Multiplying equation (1.2) by -∇ • |∇ρ| 2 ∇ρ e a/τ , and integrating in a and x gives:

t∈[0,T ) ρ(t) 2 , sup t∈[0,T ) Q(t)
- A 0 Ω ∂ρ ∂t ∇ • |∇ρ| 2 ∇ρ e a/τ dxda - A 0 Ω ∂ρ ∂a ∇ • |∇ρ| 2 ∇ρ e a/τ dxda = A 0 Ω µρ ∇ • |∇ρ| 2 ∇ρ e a/τ dxda - A 0 Ω ∇ • D + d ∇ρ ∇ • |∇ρ| 2 ∇ρ e a/τ dxda. (3.47)
Making a double integration by part, and following a straightforward computation procedure, we get

A 0 Ω ∇ • D + d ∇ρ ∇ • |∇ρ| 2 ∇ρ e a/τ dxda = A 0 Ω D + d |∇ρ| 2 |∇ 2 ρ| 2 + 2H(∇ρ, ∇ 2 ρ) e a/τ dxda + E, (3.48)
where

|∇ 2 ρ| 2 = ∂ 2 ρ ∂x 2 1 2 + 2 ∂ 2 ρ ∂x 1 ∂x 2 2 + ∂ 2 ρ ∂x 2 2 2 , (3.49) H(∇ρ, ∇ 2 ρ) = ∂ρ ∂x 1 ∂ 2 ρ ∂x 2 1 + ∂ρ ∂x 2 ∂ 2 ρ ∂x 1 ∂x 2 2 + ∂ρ ∂x 1 ∂ 2 ρ ∂x 1 ∂x 2 + ∂ρ ∂x 2 ∂ 2 ρ ∂x 2 2 2 , ( 3
|E| ≤ C 1 A 0 Ω |∇P | + |∇Q| |∇ρ| 3 |∇ 2 ρ| e a/τ dxda ≤ d 4 A 0 Ω |∇ρ| 2 |∇ 2 ρ| 2 e a/τ dxda + C 2 1 d A 0 Ω |∇P | + |∇Q| 2 |∇ρ| 2 e a/τ dxda; (3.53)
in order to get the last expression in (3.53), we used

U V ≤ d 4 U 2 + 1 d V 2 with V = C 1 |∇P | + |∇Q| |∇ρ|.
Concerning the other terms of (3.47)

d ∇ρ 4 4 dt + 1 4 χ(A) Ω |∇ρ(•, A, •)| 4 e A/τ dx + A 0 Ω D + d |∇ρ| 2 |∇ 2 ρ| 2 e a/τ dxda + 2 A 0 Ω D + d H(∇ρ, ∇ 2 ρ) e a/τ dx ≤ C 3 4τ 4 Ω |∇Q| 4 dx + 1 4τ ∇ρ 4 4 + C 2 ∇ρ 4 4 + 1 + d 4 A 0 Ω |∇ρ| 2 |∇ 2 ρ| 2 e a/τ dxda + C 2 1 d A 0 Ω |∇ρ| 4 e a/τ dxda + C 4 ∇Q 4 4 + ∇ρ 4 4 , (3.57)
and passing the fourth term of the right hand side in the left hand side we can deduce

d ∇ρ 4 4 dt ≤ C 5 d ∇ρ 4 4 + ∇Q 4 4 + 1 , (3.58)
for a constant C 5 .

Multiplying equation (1.2) by -∇ • |∇ρ| 2 ∇ρ e 4a/τ , and making the same operations as previously, we obtain an inequality which is (3.57) with e a/τ replaced by e 4a/τ and ∇Q 4 4 replaced by A 0 Ω |∇ρ| 4 e 4a/τ dxda. From this, we can deduce

d A 0 Ω |∇ρ| 4 e 4a/τ dxda dt + χ(A) Ω |∇ρ(•, A, •)| 4 e 4A/τ dx ≤ C 6 d A 0 Ω |∇ρ| 4 e 4a/τ dxda + ∇Q 4 4 + 1 . (3.59)
On another hand, computing the gradient of (1.1), and multiplying by ∇Q|∇Q| 2 yields 1 4

d ∇Q 4 4 dt = 1 τ Ω ∂ξ ∂Q Q|∇Q| 4 dx + Ω 1 -ξ τ |∇Q| 4 dx + Ω A 0 ∇µ ρ e a/τ da • ∇Q|∇Q| 2 dx + Ω A 0 µ∇ρ e a/τ da • ∇Q|∇Q| 2 dx + χ(A) Ω ∇ρ(•, A, •)e A/τ • ∇Q|∇Q| 2 dx. (3.60)
Because of the regularity of ξ and µ and of Lemma 3.6, since, applying Young's inequality,

Ω A 0 µ∇ρ e a/τ da • ∇Q|∇Q| 2 dx ≤ Ω A 0 µ|∇ρ| e a/τ da 4 dx 1/4 Ω |∇Q| 3 4/3 dx 3/4 ≤ μ A 0 Ω |∇ρ| 4 e 4a/τ dxda 1/4 ∇Q 3 4 , (3.61)
and

χ(A) Ω ∇ρ(•, A, •) e A/τ • ∇Q|∇Q| 2 dx ≤ χ(A) Ω |∇ρ(•, A, •)| 4 e 4A/τ dx 1/4 ∇Q 3 4 ≤ χ(A) C ′ 6 d A 0 Ω
|∇ρ| 4 e 4a/τ dxda 1/4

+ ∇Q 4 + 1 ∇Q 3 4 , (3.62)
we deduce from (3.60)

d ∇Q 4 4 dt ≤ C 4 ( ∇Q 4 4 + 1) + C 5 d A 0 Ω
|∇ρ| 4 e 4a/τ dxda 1/4

+ ∇Q 4 + 1 ∇Q 3 4 . (3.63) Inequalities (3.58), (3.59), (3.63)
and the assumptions on Q 0 and ρ 0 give

∇Q 4 + ∇ρ 4 + A 0 Ω |∇ρ| 4 e 4a/τ dxda 1/4 ≤ C d , (3.64)
and finally the Lemma.

Lemma 3.8 The solution (Q, ρ) given by Theorem 1.1 satisfies

∇Q(t) 2 2 ≤ c d , ∇ρ(t) 2 2 ≤ c d for any 0 ≤ t ≤ T, (3.65)
s 0 ∆ρ(t) 2 2 dt = s 0 A 0 Ω |∆ρ(t)| 2 e a/τ dxdadt ≤ c d 2 , for any 0 ≤ s ≤ T, (3.66)
s 0 ∂Q ∂t (t) 2 2 dt ≤ c, s 0 ∂ρ ∂t (t) 2 2 dt ≤ c d 2 and s 0 ∂ρ ∂a (t) 2 2 dt ≤ c d 2 for any 0 ≤ s ≤ T, (3.67)
for a constant c (which does not depend on d).

Proof. Multiplying (1.2) by (-∆ρ)e a/τ and integrating in a and x gives

d ∇ρ 2 2 dt + χ(A)e A/τ Ω ∇ρ(•, A, •) 2 dx + 2 A 0 Ω µ ∇ρ 2 e a/τ dxda + 2 A 0 Ω (D + d) ∆ρ 2 e a/τ dxda = A 0 Ω ∇µ • ∇ρ ρ e a/τ dxda - A 0 Ω ∇ D(M, Q, P ) • ∇ρ∆ρ e a/τ dxda + 1 τ ∇ρ 2 2 + Ω ξ 2 + ( ∂ξ ∂Q Q) 2 |∇Q| 2 dx ≤ c 1 ∇ρ 2 + ∇ρ 2 2 + ∇Q 2 2 + A 0 Ω ∂D ∂M ∇M + ∂D ∂P ∇P + ∂D ∂Q ∇Q ∇ρ∆ρ e a/τ dxda ≤ c 1 ∇ρ 2 + ∇ρ 2 2 + ∇Q 2 2 + 1 d A 0 Ω ∂D ∂M

Existence and uniqueness of the solution

Once the a priori estimates are set, the proof of existence is classical and in the spirit of Ladyzenskaja, Solonnikov and Ural'ceva [START_REF] Ladyzenskaja | Linear and Quasi-linear Equation of Parabolic Type[END_REF]. It essentially consists in linearization and passing to the limit using the estimates.

In the following L p (Ω) and L p ([0, A)× Ω, e a/τ dadx) are the functional spaces associated with the norms (3.1), W k,p (Ω) and W k,p ([0, A) × Ω, e a/τ dadx) are the Sobolev spaces composed of functions whose derivatives up to order k are in L p (Ω) or L p ([0, A) × Ω, e a/τ dadx). L p ([0, T ) × Ω) and L p ([0, T ) × [0, A) × Ω, e a/τ dadxdt) are the spaces of functions having finite norm Finally, for a functional space W , L ∞ (0, T ; W ) stands for the functions whose norm in W is finite for any t ∈ [0, T ].

Linearization

We linearize the system (1.1)-(1.6) Then using classical results on pde and ode, we give an existence and uniqueness result for the solution to this linearized system. We set Q 0 = Q 0 and ρ 0 = ρ 0 and for n ∈ N * , we consider (Q n , ρ n ) solution to:

∂Q n ∂t = 1 -ξ τ Q n + A 0 ρ n (•, a, •)e a/τ µ(•, a, •) da + χ(A)ρ n (•, A, •)e A/τ , on [0, T ) × Ω, (4.2)
∂ρ n ∂t + ∂ρ n ∂a = -µρ n + ∇ • D(M n-1 , Q n-1 , P n-1 ) + d ∇ρ n , on [0, T ) × [0, A) × Ω, (4.3)
ρ n (•, 0, •) = ξ τ Q n , on [0, T ) × Ω, (4.4
)

ρ n (0, •, •) = ρ 0 , on [0, A) × Ω, (4.5)
Q n (0, •) = Q 0 , on Ω, (4.6) ∂ρ n ∂ -→ ν = 0, on [0, T ) × [0, A) × ∂Ω, (4.7)
where for n ∈ N,

P n (t, x) = A amin ρ n (t, a, x)e a/τ da, (4.8)
and M n is solution to

∂M n ∂t = 1 P max -p max H r P n -p max P max -p max H r (1 -M n ) - 1 p min -P min H r p min -P n p min -P min H r (M n ), M n (0, •) = M 0 . (4.9)
Theorem 4.1 Under assumptions (1.8),(1.9),(1.10) and (1.13), if

ρ 0 ≥ 0 ∈ L 1 ∩W 2,2 ∩W 1,4 ([0, A)× Ω, e a/τ dadx) satisfies (1.15) and if Q 0 ≥ 0 ∈ L 1 ∩ W 2,2 ∩ W 1,4 (
Ω) then for any n ∈ N, there exists a unique solution

(Q n , ρ n ) ∈ C 0 b 0, T ; (L 1 ∩ W 1,2 (Ω)) × (L 1 ∩ W 1,2 ([0, A)
× Ω, e a/τ dadx)) to system (4.2) -(4.7) coupled with (4.8) and (4.9). Moreover, Q n ≥ 0, ρ n ≥ 0 and (Q n , ρ n ) satisfies estimates (3.2), (3.7), (3.41), (3.46), (3.65), (3.66) and (3.67) with constants independent of n.

The proof of this Theorem uses only classical pde and ode arguments. Hence we only sketch it. Proof. The proof consists in an induction procedure. Because of the assumptions on (Q 0 , ρ 0 ) and the definition of (Q 0 , ρ 0 ), the Theorem is true for n = 0. Then, if the Theorem is true for n -1, by regularization arguments, we can get that that M n-1 exists and is unique on [0, T ) × Ω and that M

n-1 ∈ C 0 b (0, T ; W 1,2 (Ω)) ∩ C 1 b (0, T ; L ∞ (Ω)
). Hence we deduce that, for each l ∈ N * , there exists a unique solution (

Q n,l , ρ n,l ) ∈ C 0 b 0, T ; (L 1 ∩ W 1,2 (Ω)) × (L 1 ∩ W 1,2 ([0, A) × Ω, e a/τ dadx)) to ∂Q n,l ∂t = 1 τ Q n,l - ξ τ Q n,l-1 + A 0 ρ n,l (•, a, •)e a/τ µ(•, a, •) da + χ(A)ρ n,l (•, A, •)e A/τ , (4.10
)

∂ρ n,l ∂t + ∂ρ n ∂a = -µρ n,l + ∇ • D(M n-1 , Q n-1 , P n-1 ) + d ∇ρ n,l , (4.11)
ρ n,l (•, 0, •) = ξ τ Q n,l-1 , ρ n,l (0, •, •) = ρ 0 , ∂ρ n,l ∂ -→ ν ∂Ω = 0, (4.12) Q n,l (0, •) = Q 0 , (4.13)
where Q n,0 is defined as Q n,0 = Q n-1 . This deduction involves first a classical semi-group or Galerkin routine in order to deduce that there exists a unique solution ρ n,l to (4.11) -(4.12) as soon as the existence of (Q n,l-1 , ρ n,l-1 ) is achieved. These routines are explained in Lions and Magenes [START_REF] Lions | Problèmes aux limites non homogènes et applications[END_REF], Ladyzenskaja, Solonnikov and Ural'ceva [START_REF] Ladyzenskaja | Linear and Quasi-linear Equation of Parabolic Type[END_REF], or -in a context close to our-in Langlais [START_REF] Langlais | A nonlinear problem in age-dependent population diffusion[END_REF].

Once the existence of ρ n,l is established, the existence and uniqueness of Q n,l follows. Now, following the way leading to (3.17) and (3.22) we deduce that (Q n,l , ρ n,l ) satisfies

d ρ n,l 2 2 dt ≤ 1 τ 2 Q n,l-1 2 2 + 1 τ ρ n,l 2 2 , (4.14) d Q n,l 2 2 dt ≤ c 1 Q n,l 2 + c ′ 1 Q n,l-1 2 t 0 Ω (Q n,l ) 2 (t -a, •) dx 1/2 da + c 3 ρ 0 2 + c 4 H(t -A) Ω (Q n,l ) 2 (t -A, •) dx 1/2 Q n,l 2 , (4.15)
which is enough to deduce that ( ρ n,l 2 + Q n,l 2 ) is bounded. As a consequence of this bound, we get that, for a subsequence still denoted l,

(Q n,l , ρ n,l ) -→ (Q n , ρ n ) in L ∞ 0, T ; (L 2 (Ω))×(L 2 ([0, A)×
Ω, e a/τ dadx)) weakly- * , where (Q n , ρ n ) is solution to (4.2) -(4.7). Finally, the estimates are led in the same way as in section 3. The uniqueness follows directly from the linear character of (4.2) -(4.7). Hence, the Theorem is true for n.

The induction procedure is then straightforward to end the proof of the Theorem.

Existence

From estimates (3.65), (3.66) and (3.67) we deduce that the sequence (

Q n , ρ n ) is bounded in W 1,2 ([0, T ) × Ω) × W 1,2 ([0, T ) × [0, A) × Ω).
Hence, up to a subsequence still denoted n, we have

(Q n , ρ n ) -→ (Q, ρ) in W 1,2 ([0, T ) × Ω) × W 1,2 ([0, T ) × [0, A)
× Ω, e a/τ dadxdt) weakly, and then, in L 2 ([0, T ) × Ω) × L 2 ([0, T ) × [0, A) × Ω, e a/τ dadxdt) strongly. From this we can also deduce that P n -→ P in L 2 ([0, T ) × Ω) strongly, with P defined from ρ by (1.11). In view of (4.9), we can deduce that (M n ), (∂M n /∂t) and, taking the gradient of (4.9), (∇M n ) are bounded in L 2 ([0, T ) × Ω). Then extracting again a subsequence still denoted n, we deduce M n -→ M strongly, where M is solution to (1.12). Using now the regularity of D, we obtain D(M n-1 , Q n-1 , P n-1 ) -→ D(M, Q, P ) in L 2 ([0, T ) × Ω) strongly.

The regularity of trace operators gives χ(A)ρ n (•, A, •) -→ χ(A)ρ(•, A, •), ρ n (•, 0, •) -→ ρ(•, 0, •), ρ n (0, •, •) -→ ρ(0, •, •), Q n (0, •) -→ Q n (0, •) weakly, and using (3.66), ∂ρ n /∂ -→ ν |∂Ω -→ ∂ρ/∂ -→ ν |∂Ω weakly.

Then passing to the limit in (4.2) -(4.9), we obtain that (Q, ρ) is solution to (1.1) -(1.6) coupled with (1.11) and (1.12). Once this existence result is established, using regularizations and truncations, we can start the computations of section 3 giving the additional regularity and the non negativity of Q and ρ.

It now remains to prove the uniqueness of the solution.

Uniqueness

Consider (Q, ρ) with associated P and M and ( Q, ρ) with associated P and M two solutions of (1.1) -(1.6). They both satisfy the estimates and the difference ( Q, ρ) = (Q -Q, ρ -ρ) satisfies for a constant k 4 From the last three inequalities we deduce, for a constant K

∂ Q ∂t = 1 -ξ τ Q + A 0 ρ(•,
≤ k 1 ρ 2 2 + Q 2 2 + k 2 ∇ρ ∞ ρ 2 + Q 2 ∇ρ 2 ≤ k 1 ρ 2 2 + Q 2 2 + k 2 2 d ∇ρ 2 ∞ ρ 2 + Q 2 2 + d 4 ∇ρ 2
Q 2 2 + ρ 2 2 + A 0 Ω
|ρ| 2 e 2a/τ dxda ≤ K Q|t=0

. 8 )

 8 Now, we explain in what sense the model (1.1)-(1.6

ρ 2

 2 e a/τ da + ρ 2 e a/τ a=A a=0 ,

(3. 21 )

 21 Using(3.19),(3.21) and(1.15) in(3.18), for four non negative constant c 1 , c 2 , c 3 and c 4 we have

Corollary 3 . 3

 33 If the assumptions of Theorem 1.1 are realized and if ρ 0 = 0 and Q 0 = 0 then the solution (Q, ρ) given by the Theorem satisfies ρ = 0 and Q = 0. (3.28)

Lemma 3 . 4

 34 If the assumptions of Theorem 1.1 are true and if there exists a solution (Q, ρ) to system (1.1)-(1.6), then Q ≥ 0 and ρ ≥ 0.

Lemma 3 . 7 ∇Q 4 + ∇ρ 4

 3744 a/τ da| which is finite by assumption. Then, (3.45) and(3.42) give that Q ∞ is finite, and as a consequence, (1.2) and (1.3) finally give the bound on ρ ∞ , ending the proof. The solution (Q, ρ) given by Theorem 1.1 satisfies C (C does not depend on d).

  W k,p ([0, T ) × Ω) and W k,p ([0, T ) × [0, A) × Ω, e a/τdadxdt) are their associated Sobolev spaces.

  a, •)e a/τ µ(•, a, •) da + χ(A)ρ(•, A, •)e A/τ , + ∇ • D(M, Q, P ) + d ∇ρ -∇ • D( M, Q, P ) -D(M, Q, P ) ∇ρ , (4M, Q, P ) + d)|∇ρ| 2 e a/τ dxda , Q, P ) -D(M, Q, P ) ∇ρ • ∇ρ e a/τ dxda

  |t=0 | 2 e 2a/τ dxda = 0, (4.23) giving Q = ρ = 0 and then the uniqueness of the solution to (1.1) -(1.6).

  ∇µ ρ |∇ρ| 2 e a/τ dxda ≤ C 2 ∇ρ 4 4 + 1 , (3.56) using the regularity of µ and the estimate on sup(ρ) given by Lemma 3.6. The regularity of ξ and the estimate on sup(Q) give |ξ 4 + ( ∂ξ ∂Q Q) 4 | ≤ C 3 for a constant C 3 .

	In a similar way,
	-	0	A	Ω	∂ρ ∂a	∇ • |∇ρ| 2 ∇ρ e a/τ dxda =	1 4	0	A	Ω	∂(|∇ρ| 4 e a/τ ) ∂a	dxda -	1 4τ	∇ρ 4 4
				=	1 4	χ(A)	Ω	|∇ρ(•, A, •)| 4 e A/τ dx -	1 4τ 4 Ω	ξ 4 + (	∂ξ ∂Q	Q) 4 |∇Q| 4 dx -	1 4τ	∇ρ 4 4 , (3.55)
	and
		A		µρ ∇ • |∇ρ| 2 ∇ρ e a/τ dxda
	0		Ω
				= -	A		µ|∇ρ| 4 e a/τ dxda -	A
							0		Ω					0	Ω
															Hence
	(3.47) yields
	1
	4
	integration by parts, we get				, since	∂|∇ρ| 4 ∂t	= 2	∂|∇ρ| 2 ∂t	|∇ρ| 2 = 4	∂∇ρ ∂t	• ∇ρ |∇ρ| 2 , making an
										-	0	A	Ω	∂ρ ∂t	∇ • |∇ρ| 2 ∇ρ e a/τ dxda =	1 4	d ∇ρ 4 4 dt	.	(3.54)

  . Now transferring the last term of the right hand side in the left hand side and using the estimate on ∇ρ 4 and ∇Q 4 which also give an estimate on ∇M 4 and ∇P 4 , we get, for a constant c 2 We could also prove that ∇Q ∞ and ∇ρ ∞ are bounded.

												∇M +	∂D ∂P	∇P +	∂D ∂Q	∇Q	2	∇ρ	2 e a/τ dxda
												+	d 4	0	A	Ω	∆ρ	2 e a/τ dxda, (3.68)
	for a constant c 1 d ∇ρ 2 2 dt	+	5d 4	∆ρ 2 2 ≤ c 2 ∇ρ 2 2 + ∇Q 2 2 +	1 d	,	(3.69)
	In a similar way, we can also get
	d(	A 0	Ω |∇ρ| 2 e 2a/τ dxda) dt	≤ c 3		0	A	Ω	|∇ρ| 2 e 2a/τ dxda + ∇Q 2 2 +	1 d	,	(3.70)
	and
			d ∇Q 2 2 dt	≤ c 4	0	A	Ω	|∇ρ| 2 e 2a/τ dxda + ∇Q 2 2 .	(3.71)
	From the three last inequalities we get (3.65). Integrating (3.69) from 0 to s gives (3.66). The
	estimate on ∇Q(t) 2 is then obtained as the estimate on ∇Q(t) 4 . Estimate (3.67) is finally a
	direct consequence of equations (1.1) and (1.2).
	Remark 3.1

  2 , (4.19) for constants k 1 and k 2 . Passing the term d 4 ∇ρ 2 2 in the left hand side gives, for a constant k 3

						d ρ 2 2 dt	≤	k 3 d	ρ 2 2 + Q 2 2 .	(4.20)
	Making the same, but multiplying (4.17) by ρe 2a/τ yields
	d	A 0	Ω |ρ| 2 e 2a/τ dxda dt	+ χ(A)e 2A/τ	Ω	ρ(•, A, •) dx ≤	k 3 d	0	A	Ω	|ρ| 2 e 2a/τ dxda + Q 2 2 . (4.21)
	Multiplying (4.16) by Q gives
				d Q 2 2 dt	≤	k 4 d		0	A	Ω	|ρ| 2 e 2a/τ dxda + Q 2 2 ,	(4.22)

* This work has been partially supported by "La Région Bretagne", F-35031 Rennes, Program 1042: "Renouvellement


Katsunori Saito
email: k-saito@ac.daido.co.jp

Teppei Okumura

Kazumi Yamamoto

Large Scale Cold Crucible Levitation Melting Furnace for Titanium

Keywords:

We, Daido Steel Co., Ltd., began titanium alloy production at our first Levitation Induction melting Furnace, we named LIF, in our Hoshizaki Plant in 2007. The LIF is a cold crucible semi-levitation melting furnace (CCLM) that has bottom tapping nozzle. A 1.8 metric tons ingot has been made by a combination of LIF and Vacuum Arc Remelting. The LIF performs completely uniform melting and enables to melt high-melting point metals.Therfore, we can produce high grade titanium alloys with high melting point metals such as pure tantalum or pure niobium. We has produced many kinds of titanium alloys for 7 years by LIF. However there are some issue in melting and bottom tapping process. In this article, we show those problems.

Introduction

Titanium alloys are being widely used nowadays as they have advantages of having low density, high corrosion resistance, high specific strength, and biocompatibility. They are used biological implant, automobile, aircraft, glasses, golf clubs and chemical plants. They are expected as suitable materials for engine parts of automobile such as connecting rod, and exhaust valves employing their high specific strength and good performance at elevated temperature. [START_REF] Shibata | Proceedings of 2 nd Electromagnetic Processing of Materials[END_REF] It has been achieved to use them for parts of the aircraft in recent years. Since titanium alloys are highly reactive metals, they react readily with oxygen during melting and casting. Their ductility and toughness decrease when contents of oxygen increase too much. Therefore, a cold crucible is used for titanium melting process instead of a refractory crucible. CCLM furnace is a type of induction melting furnace that is replaced a conventional refractory crucible with a watercooled copper crucible. High frequency coil current induces magnetic flux density and eddy current on the surface of the metal. Eddy current melt the metal by Joule heat. On CCLM, Lorentz force generated by the eddy current and the magnetic flux density confine the metal and keep noncontact form the water-cooled copper crucible. That is well known as a titanium melting furnace with melting has several advantages, such as no restriction for raw material's shape and melting atmosphere compare to vacuum arc melting (VAR), electron beam melting (EB), and plasma arc melting (PAM). Tilting tap is conventional method to tap molten metal and use well in small scale of cold crucible. However, for the cold crucible, temperature of molten metal rapidly drops when it contact to a crucible during tapping, and it caused decrease of ingot yield. In addition, when large scale of cold crucible is designed, it is difficult to design water-cooled and segmented crucible and co-axial conductor that is suitable for tilting tap. So the bottom tapping method was developed. Bottom tapping has many advantages for cold crucible because no temperature drop, do not require watercooled cable nor co-axial conductor which has difficulty for high frequency furnace, do not require tundish and so on. [START_REF] Kusamichi | [END_REF]3] Therefore, a new process of manufacturing several hundred kilograms ingots of titanium alloys expected to construct with utilizing large scale cold crucible with bottom tapping nozzle. The shape of the crucible and bottom-tapping nozzle was designed properly through the simulations and some experiments. Titanium alloy melting and bottom tapping were successfully performed with this process.

Features of LIF

There are two processes in the CCLM. Full levitation melting is a process to float the molten metal completely, and semi levitation melting is a process to maintain the molten metal on the solidified layer of the metal (skull). Full levitation melting can melt metals of high melting point and allows metals to be melted without contamination, but the power-consumption is very large. So, semi levitation melting was selected for the development of LIF. The skull formed at bottom of the crucible itself supports molten metal since the Lorentz force is not strong enough to levitate whole metal, including lower part of crucible. Still the Lorentz force supports metal confinement without touching the sidewall, the metal temperature is kept high. Fig. 1 shows LIF. The specification of LIF is listed in table 1. The watercooled copper crucible is composed of 60 segments. The water-cooled bottom tapping nozzle, which is composed of copper segments, is located at the bottom of the crucible Coils are wound around the crucible and the nozzle. These are connected to different high-frequency power supplies. Electromagnetic levitation force is described with the equation ( 1) when eddy current flows near surface of a conductor. [4,[START_REF] Tadokoro | Proceedings of the 3 th International Symposium on Electromagnetic Processing of Materials[END_REF] It is derived from equation ( 1) that the lower the frequency is, the stronger the electromagnetic levitation force is if consumed power is constant. On the other hand, induced eddy current concentrated at the surface of the conductor leads high heating effect when frequency becomes higher. The upper coil for LIF induces strong electromagnetic force for levitating and stirring molten metal and the lower coil for bottom tapping nozzle supplies heat to metal in nozzle rather than confinement force. [START_REF] Okamura | Processing of 5 th International Symposium on Electromagnetic Processing of Materials[END_REF] Melting target LIF has melted titanium alloy since 2007. Table-2 shows amount of commercial melting of titanium alloy by additional alloy. In most cases we use popular alloy (Al-V and Al-Nb etc.) in titanium alloy melting. But in rare cases we use high melting-point alloy (pure molybdenum and pure niobium etc.). They have not occurred high density inclusion.

   f W F 2  ( 1

Challenges in commercial production

We have some challenges in commercial production.

Difficulty bottom tapping

We cannot be tapped molten metal from tapping nozzle in Specially titanium alloy (Ti-20%Nb-5%Cr-4%Zr) first tapping.

Turbulence flow in bottom tapping

We cannot cast ingot in a mold when turbulence occurs in bottom tapping. Photo-1 shows turbulence occurs in bottom tapping. Photo-2 shows top of mold after casting.

Difficulty bottom tapping

In the case of Ti-20%Nb-5%Cr-4%Zr (TNCZ), we cannot be tapped from bottom tapping nozzle. We think material characteristics influence it. So we examine by using material characteristics. We cannot be tapped from bottom tapping nozzle in case of TNCZ, because the skull is thick. The reason is molten metal temperature is low. So we think that Power ratio applied to the molten metal in case of TNCZ is lower than Ti-4%Al-22%V on immediately before bottom tapping. So we compare TNCZ with Ti-4%Al-22%V on melting. Fig. 2 shows power outlet and power efficiency to the metal in TNCZ. Fig. 3 shows power outlet and power efficiency to the metal in Ti-4%Al-22%V. Power distribution to the metal drop to a lower level as melting goes on and it becomes 0% on immediately before bottom tapping. But, Result of comparison between TNCZ and Ti-4%Al-22%V, there is no difference. Fig. 4 shows total power efficiency on melting. It means crucible loss energy is larger in TNCZ than Ti-4%Al-22%V. So it is shown that area in contact with the crucible on TNCZ is larger than Ti-4%Al-22%V. Therefore we derived a contact area with the crucible from protuberance height. Fig. 5 shows molten metal shape in LIF. Gravity acting on the molten metal=Force to lift the molten

S f W gh     2  (1)
γ:specific gravity, g:gravitational acceleration, h: protuberance height W: Power outlet, η:Power distribution ratio to the metal μ:Magnetic permeability, ρ:Specific resistance, f:Power supply frequency, S:Surface area on protuberance, d:diameter of crucible

We put the case that protuberance shape is spheroid.

h d dh h d d h d f W gh 2 4 tanh 4 4 2 4 2 2 2 2 1 2 2 2 2 2 2             (4)
Obtained the answer by that it will assign parameter in case of table-3 to (4).

Conclusion

LIF can melt high melting point metal (pure niobium, pure tantalum etc.). There are many technical issues in bottom tapping, so we need to any research. We think protuberance height and power distribution to the metal have a great influence.

  ) F : electromagnetic levitation force [N] W : consumed power [W] μ : permeability [A/m] f : frequency [Hz] ρ : specific resistivity [Ωm]

Fig. 1 :

 1 Fig.1: Schematic diagram of LIF with bottom tapping nozzle

  Fig.2 TNCZ Melting Fig.3 Ti-4%Al-22%V Melting

Fig. 5

 5 Fig.4 Total power efficiency on melting

  d

Table - 2

 - Amount of commercial melting of titanium alloy

	Additional master alloy	Melting weight(MT)
	Aluminum and Vanadium	2250
	Aluminum and Niobium	348
	Pure Molybdenum	5
	Pure Niobium	2
	Pure Tantalum	2

Table - 1

 - Specifications of the LIF with bottom tapping nozzle

	Melting power	2400kW	Melting crucible	Φ670mm×H700mm	Casting mold	Φ440mm×H3600mm
	Frequency	800Hz	Number of segment	60	Melting chamber Φ2600mm×H2000mm
	Melting weight	500kg-Ti	Cooling water flux	192m 3 /h	Casting chamber Φ2000mm×H6000mm

Table - 3

 - Calculation conditionThe calculated results were obtained answers in Table-4. It shows that when protuberance height is small, the contact area grows. So we think molten metal temperature is lowered by heat removal amount increases, it estimated that the skull thickness is increased.

			Ti-4%Al			Ti-4%Al
	Listing	TNCZ		Listing	TNCZ
			-22%V			-22%V
	specific gravity γ(kg/m 3 )	5200	4690	Magnetic permeability μ(H/m)	1.0	1.0
	gravitational acceleration			Specific resistance
	g(kg・m/s 2 )	9.8	9.8	(Ωm) ρ×10 -6	1.29	1.48
				Power supply
	Power outlet W(kW)	2250	2250		820	820
				frequency f(Hz)
	Power distribution ratio			diameter of crucible
	to the metal η(%)	10.8	15.3	(mm)	600	600
	Table-4 Calculation result
				Ti-4%Al
	Listing		TNCZ
				-22%V
	protuberance height (mm)	168	220
	Contact area with crucible ratio (%)	47.4	44.3


Baokuan Li

Qiang Wang

Kang Deng

A Multiscale Three-Dimensional Numerical Simulation of Electroslag Remelting Process

Keywords: electroslag remelting, electromagnetism, two-phase flow, solidification, macrosegregation, microstructure

A three-dimensional, transient and multiscale mathematical model has been presented, which is able to reveal the coupled physical fields in the electroslag remelting (ESR) process. The mold is assumed to be conductive rather than insulated. The electric current, Lorentz force and Joule heating fields were demonstrated by Maxwell's equations. The volume of fluid (VOF) approach was implemented for the metal droplet behavior. The solidification was modeled by an enthalpy-porosity formulation, in which the mushy zone was treated as a porous medium with porosity equal to the liquid fraction. Besides, the solute distribution was revealed by the continuum mixture model. The macro-model was linked to a meso-model through the flow and temperature fields. A regular network of square cells with a much finer scale was drawn for the grain structure with the cellular automaton (CA) technique. The continuous nucleation, which is based on the Gaussian distribution, was implemented to describe the heterogeneous nucleation. The growth kinetics of the dendritic tip was taken into account by the Kure-Giovanola-Trivedi (KGT) model. Moreover, a moving mesh was carried out to account the growing of ingot. The electric current flows to the mold lateral wall especially in the slag layer. A large amount of Joule heating around the metal droplet varies as it falls. The hottest region appears under the outer radius of the electrode tip, close to the slag/metal interface instead of the electrode tip. The metal pool becomes deeper with more power. The maximal temperature increases from 1951 K to 2015 K, and the maximum metal pool depth increases from 34.0 mm to 59.5 mm with the applied current ranging from 1000 A to 2000 A. The vertical columnar grains appear at the bottom of the ingot and an inverse V-shaped grain structure is observed at the upper part.

Introduction

ESR is an advanced technology of ingot production which is used in critical applications such as aeronautics, power generation, medicine, and tooling [START_REF] Hernandez-Morales | Review of mathematical models of fluid flow, heat transfer, and mass transfer in electroslag remelting process[END_REF]. The passage of an alternating current (AC) from the electrode to the baseplate creates a Joule heating in the highly resistive calcium fluoride-based slag, which is sufficient to melt the electrode. The interaction between the self-induced magnetic field and the AC gives rise to the Lorentz force. Metal droplets sink through the less dense slag to form a liquid metal pool in the water-cooled mold. Heat loss to the mold solidifies the metal, forming an ingot and maintaining a liquid metal pool throughout the process. The redistribution of the compositions in the metal pool creates the macrosegregation in the final ingot [START_REF] Ludwig | Modeling of multiscale and multiphase phenomena in materials processing[END_REF].

Given the complicated phenomena involved, and the difficulty and expense of performing experiments on a real apparatus, numerical simulations constitute an attractive approach. Dilawari et al. [START_REF] Dilawari | Heat Transfer and fluid flow phenomena in electroslag refining[END_REF] and Jardy et al. [START_REF] Jardy | Magnetohydrodynamic and thermal behavior of electroslag remelting slags[END_REF] developed a two-dimensional (2D) axisymmetric model, which accounts for the steady state Lorentz force, flow field and temperature distribution. In recent years, a transient comprehensive 2D axisymmetric model has been established by Weber et al. [START_REF] Weber | A comprehensive model of the electroslag remelting process: description and validation[END_REF]. Moreover, the transports of the solutes as well as the inclusions during the process were illustrated by Kelkar et al. [START_REF] Kelkar | Pro. of the 2013 International Symposium on Liquid Metal Processing & Casting[END_REF]. Yanke et al. [START_REF] Yanke | Pro. of the 2013 International Symposium on Liquid Metal Processing & Casting[END_REF] focused on the formation of the solidified slag skin. Kharicha et al. [START_REF] Kharicha | Pro. of the 2013 International Symposium on Liquid Metal Processing & Casting[END_REF] and Hugo et al. [START_REF] Hugo | Pro. of the 2013 International Symposium on Liquid Metal Processing & Casting[END_REF] investigated the impact of the mould current on the Joule heating distribution. Fezi et al. [START_REF] Fezi | Pro. of the 2013 International Symposium on Liquid Metal Processing & Casting[END_REF] represented the effect of the current on the segregation level in the ESR ingot. Nastac et al. [START_REF] Nastac | A multiscale transient modeling approach for predicting the solidification structure in VAR-Processed alloy 718 ingots[END_REF] adopted a stochastic numerical approach to model the formation of grain structure and secondary phases during the solidification of ESR ingot.

The present work aims to develop a transient 3D macro-model, which is able to show the varying of the physical fields in the ESR process. A one-way coupled meso-model is used to simulate grain structure development. A mesh model is built based on a real scale furnace, which was used to perform some experiments. The comparison of the experimental observations and simulated results is carried out.

Model Description Electromagnetism

The magnetic Reynolds number, which expresses the ratio of the magnetic convection to magnetic diffusion, remains very low in the ESR process. The magnetic diffusion therefore dominates the electromagnetic phenomenon. The quasistatic magnetic equations are solved, because the changing of the magnetic diffusion is fast enough. Moreover, we use the phasor notation to account for the periodic time dependence.

The magnetic field intensity continuous at the inlet and bottom, and is related to the current on the top surface of slag and the mold lateral wall. Finally, we can obtain the time-averaged Joule heating and Lorentz force fields.

Fluid Flow

The flow is modeled by the continuity equation and time-averaged Navier-Stokes equations. The source term t F  represents the thermal buoyancy determined by the Boussinesq approximation, and

s F 
is the solutal buoyancy induced by the solute redistribution. We treat the mushy zone as a porous medium where the porosity gradually decreases from 1 to 0 as the metal solidified, and thus

d F 
is the turbulent damping force solved by the Darcy's law. The permeability of the mushy zone, which is determined by the Kozeny-Carman equation, is related to secondary dendrite arm spacing. The movement of the slag and the metal is weakly turbulent. The RNG k-ε turbulence model, which is able to capture the behavior of flows with lower Reynolds number, therefore is employed to calculate the turbulent viscosity. An enhanced wall function is employed to work with the RNG k-ε turbulence model, particularly since liquid metal has a low Prandtl number. The two-phase flow is tracked with the VOF approach. The continuum surface force model is implemented to consider the surface tension between the metal and slag.

An iterative procedure is used to include the varying melt rate related to the Joule heating, and the changing melt rate is imposed on the inlet. A no-slip condition is applied to the lateral wall and bottom. A zero shear stress is adopted on the top surface of the slag.

Heat and Mass Transfer

The energy conservation equation of the enthalpy-formulation is employed to demonstrate the temperature distribution. The enthalpy of the metal and the slag is computed as the sum of the sensible enthalpy and the latent heat released in the mushy zone.

Local solute mass fraction is described by the species conservation equation. The convection diffusion is included, while the molecular diffusion is neglected. Besides, the elements motions are irrelevant. The liquidus and solidus temperatures are a function of the local solute concentration. The liquid fraction is updated by the liquid/solid interface temperature. The lever rule is used to describe the liquid/solid interface temperature, and the interfacial transfer of species at the liquid/solid interface as well as the species movement in the solid phase.

In order to consider the melting process in a more simplified way, the temperature of the metal at the inlet is given by a parabolic profile. This parabolic profile has an approximate 30 K superheat and a peripheral boundary temperature close to the metal liquidus temperature. Equivalent heat transfer coefficients are applied to the lateral wall and bottom. The concentration of each species at the inlet is equal to the nominal concentration in the AISI 201 stainless steel. CA Model A probabilistic approach based on the instantaneous nucleation and Gaussian distribution of nucleation sites is introduced in the present work. The grain growth is described by the KGT model depending on the local undercooling. The total undercooling is generally the sum of four contributions, which are associated with solute diffusion, thermal diffusion, attachment kinetics and solid-liquid interface curvature. Only the solute diffusion undercooling is considered, because the last three terms are negligibly small. The SDAS is a function of the temperature gradient, solidification rate and temperature difference.

Results and Discussions

Electromagnetic Fields Fig. 1 illustrates the simulated current density and Joule heating density fields at 1400 s with the current of 2000 A. The electric current flows downward from the inlet to the two sides. The maximal current density is found just near the periphery of the inlet and the minimum is located at the outer side of the top slag. The skin effect is not observed, because the depth of skin effect is approximately 85 mm larger than the ingot radius. Most Joule heating is generated in the slag due to the smaller electrical conductivity. Moreover, the distribution of the Joule effect is similar to that of the current density. Fig. 2 displays the calculated Lorentz force and phase distribution. The inward Lorentz force generates a pinch effect. Thus, the metal moves toward the middle and forms a droplet. According to the Biot-Savart law, the force at the outer edge is larger than that at the middle. MHD Thermosolutal Convection Fig. 3 represents the simulated temperature distribution. The highest temperature region is just at the same place with the maximal Joule effect mentioned above. The slag becomes colder from the top and a lower temperature region is observed at the outer edge of the bottom slag. With more heat, the slag is hotter than the metal. The colder metal droplet therefore is heated by the slag during its falling. Due to the forced cooling in the ESR process, a large thermal gradient is achieved leading to a faster freezing speed.

Two pairs of vortices are found in the slag layer as shown in Fig. 4. The heat extracted by the cooling water results in a descent of the slag at the wall due to buoyancy. This causes a stable clockwise circulation, which can be seen on the right side of the figure. The downward flow at the lateral wall turns the corner to form a radial jet along the slag/metal interface. The strength of the jet diminishes as it moves inward and finally turns upward. At the center of the slag, a counterclockwise circulation is caused by the Lorentz force and the falling metal droplets. The heat extracted by the cooling water at the bottom and lateral wall, which drives the freezing of metal. A shallow metal pool therefore is formed. The flow in the metal pool is different from that in the slag. Due to the forced cooling, the thermal buoyancy dominates the flow pattern. The clockwise cell in the pool results in the motion of heavy cold metal downward along the curved solidification front. Ni from the bottom to the upper part of the metal pool, giving rise to a wider zone with the ingot growing. In addition, a lower concentration region at the both sides of the bottom is observed, where the metal freezes faster at early times and the Ni is depleted. Grain Structure Fig. 6 shows the comparison of the grain morphology between the experiment and simulation. The calculated growth direction and size of the grain, as well as the metal pool profile agree closely with the experimental observation. The cooling at the bottom and lateral wall considerably influences the grain growth. The nucleation simultaneously occurs at the chill bottom and lateral wall yielding a thin layer of the equiaxed grain, and then the equiaxed grain transforms to the columnar grain. The growth direction of the columnar grain is approximately 45° with respect to the vertical axis referred to an inverse V-shaped grain structure. Moreover, the columnar grain can extend to the ingot center because of the large undercooling.

Fig. 1 Fig. 3

 13 Fig. 1 Current density and Joule heating fields. Fig. 2 Phase distribution and Lorentz force field.

Fig. 5

 5 Fig. 5 Ni concentration distributions.

Fig. 6

 6 Fig. 6 Grain structures of experiment and simulation. Conclusions The present work represents a transient 3D macro-model of the ESR process. The solution of the mass, momentum, energy, and species conservation equations are simultaneously implemented by the finite volume method with full coupling of the Joule heating and Lorentz force through solving the Maxwell's equations. Besides, a transient asymmetric meso-model is developed using the CA technique with a steady state Lorentz force and Joule heating fields. The electric field, Lorentz force field, Joule heating distribution, two-phase flow, temperature distribution, solidification sequence and the macrosegregation are clearly clarified. Moreover, two nondimensional criterions are first proposed to distinguish the effects of the Lorentz force, the thermal buoyancy and the solutal buoyancy. The evolution of the grain morphology with time is demonstrated. The present model is able to give a deeper understanding of the ESR process. Acknowledgements The authors' gratitude goes to National Natural Science Foundation of China (No. 51210007) and Fundamental Research Funds of People's Republic of China for the Central Universities (N140206002). References [1] B. Hernandez-Morales, A. Mitchell (1999), Review of mathematical models of fluid flow, heat transfer, and mass transfer in electroslag remelting process, Ironmak. Steelmak. 26, 423-438. [2] A. Ludwig, A. Kharicha, M. Wu (2014), Modeling of multiscale and multiphase phenomena in materials processing, Metall. Mater. Trans. B 45B, 36-43. [3] A.H. Dilawari, J. Szekely (1978), Heat Transfer and fluid flow phenomena in electroslag refining, Metall. Mater. Trans. B 9B, 77-87. [4] A. Jardy, D. Ablitzer, J.F. Wadier (1991), Magnetohydrodynamic and thermal behavior of electroslag remelting slags, Metall. Mater. Trans. B 22B, 111-120. [5] V. Weber, A. Jardy, B. Dussoubs, D. Ablitzer, S. Rybéron, V. Schmitt, S. Hans, H. Poisson (2009), A comprehensive model of the electroslag remelting process: description and validation, Metall. Mater. Trans. B 42B, 271-280. [6] K.M. Kelkar, S.V. Patankar, S.K. Srivatsa, R.S. Minisandram, D.G. Evans, J.J. deBarbadillo, R.H. Smith, R.C. Helmink, A. Mitchell, H.A. Sizek (2013), Pro. of the 2013 International Symposium on Liquid Metal Processing & Casting, Austin (USA), 3-12. [7] J. Yanke, M.J.M. Krane (2013), Pro. of the 2013 International Symposium on Liquid Metal Processing & Casting, Austin (USA), 71-78. [8] A. Kharicha, E.K. Sibaki, M. Wu, A. Ludwig (2013), Pro. of the 2013 International Symposium on Liquid Metal Processing & Casting, Austin (USA), 95-99. [9] M. Hugo, B. Dussoubs, A. Jardy, J. Escaffre, H. Poisson (2013), Pro. of the 2013 International Symposium on Liquid Metal Processing & Casting, Austin (USA), 79-85. [10] K. Fezi, J. Yanke, M.J.M. Krane (2013), Pro. of the 2013 International Symposium on Liquid Metal Processing & Casting, Austin (USA), 151-158. [11] L. Nastac (2014), A multiscale transient modeling approach for predicting the solidification structure in VAR-Processed alloy 718 ingots, Metall. Mater. Trans. B 45B, 44-50.

Acknowledgements

The authors' gratitude goes to National Natural Science Foundation of China (No. 51210007) and Fundamental Research Funds of People's Republic of China for the Central Universities (N140206002).