Lilian Sofia

Sepulveda Salcedo
email: lssepulveda@uao.edu.co

Michel De Lara

Robust Viability Analysis of a Controlled Epidemiological Model

Keywords: epidemics control, viability, uncertainty and robustness, Ross-Macdonald model, dengue programming equation . . . . . . . . . . . . . . . . . 11 . . . . . . . . 12 3.2 Numerical computation of robust viability kernels . . . . . . . 13 3.3 Robust viability analysis . . . . . . . . . . . . . . . . . .

Managing infectious diseases is a world public health issue, plagued by uncertainties. In this paper, we analyze the problem of viable control of a dengue outbreak under uncertainty. For this purpose, we develop a controlled Ross-Macdonald model with mosquito vector control by fumigation, and with uncertainties affecting the dynamics; both controls and uncertainties are supposed to change only once a day, then remain stationary during the day. The robust viability kernel is the set of all initial states such that there exists at least a strategy of insecticide spraying which guarantees that the number of infected individuals remains below a threshold, for all times, and whatever the sequences of uncertainties. Having chosen three nested subsets of uncertainties -a deterministic one (without uncertainty), a medium one and a large one -we can measure the incidence of the uncertainties on the size of the kernel, in particular on its reduction with respect to the deterministic case. The numerical results show that the viability kernel without uncertainties is highly sensitive to the variability of parameters -here the biting rate, the probability of infection to mosquitoes and humans, and the proportion of female mosquitoes per person. So the robust viability kernel is a possible tool to reveal the importance of uncertainties regarding epidemics control.

Introduction

Managing infectious diseases is a world public health issue. The joint dynamics of infectious agents, vectors and hosts, and their spatial movements make the control or eradication problems difficult. On top of that, uncertainties abound. Despite advances in epidemiological surveillance systems, and data on infected, recovered people, etc., there remain inaccuracies and errors. Factors such as ambient temperature, host age, social customs also contribute to uncertainty. In this paper, we focus on the impact of uncertainty on the viable control of a Ross-Macdonald epidemiological model.

Our approach departs from the widespread approach in mathematical epidemiology, where epidemic control aims at driving the number of infected humans to zero, asymptotically. Indeed, many studies on mathematical modeling of infectious diseases consist of analyzing the stability of the equilibria of a differential system (behavioral models such as SIR, SIS, SEIR [START_REF] Brauer | Mathematical models in population biology and epidemiology[END_REF]). Those studies focus on asymptotic behavior and stability, generally leaving aside the transient behavior of the system, where the infection can reach high levels. In many epidemiological models, a significant quantity is the "basic reproductive number" R 0 which depends on parameters such as the transmission rate, the mortality and birth rate, etc. Numerous works (see references in [START_REF] Hethcote | The mathematics of infectious diseases[END_REF][START_REF] Diekmann | Mathematical Epidemiology of Infectious Diseases[END_REF]) exhibit conditions on R 0 such that the number of infected individuals tends towards zero. With this tool, different (time-stationary) management strategies of the propagation of the infection -quarantine, vaccination, etc. -are compared with respect to how they modify R 0 , that is, with respect to their capacity to drive the number of infected individuals towards zero, focusing on asymptotics. However, during the transitory phase, the number of infected can peak at high values.

By contrast, we focus both on the transitory and the asymptotic regimes, where we aim at avoiding that the number of infected individuals peaks at high values. As a consequence, our approach makes no reference the concept of basic reproductive number R 0 .

In [START_REF] Lara | Viable control of an epidemiological model[END_REF], we used viability theory to analyze the problem of maintaining the number of infected individuals below a threshold, with limited fumigation capacity. The setting was deterministic, without any uncertainty in the model parameters. We said that a state is viable if there exists at least one admissible control trajectory -time-dependent mosquito mortality rates bounded by control capacity -such that, starting from this state, the resulting proportion of infected individuals remains below a given infection cap for all times. We defined the so-called viability kernel as the set of viable states. We obtained three different expressions of the viability kernel, depending on the couple control capacity-infection cap.

Several studies have applied the deterministic viable control method to managing natural resources, for example [START_REF] Béné | A viability analysis for a bio-economic model[END_REF], [START_REF] Bonneuil | Population viability in three trophiclevel food chains[END_REF] and [START_REF] Bonneuil | Viable populations in a prey-predator system[END_REF]; different examples can be found in [START_REF] Lara | Sustainable Management of Natural Resources[END_REF] in the discrete time case. Yet, few studies have undertaken a robust approach to these issues [START_REF] Béné | Contribution values of biodiversity to ecosystem performances: A viability perspective[END_REF][START_REF] Regnier | Robust viable analysis of a harvested ecosystem model[END_REF].

In this paper, we analyze what happens to the viability kernel when additional uncertain factors, that affect the dynamics of the disease, are considered. We make use of the so-called robust viability approach [START_REF] Lara | Sustainable Management of Natural Resources[END_REF], looking for vector fumigation policies able to maintain the proportion of infected individuals below a given infection cap, for all times and for all scenarios of uncertainties.

Comparison of deterministic and robust viable states can contribute to shed light on the distance between the outcomes of these two extreme approaches: ignoring uncertainty versus hedging against any risk (totally risk averse context [START_REF] Béné | Contribution values of biodiversity to ecosystem performances: A viability perspective[END_REF]). In robust viability, constraints must be satisfied even under uncertainties related to an unlikely pessimistic scenario. By contrast, reducing uncertainties to zero amounts at dressing the problem as deterministic [START_REF] Aubin | Viability theory. Systems & Control: Foundations & Applications[END_REF].

The paper is organized as follows. In Section 2, we introduce the controlled Ross-Macdonald model, discuss uncertainties and set the robust viability problem. Then, we introduce the robust viability kernel and present the dynamic programming equation to compute it. In Section 3, we provide a numerical application in the case of the dengue outbreak in 2013 in Cali, Colombia. We evaluate the impact of different sets of uncertainties on the robust viability kernel. We conclude in Section 4.

The robust viability problem

In §2.1, we introduce the Ross-Macdonald model, we discuss uncertainties and we set the robust viability problem. Then, in §2.2 we introduce the robust viability kernel. In §2.3, we present the dynamic programming equation to compute the robust viability kernel.

Ross-Macdonald model with uncertainties

To address the robust viability problem, we work with a controlled Ross-Macdonald model with uncertainties. We suppose that both controls and uncertainties change values only at discrete time steps, then remain stationary between two consecutive time steps.

Controlled Ross-Macdonald model

Different types of Ross-Macdonald models have been published [START_REF] Smith | Revisiting the basic reproductive number for malaria and its implications for malaria control[END_REF]. We choose the one in [START_REF] Anderson | Infectious Diseases of Humans: Dynamics and Control[END_REF], where both total populations (humans, mosquitoes) are normalized to 1 and divided between susceptibles and infected. The basic assumptions of the model are the following.

i) The human population (N h ) and the mosquito population (N m ) are closed and remain stationary.

ii) Humans and mosquitoes are homogeneous in terms of susceptibility, exposure and contact.

iii) The incubation period is ignored, in humans as in mosquitoes.

iv) Mortality induced by the disease is ignored, in humans as in mosquitoes.

v) Once infected, mosquitoes never recover.

vi) Only susceptibles get infected, in humans as in mosquitoes.

vii) Gradual immunity in humans is ignored.

Time is continuous, denoted by s ∈ R + . Let m(s) ∈ [0, 1] denote the proportion of infected mosquitoes at time s, and h(s) ∈ [0, 1] the proportion of infected humans at time s. Therefore, 1 -m(s) and 1 -h(s) are the respective proportions of susceptibles.

The Ross-Macdonald model, used in [START_REF] Lara | Viable control of an epidemiological model[END_REF], is the following differential system

dm ds = αp m h(1 -m) -δm , (1a)
dh ds = αp h N m N h m(1 -h) -γh , (1b)
where the parameters α, p m , p h , ξ = Nm N h , δ and γ are given in For notational simplicity, we put

A m = αp m , A h = αp h N m N h = αp h ξ . (2)
We turn the dynamical system (1) into a controlled system by replacing the natural mortality rate δ for mosquitoes in (1) by a piecewise continuous function, called control trajectory

u(•) : R + → R , s → u(s) . (3)
Therefore, the controlled Ross-Macdonald model is

dm ds =A m h(s) 1 -m(s) -u(s)m(s) , (4a)
dh ds =A h m(s) 1 -h(s) -γh(s) . (4b)
As is easily seen, the state space [0, 1] 2 is invariant by the dynamics (4).

Controlled Ross-Macdonald model sampled at discrete time steps

Time steps, separated by a period of one day, are denoted by

t = t 0 , t 0 + 1, . . . , T -1, T , (5)
where t 0 ∈ N is the initial time (day) and T ∈ N, T ≥ t 0 + 1 is the horizon.

Any interval [t, t + 1[ represents one day. Working with continuous time controls and uncertainties would make the mathematical framework more delicate, as well as the numerical calculation of the robust viability kernel. Thus our approach will be a mix, where the state follows a differential equation in which controls and uncertainties remain stationary between two consecutive time steps.

From now on, we suppose that mosquito mortality rates induced by fumigation remain stationary during every time period (day), that is,

∀t = t 0 , t 0 + 1, . . . , T -1 , u(s) = u(s ) , ∀{s, s } ⊂ [t, t + 1[ . (6)
Let us denote by Φ(M, H, u, A M , A H ) the solution, at time s = 1, of the differential system (4) with initial condition m(0), h(0) = (M, H). We obtain the following sampled and controlled Ross-Macdonald model

M (t + 1), H(t + 1) = Φ M (t), H(t), u(t), A M (t), A H (t) (7)
where

• time t runs from t 0 to T -1, with time step one day, as in ( 5),

• the state vector M (t), H(t) ∈ [0, 1] 2 represents the proportion of mosquitoes and the proportion of humans that are infected at the beginning of day t,

• the control variable u(t) represents the mosquito mortality rate applied during all day [t, t + 1[, due to application of chemical control over the mosquitoes,

• the vector A M (t), A H (t) ∈ R 2 + -whose components are infection rates for mosquito and human, respectively -carries all the uncertainties; we discuss them below.

Uncertainties

Uncertain factors in the dynamics of dengue transmission are the number of larval breeding sites, the rapid growth and urbanization of humans population, seasonal variability of mosquitoes population correlated with environmental factors such as rainfall, etc. The interactions between all these factors yield significant differences in parameters values -such as mosquito bite rate, the proportion of mosquitoes females, probabilities of infection in humans and mosquitoes -involved in the formulation of compartmental models of the process of transmission of an infectious diseases as dengue.

In the Ross-Macdonald model (1), we consider the following parameters as uncertain: bite rate (α), probability of mosquito infection (p m ), probability of infection of a human (p H ) and number of female mosquitoes per person ( Nm N h ). As detailed in Appendix A.3, we have decided to incorporate uncertainty in [START_REF] Brauer | Mathematical models in population biology and epidemiology[END_REF] through the aggregate parameters A m = αp m and A h = αp h M H in (2). Thus, the couples

A M (t), A H (t) ∈ R 2 + , ∀t = t 0 , . . . , T -1 (8)
represent uncertainty variables that affect the population of mosquitoes and humans population, respectively. We define a scenario of uncertainties as a sequence, of length T -t 0 of uncertainty couples:

A M (•), A H (•) = A M (t 0 ), A H (t 0 ) , . . . , A M (T -1), A H (T -1) ∈ (R 2 + ) T -t 0 . (9)

The robust viability kernel

Viable management in the robust sense evokes a pessimistic context. Here, it will mean that we want to satisfy the constraint of keeping the number of infected humans below a threshold imposed, whatever the uncertainties. We now give precise mathematical statement of this problem.

Control constraints

Let (u, u) be a couple such that

0 ≤ u ≤ u ≤ 1 . ( 10
)
We impose that the control variable u(t) -the mosquito mortality rate applied during all day [t, t + 1[ in ( 7) -satisfies the constraints

u ≤ u(t) ≤ u , ∀t = t 0 , . . . , T -1 . (11)

Uncertainty constraints

The two terms A M (t) and A H (t) in ( 7) encapsulate all the uncertainties affecting each population, respectively. We suppose that they can take any value in a known set S of uncertainties:

A M (t), A H (t) ∈ S ⊂ R 2 + , ∀t = t 0 , . . . , T -1 . ( 12
)
We could have allowed the set S to differ from one time to another (hence being a sequence of S(t) for t = t 0 , . . . , T -1). However, for simplicity reasons, we will only consider the stationary case.

State constraints

Let H be a real number such that

0 < H < 1 , (13)
which represents the maximum tolerated proportion of infected humans, or infection cap. Thinking about public health policies set by governmental entities, we impose the following constraint: the proportion H(t) of infected humans must always remain below the infection cap H. Therefore, we impose the state constraint

H(t) ≤ H , ∀t = t 0 , . . . , T -1, T . (14)
For example, the Municipal Secretariat of Public Health of Cali, Colombia, establishes a so-called "endemic canal" (canal endémico) with an upper bound on infected individuals.

Strategies

To define the robust viability kernel, we need the notion of strategy. A control strategy u is a sequence of mappings from states (M, H) towards controls u as follows:

u = (u t ) t=t 0 ,...,T -1 with u t : [0, 1] 2 → R . (15)
A control strategy u is said to be admissible if

u t : [0, 1] 2 → [u, u] . (16)
For any scenario A M (•), A H (•) as in [START_REF] Lara | Sustainable Management of Natural Resources[END_REF], and any initial state

M (0), H(0) = (M 0 , H 0 ) ∈ [0, 1] 2 , (17a)
a control strategy u as in [START_REF] Moré | The Levenberg-Marquardt algorithm: Implementation and theory[END_REF] produces -through the dynamics (7) -a state trajectory M (•), H(•) by the closed-loop dynamics

M (t+1), H(t+1) = Φ M (t), H(t), u t M (t), H(t) , A M (t), A H (t) , (17b)
for t = t 0 , . . . , T -1, and a control trajectory u(•) by

u(t) = u t M (t), H(t) , t = t 0 , . . . , T -1 . ( 17c
)

Robust viability kernel

Now, we are ready to lay out the definition of the robust viability kernel.

Definition 1. The robust viability kernel (at initial time t 0 ) is

V S (t 0 ) =                    (M 0 , H 0 ) ∈ [0, 1] 2
there is at least one admissible strategy u as in [START_REF] Regnier | Robust viable analysis of a harvested ecosystem model[END_REF] such that, for every scenario

A M (•), A H (•) ∈ S T -t 0 , the state trajectory M (•), H(•)
given by [START_REF] Scott | Longitudinal studies of Aedes aegypti (diptera: Culicidae) in Thailand and Puerto Rico: Blood feeding frequency[END_REF] satisfies the state constraint (14)

                   . ( 18
)
The states belonging to the robust viability kernel are called viable robust states.

We use the notation V S (t 0 ) to mark the dependency of the robust viability kernel of ( 18) with respect to the set S of uncertainties in [START_REF] Hethcote | The mathematics of infectious diseases[END_REF]. Indeed, in the sequel, we will study how V S (t 0 ) varies with S.

Dynamic programming equation

In [START_REF] Scott | Longitudinal studies of Aedes aegypti (diptera: Culicidae) in Thailand and Puerto Rico: Population dynamics[END_REF], we observe that the set of scenarios -with respect to which the robust viability kernel is defined -is the rectangle S T -t 0 . This is a strong assumption that corresponds to stagewise independence of uncertainties. Indeed, when the "head" uncertainty trajectory A M (t 0 ), A H (t 0 ) , . . . , A M (t), A H (t is known, the domain where the "tail" uncertainty trajectory A M (t + 1), A H (t + 1) , . . . , A M (T -1), A H (T -1) takes its value does not depend on the head uncertainty trajectory, as it is the rectangle S T -t+1 . This rectangularity property of the set of scenarios makes it possible to compute the robust viability kernel V S (t 0 ) of ( 18) by dynamic programming. The proof of the following result is an easy extension of a proof to be found in [START_REF] Lara | Sustainable Management of Natural Resources[END_REF].

Let us introduce the so-called constraints set

A = {(M, H)|0 ≤ M ≤ 1 , 0 ≤ H ≤ H} = [0, 1] × [0, H] (19)
and let 1 A (M, H) denote its indicator function

1 A (M, H) = 1 if (M, H) ∈ A , 0 if (M, H) / ∈ A . ( 20
)
Proposition 2. The robust viability kernel V S (t 0 ) of ( 18) is given by

V S (t 0 ) = {(M, H) ∈ [0, 1] 2 | V t 0 (M, H) = 1} , (21)
where the function V t 0 is solution of the following backward inductioncalled dynamic programming equation -that connects the so-called value functions

V T (M, H) = 1 A (M, H) , ∀(M, H) ∈ [0, 1] 2 , ( 22a
) V t (M, H) = 1 A (M, H) sup u∈[u,u] inf (A M ,A H )∈S V t+1 Φ(M, H, u, A M , A H ) , (22b) ∀(M, H) ∈ [0, 1] 2 , ( 22c
)
where t runs down from T -1 to t 0 .

We are now equipped with the concept of robust viability kernel and we dispose of a method to compute it.

Numerical results and viability analysis

Now, after having set the theory in Section 2, we will provide a numerical application in the case of the dengue outbreak in 2013 in Cali, Colombia. We introduce three nested sets of uncertainties in §3.1, we study their impact on the robust viability kernels computed in §3.2, and, finally, we discuss in §3.3 the results thus obtained.

Three nested sets of uncertainties

With the purpose of evaluating the sensitivity of the size and shape of the robust viability kernel V S (t 0 ) of ( 18) with respect to the set S of uncertainties in [START_REF] Hethcote | The mathematics of infectious diseases[END_REF], we consider the following three possible cases -S, S and S -for the set S:

L) Low case (deterministic) S = ÅM × ÅH , (23a)
M) Middle case

S = A M , A M × A H , A H , (23b)
H) High case

S = A M , A M × A H , A H . (23c)
We make the additional assumption that

S ⊂ S ⊂ S ⊂ R 2 + , (24)
that is,

0 ≤ A M ≤ A M ≤ ÅM ≤ A M ≤ A M , (25a) 0 ≤ A H ≤ A H ≤ ÅH ≤ A H ≤ A H . (25b)
By Definition 1, we easily obtain that

V S (t 0 ) ⊂ V S (t 0 ) ⊂ V S(t 0 ) , (26)
as less and less initial states can comply with the constraints in [START_REF] Scott | Longitudinal studies of Aedes aegypti (diptera: Culicidae) in Thailand and Puerto Rico: Population dynamics[END_REF] as the set S of uncertainties increases.

Numerical computation of robust viability kernels

We discuss the choice of the following numerical values in Appendix A. For the numerical calculation of the robust viability kernel in [START_REF] Scott | Longitudinal studies of Aedes aegypti (diptera: Culicidae) in Thailand and Puerto Rico: Population dynamics[END_REF], we use the dynamic programming equation ( 22) with over 60 days, that is,

t 0 = 0 , T = 60 days. ( 27
)
Following the recommendation of To remain below the "endemic canal" (canal endémico) established by the Municipal Secretariat of Public Health (Cali, Colombia), we take for the infection cap the value

H = 0.00001 = 0.001% . ( 28
)
With this value, we express that we do not want that more than 0.001% of the total population of humans be infected. We will also increase this value to H = 0.0001 = 0.01%, to test how the robust viability kernel is impacted.

In general, the image Φ(M, H, u, A M , A H ) of the dynamics does not fall exactly on one element of the 70 × 70 grid -corresponding to the cells of the state space [0, 1] × [0, H] -over which the numerical value function V t+1 is defined. We have taken a conservative stand: we put V t+1 Φ(M, H, u, A M , A H ) = 1 if and only if V t+1 = 1 for all points in the grid that surround the image Φ(M, H, u, A M , A H ).

The (discrete) robust viability kernel is defined as the set of points of the 70 × 70 grid where the indicator function V t 0 M, H is equal to 1.

Robust viability analysis

Now, as we are able to compute robust viability kernels as just seen in §3.2, we will compare them under the three nested sets of uncertainties introduced in §3.1.

Comparison between deterministic viability kernels in discrete time and in continuous time

We consider the low case (deterministic), for which we take S = {0.07660} × {0.0722} day -1 × day -1 .

(30)

These numbers were obtained in (38a) from the parameters adjusted to the 2013 dengue outbreak in Cali, Colombia, as described in [START_REF] Lara | Viable control of an epidemiological model[END_REF]. Details can be found in Appendix A.

In Figures 1 and2, we present both the deterministic viability kernel in discrete time -that is, obtained when S = S in (23a) and computed with the dynamic programming Algorithm 1-and the deterministic viability kernel in continuous time -that is, the one given by a mathematical formula in [START_REF] Lara | Viable control of an epidemiological model[END_REF]. The solid continuous line corresponds to the discrete time case and the broken starred line corresponds to the continuous time case. Figure 1: Viability kernels for the deterministic case: obtained by the dynamic programming Algorithm 1 in discrete time ( S = {0.07660} × {0.0722}, dotted line) and given by a mathematical formula in [START_REF] Lara | Viable control of an epidemiological model[END_REF] in continuous time (broken starred line), for the infection cap H = 0.00001 = 0.001% Figure 2: Viability kernels for the deterministic case: obtained by the dynamic programming Algorithm 1 in discrete time ( S = {0.07660} × {0.0722}, dotted line) and given by a mathematical formula in [START_REF] Lara | Viable control of an epidemiological model[END_REF] in continuous time (broken starred line), for the infection cap H = 0.0001 = 0.01%

We observe that the viability kernel for the deterministic case in discrete time is almost identical to the deterministic viability kernel in continuous time. Therefore, neither the sampling of controls and uncertainties -as described in §2.1 -nor the discretization process -described in §3.2lead to a degradation of the theoretical viability kernel.

Comparison between robust viability kernels as the set S of uncertainties increases

In Figures 3 and4, we display the three robust viability kernels corresponding to the three cases presented in §3.1: low case (deterministic) with S = S, middle case with S = S and high case with S = S.

For the middle case, we take

S = [0, 5] × [0, 25] day -1 × day -1 . (31)
These numbers were obtained from the ranges for the parameters adjusted to the 2013 dengue outbreak in Cali, Colombia, as described in [START_REF] Lara | Viable control of an epidemiological model[END_REF]. Details can be found in Appendix A.3.

For the high case, we doubled the right ends of each interval in (31), giving:

S = [0, 10] × [0, 50] day -1 × day -1 . (32)

Discussion

In Figures 3 and4 -both in the middle case (middle dashed line) and more markedly in the high case (left hand side dotted line) -we observe the following: a large part of the initial states (M 0 , H 0 ) that are identified as viable in the deterministic case (right hand side continuous line) are no longer viable when taking into account uncertainties.

In addition, as we expand the set S where uncertainties take their values, the gap with the deterministic case is larger and larger. In the high case, when we double the extremites of S, we observe in Figures 3 and4 that all the initial states that are below the right hand side continuous line (boundary of the deterministic viability kernel) and above the lower left hand side dotted line no longer belong to the viability kernel. Thus, enlarging the set of uncertainties can have a strong impact on the viability kernel. When we have much variability in the values of the uncertain variables (A m , A H ), the robust viability kernel is greatly reduced.

Therefore, we cannot guarantee compliance with the constraint imposed on the proportion of infected humans for all uncertainties, if we take initial conditions belonging to the deterministic kernel. Picking initial conditions within the deterministic kernel would yield optimistic conclusions with respect to infection control, when uncertainties affect the epidemics dynamics.

These observations depend on the temporal horizon and on the structure of the set of scenarios. Recall that, in §2.3, we stressed the importance of having a rectangular set of scenarios to obtain a dynamic programming equation, making it possible to recursively compute robust viability kernels. This rectangularity property of the set of scenarios corresponds to stagewise independence of uncertainties: whatever the value of the uncertainty at time t, the next uncertainties can take any values (within the rectangle). This rectangularity property contributes to having small robust viability kernels. Indeed, scenarios can display arbitrary evolutions, switching from one extreme to another between time periods (days). High and low values for the uncertainties alternate, submitting the epidemics dynamics to a strong stress, and thus narrowing the possibility of satisfying the state constraints for all times. Such contrasted scenarios deserve the label of worst-case scenarios (they represent a small fraction of the extreme points of the rectangular set of scenarios). This is why amplifying the distance between extreme uncertainties shrinks the robust viability kernel.

Conclusion

In [START_REF] Lara | Viable control of an epidemiological model[END_REF], the two authors obtained a neat expression of the deterministic viability kernel for a controlled Ross-Macdonald model. However, uncertainties abound and we wanted to assess their importance regarding epidemics control. The numerical results show that the viability kernel without uncertainties is highly sensitive to the variability of parameters -here the biting rate, the probability of infection to mosquitoes and humans, and the proportion of female mosquitoes per person. So, a robust viability analysis can be a tool to reveal the importance of uncertainties regarding epidemics control.

A Appendix. Fitting an epidemiological model for dengue

Here, we present how we identify parameters for the Ross-Macdonald model (1), and then obtain ranges for the uncertain aggregate parameters.

A.1 Parameters and daily data deduced from health reports

We introduce the vector of parameters

θ = α, p h , p m , ξ, δ ∈ Θ ⊂ R 5 + , (33)
consisting of the five parameters previously defined in Table 1. The parameter set Θ ⊂ R 5 + is given by the Cartesian product of the five intervals in the third column of Table 2.

With these notations, the Ross-Macdonald model (1) now writes

dm(s; θ) ds = αp m h(s; θ) 1 -m(s; θ) -δm , dh(s; θ) ds = αp h ξm(s; θ) 1 -h(s; θ) -γh , m(s 0 ; θ), h(s 0 ; θ) = (m 0 , h 0 ) . (34)
Notice that the rate γ of human recovery does not appear in the parameter vector θ in (33). Indeed, in the data provided by the Municipal Secretariat of Public Health (Cali, Colombia), we only have new cases of dengue registered per day; there is no information regarding how many inviduals recover daily. We choose an infectiousness period of 10 days, that is, a rate of human recovery fixed at γ = 0.1 day -1 .

Under this assumption, the daily incidence data (i.e., numbers of newly registered cases reported on daily basis) provided by the Municipal Secretariat of Public Health can be converted into the daily prevalence data (i.e., numbers of infected inviduals on a given day, be they new or not). With this, we deduce values of daily proportion of infected inviduals in the form of the set

O = r, ĥr , r = r 0 , r 0 + 1, . . . , R , (36)
where r refers to day r within the observation period of R + 1 -r 0 days, and where ĥr stands for the fraction of infected inviduals at day r. Naturally, the first couple in the set O in (36) defines the initial condition h(r 0 ; θ) = ĥ0 . Unfortunately, there is no available data for the fraction of infected mosquitoes. As mosquito abundance is strongly correlated with dengue outbreaks [START_REF] Jansen | The dengue vector Aedes aegypti : what comes next[END_REF],

we have chosen a relation m(r 0 ; θ) = 3 ĥ0 at the beginning of an epidemic outburst (other choices gave similar numerical results).

A.2 Parameter estimation

To estimate a parameter vector θ ∈ Θ in (33) that fits with the data provided by the Municipal Secretariat of Public Health we apply the curve-fitting approach based on least-square method. More precisely, we look for an optimal solution to the problem min θ∈Θ R r=r 0 h(r; θ) -ĥr 2 , R = 60 days , (37) subject to the differential constraint (34). Regarding numerics, we have solved this optimization problem with the lsqcurvefit routine (MATLAB Optimization Toolbox), starting with an admissible θ (0) ∈ Θ (the exact value stands in the second column of Table 2). The routine generates a sequence θ (1) , θ (2) . . . that we stop once it is stationary, up to numerical precision. For a better result, we have combined two particular methods (Trust-Region-Reflective Least Squares Algorithm [START_REF] Sorensen | Newton's method with a model trust region modification[END_REF] and Levenberg-Marquardt Algorithm [START_REF] Moré | The Levenberg-Marquardt algorithm: Implementation and theory[END_REF]) in the implementation of the lsqcurvefit MATLAB routine. The last column of Table 2 provides estimated values for the parameters θ = α, p h , p m , ξ, δ , and Figure 5 displays the curve-fitting results. Therefore, the aggregate parameters (2) are estimated as A m = 0.076608 day -1 , A h = 0.0722633 day -1 .

(38a)

For the natural mortality rate δ = u of mosquitoes (see Table 1), we obtain δ = u = 0.0333 day -1 .

(38b)

For the mosquito mortality maximal rate u, we take u = 0.05 day -1 . (38c)

A.3 Ranges for the uncertain aggregate parameters

From the ranges for the parameters α, p h , p m , ξ, δ displayed in the third column of Table 2, we obtain ranges for the aggregate parameters (2):

A m = αp m ∈ [0, 5] day -1 , A h = αp h ξ ∈ [0, 25] day -1 . (39)

Figure 3 :Figure 4 :

 34 Figure 3: Deterministic viability kernel ( S = {0.0766} × {0.0722}, right hand side continuous line) and robust viability kernel for the middle case (S = [0, 5]×[0, 25], middle broken line) and the high case (S = [0, 10]×[0, 50], left hand side dotted line), for the infection cap H = 0.00001 = 0.001%

Figure 5 :

 5 Figure 5: Fraction of inviduals infected with dengue, obtained by adjustment of the Ross-Macdonald model (1) (smooth solid curve) versus registered daily prevalence cases (star isolated points) during the 2013 dengue outbreak in Cali, Colombia

Table 1 .

 1 The state space is the unit square [0, 1] 2 .

	Parameter	Description	Unit
	α ≥ 0	biting rate per time unit	time -1
	ξ = N m /N h ≥ 0 number of female mosquitoes per human dimensionless
	1 ≥ p h ≥ 0	probability of infection of a susceptible
		human by infected mosquito biting	dimensionless
	1 ≥ p m ≥ 0	probability of infection of a susceptible
		mosquito when biting an infected human dimensionless
	γ ≥ 0	recovery rate for humans	time -1
	δ ≥ 0	(natural) mortality rate for mosquitoes	time -1

Table 1 :

 1 Parameters of the Ross-Macdonald model[START_REF] Anderson | Infectious Diseases of Humans: Dynamics and Control[END_REF]

Table 2 :

 2 Initial values, admissible ranges, respective source references, estimated values of parameters (numerical solution of the optimization problem (37)-(34)) and units

Acknowledgments. The authors thank the French program PEERS-AIRD (Modèles d'optimisation et de viabilité en écologie et en économie) and the Colombian Programa Nacional de Ciencias Básicas COLCIENCIAS (Modelos y métodos matemáticos para el control y vigilancia del dengue, código 125956933846) that offered financial support for missions, together with École des Ponts ParisTech (France), Université Paris-Est (France), Universidad Autónoma de Occidente (Cali, Colombia) and Universidad del Valle (Cali, Colombia). We thank the Municipal Secretariat of Public Health (Cali, Colombia) for technical support and discussions. We are indebted to the editor-in-chief and to a reviewer, whose comments contributed to improve the manuscript.

We discretize the state space [0, 1] 2 , the control space (see (38b)-(38c) for the choice of numerical values)

and the sets S and S of uncertainties in (23b)-(23c). We take a partition of 70 elements for

where the proportion M of infected mosquitoes takes its values,

• the interval [0, H], where the proportion H of infected humans takes its values,

• the interval [u, u], where the control variable u takes its values,

where the uncertainties A M (t) and A H (t) take their values, respectively, as in (23b)-(23c).

To compute the robust viability kernel V S (t 0 ) by ( 22), we implement the following dynamic programming Algorithm 1.

Algorithm 1: Dynamic programming algorithm to compute the robust viability kernel by (22)


domain: math phys

text:
Tao Guang

Francisco Xu

Adrien Chinesta

Michel Leygue

Visonneau

Guangtao Xu
email: guangtao.xu@ec-nantes.fr

Francisco Chinesta
email: francisco.chinesta@ec-nantes.fr

Adrien Leygue
email: adrien.leygue@ec-nantes.fr

Michel Visonneau
email: michel.visonneau@ec-nantes.fr

PGD for Solving the Biharmonic Equation

de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

The Biharmonic problem [START_REF] Marin | The method of fundamental solutions for inverse boundary value problems associated with the two-dimensional biharmonic equation[END_REF] has been raised in many research fields [START_REF] Erturk | Numerical solutions of 2-d steady incompressible flow in a driven skewed cavity[END_REF], such as in elasticity problem which dealing with the transverse displacements of elastic plates [START_REF] Li | The trefftz method using fundamental solutions for biharmonic equations[END_REF] and in 2D flows when using the stream function [START_REF] Montlaur | Discontinuous galerkin methods for the stokes equations using divergence-free approximations[END_REF].

The Finite Element Method -FEM - [START_REF] Gudi | Mixed discontinuous galerkin finite element method for the biharmonic equation[END_REF] and the Boundary Element Method -BEM - [START_REF] Mai-Duy | A domain-type boundary-integral-equation method for twodimensional biharmonic dirichlet problem[END_REF] can be used to solve biharmonic equations. Other methods and variants exist, e.g. [START_REF] Marin | The method of fundamental solutions for inverse boundary value problems associated with the two-dimensional biharmonic equation[END_REF].

The spectral method has been widely used in the solution of Partial Differential Equations -PDE -, in particular high order PDEs, e.g. [START_REF] Mai-Duy | A spectral collo-cation method based on integrated chebyshev polynomials for two-dimensional biharmonic boundary-value problems[END_REF].

The Chebyshev spectral collocation method [START_REF] Martinez | A chebyshev collocation spectral method for numerical simulation of incompressible flow problems[END_REF][START_REF] Li | Iterative and direct chebyshev collocation spectral methods for one-dimensional radiative heat transfer[END_REF] has been traditionally used to solved biharmonic problems. Its main advantage lies in the fact that it only needs a degree of freedom per node and it exhibits exponential convergence rates.

In this paper spectral collocation schemes are combined with the PGD technique [START_REF] Chinesta | An overview of the proper generalized decomposition with applications in computational rheology[END_REF][START_REF] Chinesta | Proper generalized decomposition of multiscale models[END_REF] that allows a separated representation of the fields involved in the model, and then in our case, transform the solution of a 2D model into the solution of few 1D problems.

THE BIHARMONIC EQUATION

We consider the biharmonic equation:

2 u = f in Ω x × Ω y = ]-1, 1[ × ]-1, 1[ (1)
where

2 = ∂ 4 ∂ x 4 + 2 ∂ 4 ∂ x 2 ∂ y 2 + ∂ 4 ∂ y 4 (2)
subjected the following boundary conditions:

u = f 1 on Γ (3)
and

∂ u ∂ n = f 2 on Γ ( 4
)
where u is the solution of the biharmonic equation. In plate theory u represents the transverse displacement and in flow simulations it represents the stream function, from which the velocity can be calculated from:

v x = ∂ u ∂ y v y = -∂ u ∂ x (5)

PGD FOR BIHARMONIC EQUATION

In this section we illustrate the biharmonic problem by using the separated representation within the PGD framework.

The aim of the method is to compute N couples of functions

(X i (x),Y i (y)), i = 1, • • • , N such that X i (x), i = 1, • • • , N and Y (y), i = 1, • • • , N are defined in 1D domains. The 2D solution reads: u(x, y) = N ∑ i=1 X i (x) •Y i (y) (6)
The weak form of problem (1) writes: Find u(x, y) verifying the boundary conditions (3) and (4) such that

Ω x Ω y u * (x, y)( 2 u(x, y) -f (x, y)) dx • dy = 0 (7)
for all the functions u * (x, y) in an appropriate functional space.

We now compute the functions involved in the separated representation. We suppose that the set of functional couples

(X i (x),Y i (y)), i = 1, • • • , n with 1 ≤ n < N
are already known (they have been previously computed ) and at the present iteration we search the enrichment couple (R(x), S(y)) by applying an alternating directions fixed-point algorithm which after convergence will constitute the next functional couple (X n+1 ,Y n+1 ). Hence at the present iteration, n + 1, we assume the separated representation

u(x, y) ≈ n ∑ i=1 X i (x) •Y i (y) + R(x) • S(y) (8)
The weighting function u * (x, y) is then assumed as

u * (x, y) = R * (x) • S(y) + R(x) • S * (y) (9)
Introducing the trial and test function into the weak form it results

Ω [R * (x) • S(y) + R(x) • S * (y)] ∂ 4 R(x) ∂ x 4 • S(y) + 2 ∂ 2 R(x) ∂ x 2 • ∂ 2 S(y) ∂ y 2 + R(x) • ∂ 4 S(y) ∂ y 4 dx • dy = Ω [R * (x) • S(y) + R(x) • S * (y)] • [ f (x, y)- ∑ n i=1 ( ∂ 4 X i (x) ∂ x 4 •Y i (y) + 2 ∂ 2 X i (x) ∂ x 2 • ∂ 2 Y i (y) ∂ y 2 + X i (x) ∂ 4 Y i (y) ∂ y 4 ) dx • dy (10)
First, we suppose that R(x) is known, implying that R * (x) = 0. Thus, equation [START_REF] Chinesta | An overview of the proper generalized decomposition with applications in computational rheology[END_REF] reads

Ω y S * (y)[α Rx S(y) + 2β Rx ∂ 2 S(y) ∂ y 2 + γ Rx ∂ 4 S(y) ∂ y 4 ]dy = Ω y S * (y)[η Rx (y) -∑ n i=1 (α i Rx Y i (y) + 2β i Rx ∂ 2 Y i (y) ∂ y 2 + γ i Rx ∂ 4 Y i (y) ∂ y 4 )]dy (11)
where

α Rx = Ω x R(x) ∂ R 4 (x) ∂ x 4 dx β Rx = Ω x R(x) ∂ R 2 (x) ∂ x 2 dx γ Rx = Ω x R(x) R(x)dx α i Rx = Ω x R(x) ∂ X 4 i (x) ∂ x 4 dx β i Rx = Ω x R(x) ∂ X 2 i (x) ∂ x 2 dx γ i Rx = Ω x R(x) X i (x)dx η Rx (y) = Ω x R(x) f (x, y)dx (12)
As the weak formulation is satisfied for all S * (y), we can come back to its associated strong form:

α Rx S(y) + 2β Rx ∂ 2 S(y) ∂ y 2 + γ Rx ∂ 4 S(y) ∂ y 4 = η Rx (y) -∑ n i=1 (α i Rx Y i (y) + 2β i Rx ∂ 2 Y i (y) ∂ y 2 + γ i Rx ∂ 4 Y i (y) ∂ y 4 ) (13)
This fourth order equation will be solved by using a pseudospectral Chebyshev method. Now, from the function S(y) just computed, we search R(x). in this case, S(y) being known, S * (y) vanishes and Eq. ( 10) reads:

Ω x R * (x)[α Sy ∂ 4 R(x) ∂ x 4 + 2β Sy ∂ 2 R(x) ∂ x 2 + γ Sy R(x)]dx = Ω x R * (x)η Sy (x)dx- Ω x R * (x)[∑ n i=1 (α i Sy ∂ 4 X i (x) ∂ x 4 + 2β i Sy ∂ 2 X i (x) ∂ x 2 + γ i Sy X i (x))]dx (14)
where α Sy = Ω y S(y)S(y)dy

β Sy = Ω y S(y) ∂ S 2 (
α Sy ∂ 4 R(x) ∂ x 4 + 2β Sy ∂ 2 R(x) ∂ x 2 + γ Sy R(x) = η Sy (x) -∑ n i=1 (α i Sy ∂ 4 X i (x) ∂ x 4 + 2β i Sy ∂ 2 X i (x) ∂ x 2 + γ i Sy X i (x)) (16)
that will be solved again by using a pseudo-spectral Chebyshev method. These two steps continue repeat until reaching the fixed point. If we denote the functions R(x) at the present and previous iteration as R p (x) and R p-1 (x), respectively, and the same for the function S(y), S p (y) and S p-1 (y), the error at present iteration can be defined from:

e = Ω x ×Ω y (R p (x) • S p (y) -R p-1 (x) • S p-1 (y)) 2 dx • dy ≤ ε (17)
where ε is a small enough parameter.

After the convergence we can define the next functional couple:

X n+1 = R and Y n+1 = S.
The enrichment procedure must continue until reaching the convergence, that can be evaluated from the error E:

E = ∆ 2 u -f (x, y) f (x, y) ≤ ε (18)
with ε another small enough parameter.

PSEUDO-SPECTRAL COLLOCATION DISCRETIZA-TION

We assume the general form of a 1D fourth order differential equation:

a d 4 u dx 4 + b d 2 u dx 2 + cu = g(x) (19)
The unknown function u(x) is approximated in

Ω x =] -1, 1[ from: u(x) = i=M ∑ i=1 α i • T i (x) (20)
where M denotes the number of nodes considered on Ω x , whose coordinates are given by

x i = cos (i -1) • π M -1 , i = 1, • • • , M (21)
The interpolants T i (x) verify the Kroenecker delta property, i.e. T i (x k ) = δ ik .

At each node k, 3 ≤ k ≤ M -2 (the remaining 4 nodes will be used for enforcing the boundary consitions) the discrete equations writes:

a • i=M ∑ i=1 α i • dT 4 i dx 4 | x k + b • i=M ∑ i=1 α i • dT 2 i dx 2 | x k + c • α k = f (x k ) (22)
When we assume that the first modes of the separated representation verified the boundary conditions (3) and (4), functions R(x) and S(y) are subjected to homogeneous Dirichlet and Neumann conditions. Thus, we should enforce u(x 1 ) = u(x M ) = 0 and du dx | x 1 = du dx | x M = 0. This conditions results in:

       α 1 = 0 ∑ i=M i=1 α i • dT i (x) dx | x 1 = 0 α M = 0 ∑ i=M i=1 α i • dT i (x) dx | x M = 0 (23)

NUMERICAL EXAMPLE

Let us consider the plate problem

The boundary conditions write

u = 0 on Γ (26)
and

∂ u ∂ n = 0 on Γ (27)
The exact solution is given by

u = 1 π 4 (1 + cos(πx))(1 + cos(πy)) (28)
which is shown in Figure 1 and that serves as reference.

The solution computed by using the separated representation within the PGD framework with M = 100 nodes in each direction is shown in Fig. 2. Figure 3 depicts the main modes involved in the separated representation. The error with respect to the reference solution (exact solution) is depicted in Fig. 4, where the error was computed at each node. As the exact solution can be expressed from 3 functional couples, the error when considering more modes is in the order of 10 -12 as noticed in Figure 5.

CONCLUSION

In this work we analyzed the possibility of using separated representations for solving high order partial differential equa-

FIGURE 1 .

 1 FIGURE 1. EXACT SOLUTION

FIGURE 2 .

 2 FIGURE 2. PGD BASED SPECTRAL SOLUTION

FIGURE 4 .FIGURE 5 .

 45 FIGURE 4. ERROR OF THE PGD SOLUTION

u(x, y) = f (x, y) in Ω x × Ω y = ]-1, 1[ × ]-1, 1[ (24)


Predicting material defects in reactive polymeric flows

Francisco Chinesta

Laboratoire de Mécanique des Systèmes et des Procédés UMR CNRS -ENSAM -ESEM 151 Boulevard de l'Hôpital, F-75013 Paris, France francisco.chinesta@paris.ensam.fr RÉSUMÉ. La simulation numérique des procédés de mise en forme des mousses nécessite tenir en compte le couplage entre la cinématique de l'écoulement et les cinétiques de polymérisation et d 'expansion avec les effets thermiques associés. La réaction de polymérisation induit une augmentation de la viscosité du fluide matrice, au même temps que le gaz produit induit l'expansion du produit. Si la réaction de polymérisation est très rapide par rapport à celle de l'expansion, le matériau devient trop rigide au cours de l'expansion, ce qui peut provoquer la rupture de la structure cellular. Par contre, si l'expansion est plus rapide que la polymérisation, à la fin de l'expansion la faible rigidité ne peut pas empêcher l'effondrement de cette structure. La simulation numérique peut aider à trouver un bon compromis ainsi qu'à contrôler le procédé afin d'éviter des autres défauts qui peuvent également apparaître.

ABSTRACT. Numerical modeling of forming processes involving the flow of chemical-reacting polymeric foams is very difficult because the flow kinematics is coupled with the thermal problem as well as with the polymerization and blowing reactions. The polymerization reaction increases the material viscosity and the expansion induces the material growing, generating a porous structure. When the polymerization takes place much faster than the foam expansion, the material becomes too rigid before the conclusion of the expansion. In this case the cellular walls can be broken with the associated consequences in the mechanical properties of the conformed pieces. On the other hand, when the foam blowing finishes before the complete polymerization, the low material consistency is not enough to preserve the final geometry, and the structure collapses when the expansion is finished. The numerical simulation can be used in order to determine a good compromise as well as to control the forming process to avoid other different defects.

MOTS-CLÉS : équation convection, suivi du front matière, mousses polymères. KEYWORDS: advection equation, flow front tracking, polymeric foams.

Introduction

Reactive polymeric systems are widely used in foam blowing or in resin transfer molding (RTM) processes. Numerical modeling of reactive processes is very difficult because the flow kinematics is coupled with the chemical reacting kinetics, which depends mainly on the reaction time (time elapsed since the reactants mixing). The kinematics -chemical kinetics coupling is a key point in foam blowing, being less important for example in usual RTM processes, where the chemical formulation can be adjusted in order to finish the mold filling before a significant rise in the resin viscosity.

The polymerization reaction increases the material viscosity and the expansion induces the material growing, generating a porous structure. When the polymerization takes place much faster than the foam expansion, the material becomes too rigid before the conclusion of the expansion. In this case the cellular walls can be broken with direct consequences in the mechanical properties of the conformed pieces. On the other hand, when the expansion finishes before the complete polymerization, the low material consistency is not enough to preserve the final geometry, and the structure collapses when the expansion is finished.

The numerical modeling of forming processes involving the flow of chemicallyreacting polymeric foams requires taking into account the different problem scales. Thus, in industrial applications a macroscopic approach is suitable, whereas the macroscopic flow parameters depend on the cellular structure: porosity, size shape and orientation of the cells, cellular walls properties, … Moreover, the shape and orientation of the cells are induced by the flow during the foam expansion (induced anisotropy). We have analyzed the influence of process conditions on the cellular shape. Elongated cells were obtained by applying a traction or a shear in the material during its expansion [START_REF] Chirivella | Experimental analysis of the induced anisotropy in foam forming processes[END_REF]. This microscopic information can be introduced in a macroscopic model by means of a homogenization technique. Thus, at each point in the material domain we can associate a characteristic volume containing some cells, whose size, shape and orientation depend on the considered point. Now, from the matrix and gas viscosities, and the cellular structure, we can compute an equivalent viscosity tensor. However, some difficulties appear in the homogenization procedure: (i) sometimes de surface tension cannot be neglected, which requires a very accurate description of the inclusion geometry as well as to consider its effect on the homogenized model; (ii) the behavior of the matrix fluid becomes quickly non-Newtonian, and its non linear character makes difficult the application of usual homogenization techniques; (iii) the reaction kinetics remain uncertain for usual industrial processes; (iv) the numerical tracking of a cell (microscopic level) which grows and deforms into the material flow is today an open problem, which in spite of some promising works [START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flows[END_REF](Sussman et al., , 1999;;Sussman and[START_REF] Sussman | An efficient interface preserving level set redistancing algorithm and its applications to interfacial incompressible fluid flows[END_REF]Sussman and[START_REF] Sussman | A coupled level set and volume-of-fluid method for computing 3D and axysymmetric incompressible two-phase flows[END_REF], important difficulties persist [START_REF] Paredes | Numerical simulation of the bubbles movement inside a flow (spanish)[END_REF]; and (v) the physical mechanism of gas diffusion, the cells interaction and fusion, and the modeling of the final stage of the process, when the material is concentrated on the cell walls, are among many others, important and unsolved difficulties found in the global modeling of foam growing.

In this way a macroscopic modeling where the reaction kinetics and the viscosity evolution are obtained experimentally is suitable from the point of view of its industrial application. Thus, such a numerical modeling could provide a visualization of the mold filling (foam volume evolution and flow front position during the expansion reaction) and will be useful to optimize the mold and process design. From that simulation we can determine the number, size and optimal position of the evacuation orifices to avoid internal overpressures (reducing the cleaning necessities as well as the material losses). Moreover, an accurate flow front tracking allows us to predict the interactions between different flow fronts when the curing reaction is too advanced (welding lines are fragile) as well as to avoid filling defects due to the big air bubbles retained between the flow front and the mold walls.

Simplified flow model

In order to obtain a macroscopic simplified flow model, able to carry out fast simulations of reacting foam flows, we consider the usual balance equations:

-The equilibrium equation neglecting the inertia and mass terms results

[1] where is the stress tensor.

-The mass conservation

[2]
where is the velocity field and the foam density. The previous equation can be rewritten as [3] When the foam expansion is homogeneous and . In this case Eq. [3] becomes [4] This function can be defined from the chemical reaction kinetics or from a simple experimental test. The first procedure, as previously discussed, remains very inaccurate because it is very difficult to know all the chemical reactions and their kinetics, in order to compute accurately the evolution of the density. Thus, we prefer, as proposed in [START_REF] Chinesta | Evaluation d'un modèle numérique simplifié de l'expansion de mousses[END_REF], to proceed from a direct determination of the density evolution.

-The constitutive equation establishes the relation between the stress tensor and the strain rate tensor

[5]

where . For compressible flows .

Eq.

[5] indicates also that the foam rheology depends on the reaction time. The simplest constitutive relation consists in [START_REF] Chirivella | Experimental analysis of the induced anisotropy in foam forming processes[END_REF] where λ denotes the volumetric viscosity, the pressure field, the unit tensor and the shear viscosity. If the volumetric viscosity contribution is neglected

[7]
where denotes the deviatoric strain rate tensor. The viscosity evolution of the matrix fluid can be measured by usual techniques [START_REF] Richter | Viscosity changes during the isothermal and adiabatic network polymerization[END_REF]. The most difficult task is the determination of the foam viscosity evolution (when the foam is considered as a pseudo-homogeneous fluid). The rheological characterization is difficult on account of the impossibility to proceed with standard rheometers which are not adapted to multiphase flows where the cell size is of the same order than the diameter or the gap of capillary or cone-plate rheometers respectively. Thus, some authors propose the use of viscosity law with the form of a mixing rule that includes as control variables the concentration in reactants, the shear rate, the temperature and the gas volume fraction [START_REF] Lefebvre | Mathematical simulation of the flow of chemicallyreacting polymeric foams[END_REF]. Other authors prefer to describe the foam rheology by means of a simple power-law model [START_REF] Deshpande | Foam phenomena in sudden expansion and contractions[END_REF] which must be identified using a non-standard device. The flow consistency and the behavior indices of the powerlaw equation could depend on time to describe the material polymerization.

-We can notice that thermal effects are not taken into account in the previous model in order to minimize the process parameters which must be identified experimentally, although it is well known that the kinetics depends significantly on the temperature field (usually the determination of the reaction kinetics is made from the temperature evolution measured during the material curing).

Numerical modeling

When the foam expansion is heterogeneous, the density depends on the spatial coordinates, and the mass conservation equation is expressed by Eq. [3]. In this case the variational formulation of Eqs. [1], [3] and [START_REF] Chirivella | Experimental analysis of the induced anisotropy in foam forming processes[END_REF] results [8] and [9] The numerical discretisation of the flow model for an homogeneous foam expansion was discussed in a former work [START_REF] Chinesta | Evaluation d'un modèle numérique simplifié de l'expansion de mousses[END_REF] From the kinematics point of view the discretisation of the flow equations governing the heterogeneous foam expansion, Eqs. [8] and [9], do not introduce particular difficulties, if we can assign to each point in the foam domain and for each time the corresponding values of , and . Usually, the volumetric viscous term can be neglected, derives directly from the density , whose time evolution is obtained following the procedure described in the previous section. Finally, in some simple rheological models, the foam shear viscosity depends on the equivalent strain rate as well as on the time, due to the polymerization and blowing kinetics.

In this form, for each point in the fluid domain at any time, we need to know, for the fluid located at that point, the time elapsed from the beginning of the chemical reactions (reaction time). Knowing that time we can assign the density and viscosity values.

In the case of homogeneous expansions we assume that the reactants are dumped in the mould simultaneously, and in consequence the reaction time is the same for any point, which corresponds with the usual Eulerian time when we consider that the reaction starts at .

However, the manufacture of big pieces requires a great amount of reactants, that must be distributed in large areas. In this case, we can not neglect the time elapsed in the reactants dumping, so that, sometimes, the blowing starts before the complete reactants dumping, giving rise to a heterogeneous foam expansion. In that case, we need to label each point of the reactant fluid volume with its dumping time , and transport this field during the foam expansion. Thus, we must solve the following linear advection equation [10] For example, if we suppose that reactants are dumped over a plate of one meter, from the left to the right , in one second, the time label could be for example . If we take the origin of times ( ) when the total amount of reactants is dumped, the reaction time will be given by .

When an Eulerian technique, using a fixed mesh, is adopted to solve Eq. [10] some difficulties appear related to the field updating in elements whose filling stage is just starting. To illustrate this fact we consider an empty finite element at time , which starts its filling process. In this case, as we are using a fixed mesh strategy, the field updating in an element i.e. the computation of at time , makes use of its value at time . Because this element starts its filling process at time , the field is not defined. If the simulation is carried out without making attention to the different field initializations, great deviations appear, mainly in the neighborhood of the flow front, whose incidence cannot be neglected. For example, in the welding of several flow fronts we require active chemical kinetics to assure the material cohesion. In order to predict this kind of material defects, we need carry out simulations taking into account the hyperbolic character of advection equation, keeping a good accuracy in the flow front treatment.

The main original contribution of this work is the establishment of a new general strategy able to carry out the field transport accurately even through the flow front. This technique minimizes the numerical diffusion and implicitly transports all fields with the flow front movement. The proposed technique can be also applied to general models involving advection or advection-diffusion equations.

Transport problems

As discussed in the last section, a transport problem appears associated with the advection of the field T, which is governed by the following linear advection equation [10]. If we consider a first order discontinuous finite element discretisation then [11] where is the value of T in the element at time , fluid volume coming to the element the upstream element (where the value of T is given by ), is the fluid volume leaving the element and is the volume of .

A first problem appears suddenly if we consider the time t for which an element starts its filling process from its upstream element . To illustrate this limitation we consider the simplest situation of an incompressible fluid, i.e. .

In this case, with , we obtain [12] which establishes that at is dominated by solution existing in the element at time , at which this element is empty, and consequently the property is not defined. question is: what is the value to assign to T before the beginning of its filling process ? An accurate and consistent discrete form of the conservation equation, can be obtained by considering a discontinuous finite element discretisation of the equation which results from the addition of Eq. [10] multiplied by the fluid fraction I , whose evolution is governed by the advection equation [START_REF] Chinesta | Evaluation d'un modèle numérique simplifié de l'expansion de mousses[END_REF] [13] with Eq. [13] multiplied by the field T, as proposed in [START_REF] Chinesta | Some difficulties in the flow front treatment in fixed mesh simulations of composites forming processes[END_REF] [14] which results in [15] where is the volume of fluid coming to from the upstream element ;

Conclusions

From these results we can conclude about the accuracy and low numerical diffusion of the numerical strategy proposed in the discretisation of Eq. [10]. The results are significantly improved, in relation to standard discontinuous finite element discretisations, mainly in the neighborhood of the flow front. In this way the proposed strategy can be used to compute accurately the reaction times in welding areas in 3D foam expansion process simulation, to predict defects in the material welding.

  We consider a numerical example involving a welding line development from two reacting flows with different reaction times. Fig.1shows the foam presence function and the velocity field at three steps in the mold filling. Fig.2depicts the transport of the dumping time. The actualization of the field T is carried out by using Eqs.[15] and[16].

Figure 1 .

 1 Figure 1. Foam presence and velocity field during the expansion and welding

Figure 2 .

 2 Figure 2. Transport of the dumping time and welding line development