Olivier Pironneau
email: olivier.pironneau@sorbonne-universite.fr

Pierre-Henri Tournier
email: pierre-henri.tournier@sorbonne-universite.fr

Simulation of the 3D Radiative Transfer with Anisotropic Scattering for Convective Trails

In two earlier publications we have proposed a numerical algorithm to solve the RTE in 3D based on an integro-differential formulation and iterations using H-matrices to speed-up the computation of the integrals. In this article we show how the method can be extended to handle an important class of non-isotropic scattering for the atmosphere. An application to the Earth atmosphere is given.

Introduction

Airplane combustion engines generate CO 2 (carbon dioxide), a gas which stays in the Earth atmosphere for ∼ 500 years before being partially absorbed by the oceans.

Currently used fuels and potential future substitutes like hydrogen, synthetic or agricultural fuel combustion generate also NO x (Nitrous oxides). Typically a nitrous oxide stays in the atmosphere ∼ 100 years.

Carbon dioxide, nitrous oxides, methane, ozone and water vapor are "green house gases", (GHG), in the sense that atmospheric and past climates measurements indicate that they are most likely responsible for the additional heat received by Earth (global warming).

Air travel -around 100.000 flights/day -produces 3% of CO 2 (or 12% of what is produced by all transport vehicles). It is still small, but the contribution of airplanes to global warming is expected to rise from today's 0.024W/m 2 to 0.084W/m 2 by 2050 [START_REF] Marquart | Future development of contrail cover, optical depth, and radiative forcing: impacts of increasing air traffic and climate change[END_REF].

A cloud similar to a cirrus, called "contrail (short for condensation trail) may appear in the airplane wake if atmospheric pressure and temperature are on the left side of a threshold curve, and further left of that pressure versus temperature curve the contrail will persist [START_REF] Paoli | Contrail modeling and simulation[END_REF]. Several studies claim that these clouds have a warming effect perhaps 3 times stronger than the one caused by the combustion gases [START_REF] Lee | The contribution of global aviation to anthropo-genic climate forcing for 2000 to 2018[END_REF].

In the past few years the authors have worked to see if these claims could be validated by a numerical simulation of the fundamental equations of physics for these problems [START_REF] Bardos | Radiative transfer for the greenhouse effect[END_REF], [START_REF] Golse | Radiative transfer in a fluid[END_REF], [START_REF] Golse | Stratified radiative transfer for multidimensional fluids[END_REF], [START_REF] Golse | Stratified radiative transfer in a fluid and numerical applications to earth science[END_REF]. we approach the problem from an applied mathematics point of view, with convergence error estimation and computational efficiency in mind.

The formation of a contrail has been simulated by solving the Navier-Stokes equations (NSE) with chemistry for the engine exhaust and phase change for the ice formation in the airplane wake (see [START_REF] Paoli | Contrail modeling and simulation[END_REF]); it is the right approach to understand the mechanisms of the formation of contrails. Once the cloud is sufficiently developed one may study its effect on the absorption and scattering coefficients of the Radiative Transfer Equations (RTE) in the Earth atmosphere (see the data in [START_REF] Emde | The libradtran software package for radiative transfer calculations[END_REF]) and then solve the RTE-NSE system. The numerical simulations of RTE in one dimension is the object of intense research [START_REF] Cahalan | Bringing together the most advanced radiative transfer tools for cloudy atmospheres[END_REF], [START_REF] Emde | The libradtran software package for radiative transfer calculations[END_REF], but mostly without coupling with a temperature equation or NSE.

In [START_REF] Golse | Radiative transfer for variable 3d atmospheres[END_REF] and [START_REF] Pironneau | Reflective conditions for radiative transfer in integral form with h-matrices[END_REF] we have proposed a numerical algorithm to solve the RTE in 3D based on an integro-differential formulation and iterations using Hmatrices to speed-up the computation of the integrals. In this article we show how the method can be extended to handle an important class of nonisotropic scattering for the atmosphere.

It is known that the preferred direction of scattering in clouds is the initial direction. It is modelled by an anisotropic probability of scattering (also called phase function) from ω to ω , p ν (ω, ω ) (see Figure 1).

The Radiative Transfer Equations

When molecular viscosity and wind convection are ignored, the temperature T in a medium exposed to electromagnetic waves satisfies the RTE as explained in [START_REF] Pomraning | The equations of Radiation Hydrodynamics[END_REF]. It involves a frequency dependent radiation intensity field ω pν (θ) ω Fig. 1 This polar plot shows the probability p(θ) that a photon in the direction ω scatters in the direction ω :

pν (θ) = 1 + 1 2 ω • ω = 1 + β cos(θ)
where θ is the angle (ω, ω ). The solid curve is for β = 1 2 , the dashed curve is for β = 0.75 and the dotted circle is the isotropic case is β = 0. I ν (x, ω) at position x in the physical domain Ω and in each direction ω. For all {x, ω, ν}

∈ Ω × S 2 × R + , ω • ∇ x I ν + ρκ ν a ν I ν - 1 4π S 2 p ν (ω, ω )I ν (ω )dω = ρκ ν (1 -a ν )[B ν (T ) -I ν ], ∞ 0 S 2 ρκ ν (1 -a ν )[B ν (T ) -I ν ]dωdν = 0, (1)
where S 2 is the unit sphere,

B ν (T ) = 2 ν 3 c 2 [e ν kT -1]
is the Planck function,

, c, k are the Planck constant, the speed of light in the medium and the Boltzmann constant. The absorption coefficient κ ν := ρκ ν , where ρ is the medium density, comes from nuclear physics, but for our purpose it is seen as the percentage of radiation absorbed per unit length. The scattering albedo is a ν ∈ (0, 1) and p ν (ω, ω ) is the probability that a ray in direction ω scatters in direction ω .

With appropriate boundary conditions, existence of solution has been established by [START_REF] Porzio | Application of accretive operators theory to evolutive combined conduction, convection and radiation[END_REF] and [START_REF] Golse | Radiative transfer in a fluid[END_REF]. In the latter, convergence of the following scheme was proved:

ω • ∇I n+1 ν + κ ν I n+1 ν = κ ν 4π S 2 p ν (ω, ω )I n ν (ω )dω + κ ν (1 -a ν )B ν (T n ), ∞ 0 S 2 κ ν (1 -a ν )[B ν (T n+1 ) -I n+1 ν ]dωdν = 0.
(

) 2
By the maximum principle for the first equation and the monotony of T → B ν (T ),

I n ≺ I n-1 , T n ≺ T n-1 =⇒ I n+1 ≺ I n , =⇒ T n+1 ≺ T n ,
where a ≺ b means a(x) ≤ b(x), for all x but not a(x) = b(x) everywhere.

Hence by choosing T 0 = 0 and I 0 = 0, the positivity of T 1 and I 1 implies that T 1 ≺ T 0 , I 1 ≺ I 0 so that , by induction, a strictly increasing sequence is obtained. Similarly a decreasing sequence is obtained if we manage T 1 T 0 , I 1 I 0 . As it was observed in [START_REF] Golse | Stratified radiative transfer in a fluid and numerical applications to earth science[END_REF] that the increasing sequences converge faster we will focus on that one only.

Integral Formulation

Consider an anisotropic scattering density

p ν (ω, ω ) = 1 + β ν ω • ω .
Denote by Γ the boundary of Ω. One must find the radiation intensity I ν (x, ω) at all points x ∈ Ω, for all directions all ω ∈ S 2 and all radiation frequencies ν ∈ R + , by solving the radiative transfer equations (RTE):

ω • ∇ x I ν + κ ν I ν = κ ν (1 -a ν ) (B ν (T ) + β ν ω • K ν ) + κ ν a ν J ν , (3)
J ν (x) := 1 4π S 2 I ν dω, K ν (x) := 1 4π S 2 ωI ν (x, ω)dω , (4)
∞ 0 κ ν (1 -a ν )(J ν -B ν (T ))dν = 0 , (5)
I ν (x, ω) = R ν (x, ω)I ν (x, ω -2(n • ω)n) + Q ν (x, ω), on Σ := {(x, ω) ∈ Γ × S 2 : ω • n(x) < 0 }. (6)
In ( 6), Q ν is the radiation source and R ν is the portion of radiation which is reflected by the boundary; n(x) is the outer normal of Γ at x. κ ν > 0 and a ν ∈ [0, 1] are the absorption and scattering coefficients; in general they depend on ν and x.

The general solution of (3) is

I(x, ω) = I(x Σ (x, ω), ω)e -τx,ω 0 κ(x-ωs)ds + τx,ω 0 e -s 0 κ(x-ωs )ds S ν (x-ωs, ω)ds,
where τ x,ω is the length |x -x Σ (x, ω)|, and S ν denotes its right-hand side,

S ν (x, ω) = κ ν (x)(1-a ν (x)) (B ν (T (x)) + β ν (x)ω • K ν (x))+κ ν (x)a ν (x)J ν (x) . (7)

Solution in Absence of Reflective Boundaries

Let us denote [x,y] κ := |y-x| 0 κ(x+s(y-x))ds. When R ν = 0, the following holds:

Proposition 1 J ν (x) = S E ν (x) + J [S ν ](x), (8)
S E ν (x) = 1 4π Γ Q ν (y, y-x |y-x| ) [(y -x) • n(y)] - |y -x| 3 e -[x,y] κ dΓ (y), (9)
J [S ν ](x) = 1 4π Ω e -[x,y] κ |y -x| 2 S ν (y, y -x |y -x| )dy. ( 10
)
Averaging ( 7) on S 2 after multiplication by ω leads to

K ν (x) := 1 4π S 2 ωI(x, ω)dω = K E ν (x) + K[S ν ](x) with K E ν (x) := 1 4π S 2 ωI(x Σ (x, ω), ω)e -τx,ω 0 κ(x-ωs)ds dω = 1 4π Γ (y -x)Q ν (y, y-x |y-x| ) [(y -x) • n(y)] - |y -x| 4 e -[x,y] κ dΓ (y), K[S ν ](x) := 1 4π S 2 τx,ω 0 ωe -s 0 κ(x-ωs )ds S ν (x -ωs, ω)dsdω = 1 4π Ω (y -x)S ν (y, y-x |y-x| ) e -[x,y] κ |y -x| 3 dy. ( 11
)

Iterative Method

In [START_REF] Golse | Radiative transfer for variable 3d atmospheres[END_REF] it was shown that, in absence of K ν , i.e. β ν = 0, the following is monotone and convergent. Its extension to β ν > 0 is: Algorithm 1. Initialize J ν , K ν and T (by zero, for instance).

2. Compute S ν by [START_REF] Golse | Radiative transfer for variable 3d atmospheres[END_REF].

3. Update J ν by ( 11),( 9), [START_REF] Golse | Stratified radiative transfer in a fluid and numerical applications to earth science[END_REF] and K ν by [START_REF] Hackbusch | A sparse matrix arithmetic based on h-matrices. part i: Introduction to h-matrices[END_REF]. 4. Update T by solving (5).

Extension to Reflective Conditions (RC)

As in [START_REF] Pironneau | Reflective conditions for radiative transfer in integral form with h-matrices[END_REF], consider boundary condition [START_REF] Emde | The libradtran software package for radiative transfer calculations[END_REF]. Proposition 1 can be extended to case with RC and in the case of non-multiple reflection, it becomes:

Proposition 2
The same iterations are proposed with

J ν (x) = S E ν,1 (x)+S E ν,2 (x)+ J [S ν ](x), K ν (x) = K E ν,1 (x)+K E ν,2 (x)+ K[S ν ](x), (12)
S E ν,1 (x) = 1 4π Γ Q ν (y, y-x |y-x| ) [(y -x) • n(y)] - |y -x| 3 e -[x,y] κ dΓ (y), S E ν,2 (x) = M n=1 1 4π Γ R ν (x n , x-x n |x-x n | )Q ν (y, x n -y |x n -y| ) × [(x n -y) • n(y)] -e -[x,x n ]∪[x n ,y] κ |x n -y| (|x -x n | + |x n -y|) 2 dΓ (y). J [S ν ](x) = 1 4π Ω e -[x,y] κ |y -x| 2 + M n=1 e -[x,x n ]∪[x n ,y] κ (|x -x n | + |x n -y|) 2 R ν (x n , x-x n |x-x n | ) S ν (y, y-x |y-x| )dy. (13)
K E ν,1 (x) = 1 4π Γ (y -x)Q ν (y, y-x |y-x| ) [(y -x) • n(y)] - |y -x| 3 e -[x,y] κ dΓ (y), K E ν,2 (x) = M n=1 1 4π Γ (x -x n )R ν (x n , x-x n |x-x n | )Q ν (y, x n -y |x n -y| ) × [(x n -y) • n(y)] -e -[x,x n ]∪[x n ,y] κ |x n -y| (|x -x n | + |x n -y|) 2 dΓ (y), K[S ν ](x) = 1 4π Ω y -x |y -x| 3 e -[x,y] κ + M n=1 (x -x n ) e -[x,x n ]∪[x n ,y] κ (|x -x n | + |x n -y|) 2 R ν (x n , x-x n |x-x n | ) S ν (y, y-x |y-x| )dy. ( 14
)
where M is the number of RC boundaries and x n is the point of reflection on the boundary of the ray going from x to y via x n .

All equations are discretized using a P 1 finite element framework. Some integrals have a singular kernel, so a careful quadrature should be used.

In [START_REF] Golse | Radiative transfer for variable 3d atmospheres[END_REF], [START_REF] Pironneau | Reflective conditions for radiative transfer in integral form with h-matrices[END_REF] a strategy is explained to compute all the integrals as a matrix vector product where the matrices are hierarchical compressed Hmatrices [START_REF] Bebendorf | Approximation of boundary element matrices[END_REF], [START_REF] Bebendorf | Hierarchical Matrices[END_REF], [START_REF] Hackbusch | A sparse matrix arithmetic based on h-matrices. part i: Introduction to h-matrices[END_REF].

In the grey case (κ ν independent of ν) we need 5+M matrices. This could be taxing in computer memory, but not in computing time because the core of the method is N ln N where N is the number of finite element vertices. In the non grey case we need P(5+M) matrices where P is the number of values taken by the numerical approximation of κ ν .

The Stratified One Dimensional Approximation

Let us consider the Earth's atmosphere submitted to black body radiations from the sun at temperature T S and earth at T E . The atmosphere is thin and the Sun is far. As a first approximation the ground is locally flat and has a negligible relief. Then all quantities are only function of the altitude z and independent of x and y. This, so called, stratified approximation of the RTE is (see [START_REF] Bardos | Radiative transfer for the greenhouse effect[END_REF])

(µ∂ τ + κ ν )I ν (τ, µ) = κ ν (1 -a ν )B ν (T (τ )) + 1 2 κ ν a ν 1 -1 p ν (µ, µ )I ν (τ, µ )dµ , (15)
I ν (0, µ) = r ν I(0, -µ) + µQ + ν , I ν (Z, -µ) = µQ - ν , 0 < µ < 1 , (16)
∞ 0 κ ν (1 -a ν )B ν (T (τ ))dν = ∞ 0 κ ν (1 -a ν )J ν (τ )dν, (17)
with

Q + ν = Q E B ν (T E ), Q - ν = Q S B ν (T S ), J ν (τ ) := 1 2 1 -1 I ν (τ, µ)dµ, ( 18
)
where µ is the cosine of ω with the vertical direction and τ is the optical thickness:

τ = z 0 ρ(z )dz .

RTE with Non-Isotropic and Rayleigh Phase Function

The scattering function is p ν (µ, µ ) = 1 + β ν µ in a cloud between altitude Z m and Z M and at high altitude z > Z M it is the Rayleigh phase function

p ν (µ, µ ) = 3 16 (3 -µ 2 ) + 3 16 (3µ 2 -1)µ 2 .
It is combined into one formula with two given altitudes and possibly frequency dependent functions b ν (τ ) and

β ν (τ ) in [0, 1], p ν (µ, µ ) = b ν + β ν µ + 3 16 (1 -b ν )(3 -µ 2 + (3µ 2 -1)µ 2 ). ( 19
)
Observe that p ν (µ, µ ) ≥ 0 and 1 2 1 -1 p ν (µ, µ )dµ = 1 . Keeping [START_REF] Porzio | Application of accretive operators theory to evolutive combined conduction, convection and radiation[END_REF] as the defining equation for T , given I, the problem becomes

(µ∂ τ + κ ν )I ν (τ, µ) = S ν (τ ) := R ν (τ ) + P ν (τ )µ 2 , I ν (0, µ) = r ν I(0, -µ) + µQ + ν , I ν (Z, -µ) = µQ - ν , 0 < µ < 1 , (20)
with 21) and [START_REF] Porzio | Application of accretive operators theory to evolutive combined conduction, convection and radiation[END_REF]. 21), then R ν (τ ) and P ν (τ ). 4. Then, compute T n+1 by solving [START_REF] Porzio | Application of accretive operators theory to evolutive combined conduction, convection and radiation[END_REF].

K ν (τ ) := 1 2 1 -1 µI ν (τ, µ)dµ , L ν (τ ) = 1 2 1 -1 µ 2 I ν dµ , R ν (τ ) :=κ ν a ν ( 7 16 b + 9 16 )J ν (τ ) + βK ν (τ ) -3 16 (1 -b)L ν (τ ) + κ ν (1 -a ν )B ν (T (τ )) , P ν (τ ) := 3 16 (1 -b)(-J ν (τ ) + 3L ν (τ )) (
Algorithm 1. Initialize J 0 ν (τ, µ) = 0, K 0 ν (τ, µ) = 0, L 0 ν (τ, µ) = 0 and T 0 (τ ) = 0. 2. Compute I n+1 by (20) with S n ν . 3. Update J n+1 ν (τ, µ), K n+1 ν (τ, µ), L n+1 ν (τ, µ) by (

Implementation

Recall the definition of the exponential integrals

E p (X) := 1 0 e -X/µ µ p-2 dµ , X > 0 . ( 22
)
For i = 3, 4, 5, let

S i (ν, τ ) := 1 2 E i (κ ν τ )Q + (τ ) + (-1) i-1 2 (E i (κ ν (Z -τ )) + r ν E i (κ ν (Z + τ )))Q -(τ ), F i (τ, t) := 1 2 E 3 (κ ν |τ -t|) + rν 2 E 3 (κ ν |τ + t|). ( 23
)
Proposition 3 With (21) and (23) the quantities needed by the algorithm above are given by

J n+1 ν (τ ) =S 3 (ν, τ ) + 1 2 Z 0 F 1 (τ, t)R n ν (τ )dt + 1 2 Z 0 F 3 (τ, t)P n ν (τ )dt , K n+1 ν (τ ) =S 4 (ν, τ ) + 1 2 τ 0 F 2 (τ, t)R n ν (t)dt + 1 2 τ 0 F 4 (τ, t)P n ν (t)dt -1 2 Z τ F 2 (τ, t)R n ν (t)dt -1 2 Z τ F 4 (τ, t)P n ν (t)dt , L n+1 ν (τ ) =S 5 (ν, τ ) + 1 2 Z 0 F 3 (τ, t)R n ν (τ )dt + 1 2 Z 0 F 5 (τ, t)P n ν (τ )dt . (24)
Proof For clarity let us assume that r ν = 0. For the general case see [START_REF] Golse | Stratified radiative transfer for multidimensional fluids[END_REF]. Applying the method of characteristics shows that

I ν (τ, µ) = µe -κν τ µ Q + ν 1 µ>0 + |µ|e -κν (Z-τ ) |µ| Q - ν 1 µ<0 + 1 µ>0 τ 0 e -κν (τ -t) µ κν µ S ν (t)dt + 1 µ<0 Z τ e -κν (t-τ ) |µ| κν |µ| S ν (t)dt . ( 25
)
Therefore with k = 0 or k = 2,

1 -1 µ k I ν (τ, µ)dµ = 1 -1 |µ| k+1 e -κν τ µ Q + ν 1 µ>0 + e -κν (Z-τ ) |µ| Q - ν 1 µ<0 dµ + 1 -1 |µ| k-1 1 µ>0 τ 0 e -κν (τ -t) µ κ ν S ν (t)dt + 1 µ<0 Z τ e -κν (t-τ ) |µ| κ ν S ν (t)dt dµ = 1 0 |µ| k+1 e -κν τ µ Q + ν dµ + 0 -1 |µ| k+1 e -κν (Z-τ ) |µ| Q - ν dµ + 1 0 |µ| k-1 τ 0 e -κν (τ -t) µ κ ν S ν (t)dt + 0 -1 |µ| k-1 Z τ e -κν (t-τ ) |µ| κ ν S ν (t)dtdµ = 1 0 |µ| k+1 e -κν τ µ Q + ν + e -κν (Z-τ ) |µ| Q - ν dµ + τ 0 1 0 |µ| k-1 e -κν (τ -t) µ κ ν (R ν (t) + P ν (t)µ 2 )dµdt + Z τ 1 0 |µ| k-1 e -κν (t-τ ) |µ| κ ν (R ν (t) + P ν (t)µ 2 )dµdt = E k+3 (κ ν τ )Q + ν + E k+3 (κ ν (Z -τ ))Q - + τ 0 E k+1 (κ ν (τ -t))R ν (t)dt + Z τ E k+1 (κ ν (t -τ ))R ν (t)dt + τ 0 E k+3 (κ ν (τ -t))P ν (t)dt + Z τ E k+3 (κ ν (t -τ ))P ν (t)dt.
(26) For k = 1, the same computation gives

1 -1 µI ν (τ, µ)dµ = 1 -1 µ 2 e -κν τ µ Q + ν 1 µ>0 -e -κν (Z-τ ) |µ| Q - ν 1 µ<0 dµ + 1 -1 1 µ>0 τ 0 e -κν (τ -t) µ κ ν S ν (t)dt -1 µ<0 Z τ e -κν (t-τ ) |µ| κ ν S ν (t)dt dµ = 1 0 µ 2 e -κν τ µ Q + ν dµ - 0 -1 µ 2 e -κν (Z-τ ) |µ| Q - ν dµ + 1 0 τ 0 e -κν (τ -t) µ κ ν S ν (t)dt - 0 -1 Z τ e -κν (t-τ ) |µ| κ ν S ν (t)dtdµ (27) = 1 0 µ 2 e -κν τ µ Q + ν -e -κν (Z-τ ) |µ| Q - ν dµ + τ 0 1 0 e -κν (τ -t) µ κ ν (R ν (t) + P ν (t)µ 2 )dµdt - Z τ 1 0 e -κν (t-τ ) |µ| κ ν (R ν (t) + P ν (t)µ 2 )dµdt = E 4 (κ ν τ )Q + ν -E 4 (κ ν (Z -τ ))Q - + τ 0 E 2 (κ ν (τ -t))R ν (t)dt - Z τ E 2 (κ ν (t -τ ))R ν (t)dt + τ 0 E 4 (κ ν (τ -t))P ν (t)dt - Z τ E 4 (κ ν (t -τ ))P ν (t)dt.
(28)

Numerical Simulations

At the top of the atmosphere, we are told, the solar radiation intensity is 340W/m 2 but the atmosphere reflects 30% of it and our RTE does not handle volumic reflection. At 10 or 12km the sunlight is diffused, so we apply an intensity proportional to ω 3 . Paris being at latitude 45 0 , the radiation is divided by √ 2 . Hence an effective number is Q sun = 162W/m 2 . And yet a good portion is spent for the evaporation of water and into convection. Arbitrarily we keep only half 80W/m 2 .

Recall the scalings defined in [START_REF] Golse | Stratified radiative transfer in a fluid and numerical applications to earth science[END_REF]. The frequencies and the temperatures are scaled as follows (primes label scaled variables), ν = 10 -14 ν, T = 10 -14 k h T = 10 -14

Q + ν = 80 = Q S ∞ 0 B 0 ν 3 e ν T S -1 (10 14 dν ) = Q S 1.4744 • 10 6 π 4 15 1.209 4 .
This leads to Q S = 4 • 10 -6 .

Similarly Earth is at T E = 288/4798 = 0.06 and emits about 350W/m 2 , but some is used for the evaporation of water so we have kept only 300W/m 2 :

300 = Q E ∞ 0 B 0 ν 3 e ν T E -1 (10 14 dν ) = Q E 1.4744 • 10 6 π 4 15 0.06 4 = Q E 124.3.
This leads to Q E = 2.41. The air density is 1 -3 4 z. The earth reflective albedo is set at r ν = 0.3. Isotropic scattering is set at a ν = 0.3.

κ ν is given by the Gemini experimental program. To measure the temperature perturbation due to an increase of CO 2 (resp NO x )in the atmosphere, we set κ ν = 1 in the range (3/16, 3/14) (resp. (3/7, 3/5).

Isotropic Scattering

Two tests were performed with the stratified 1D approximation, isotropic scattering, no cloud and Gemini κ ν + Rayleigh scattering equal to 0.3 for ν > 3 and τ > 0.7. One computation is with an infrared source Q E at τ = 0 only (similar to the situation at night) and another with a sunlight source Q S at τ = Z only.

Results are on Figures 3,4. Computations are done with 3 different atmospheres. One is conform to the Gemini measurements, the other two imitate an addition of CO 2 or methane.

One can observe a small additional greenhouse effect due to addition of CO 2 but a cooling effect at high altitude. The opposite is observed for addition of NO x . Only the Sun radiation is taken as a source.

Anisotropic Scattering in a Cloud

Next, still with the stratified approximation, a cloud is added between altitude Z m = 7000m and Z M = 9000m with anisotropic scattering coefficient

β = aν 2 , b = 1, a ν = 0.3 4 (Z M -Zm) 2 (z -Z m ) + (Z M -z) + .
Above Z M Rayleigh scattering 0.31 ν>3 , replaces the anisotropic scaattering.

Results are shown on Figure 5, 6 and 7. The strong effect of the cloud is shown on Figure 6. Figure 7 shows how important it is to take all the details of κ ν into account.

Conclusion

The integro-differential formulation of the RTE and its solution by iterations on the source has been extended here to handle anisotropic scattering. The iterative part of the method is O(N ln N ), thanks to an efficient use of H-matrices.

The precision is good enough to evaluate the effect of sensitive parameters for the study of contrails. Most of the time the stratified 1D approximation should suffice, but in complex cases with high relief the 3D formulation is needed.

Fig. 10 Ground and high altitude temperature in an atmosphere receiving sunlight and Earth infrared radiation. In addition 0.3% of the sunlight is reflected in the normal direction of the ground and 0.7% if the ground is covered by snow.

Fig. 2

 2 Fig. 2 Visual representations of the H-matrices surface-to-volume on the left and volumeto-volume on the right, for the computation of the temperature in the Chamonix valley. The level of compression is shown by light green-to-red and the numbers correspond to the rank of the truncated SVD approximation of each block.

Fig. 3 Fig. T vs z for 3

 33 Fig. 3 Temperature T versus altitude z for 3 values of ν → κν . Only the earth radiation is taken as a source. (night time).

Fig. 5 TFig. 6 Fig. 7

 567 Fig. 5 T versus z computed with a cloud

  .

Fig. 9

 9 Fig.9Temperature above the city of Chamonix versus altitude in a grey case. The growth of T from the ground up is partially due to quadrature errors, partially to the parabolic effect of mountains.

  1.381 • 10 -23 6.626 • 10 -34 T = • 10 16 10 42 = 1.4744 • 10 -8 .

	T 4798	.
	Planck's function is written as
	B 0 2.998 2 Consequently the computed intensity, I ν , is the physical intensity divided by ν 3 e ν T -1 , with B 0 = 2h 2 × 6.626 • 10 -34 c 2 10 42 =
	B 0 . The scaled Sun temperature is T S = 5800/4798 = 1.209. According to
	(18),

∞ 0

Finally, in Table 1 total radiative intensities ∞ 0 J ν (τ ) are shown at τ = 0 and τ = Z in a variety of situations to see the effect of clouds, scattering, earth albedo. In the presence of CO 2 the total radiative energy at ground level (resp. τ = Z), is 9.91316 (resp. 5.72292). With NO x it is 9.55564 (resp. 5.8434). These numbers must be compared with the first number on the left in the first (resp 2 nd ) row in Table 1.

A 3D formulation with a Stratified Part

To implement in 3D the Earth albedo (reflective condition by the ground) in all generality is quite complicated. To simplify the computation we consider the reflection of the sunlight only and we use the linearity, with respect to the sources, of the equation for I ν , for a fixed temperature field.

For a portion Ω of atmosphere between the ground g = (x 1 , x 2 , x 3 = g(x 1 , x 2 )), g ≥ 0, and the lower stratosphere

The solution is I ν = I E ν + I S ν with I E ν computed with the condition I ν | g only and the source term in T , and I S ν computed with I ν | Z only and T = 0:

ν the stratified approximation is justified. Furthermore, it can be precomputed beforehand.

To account for the fact that a fraction r ν (x) of the sunlight is reflected by the Earth surface, we consider a partial reflective condition

so, neglecting scattering, a ν = 0 (see also [START_REF] Bardos | Remarks on the Radiative Transfer Equations for Climatology[END_REF])

.

Therefore, with ω • n < 0 and x on the ground,

For a flat ground, it is equal to ω

and its angle average on S 2 is 1 2 Q S B ν (T S )r ν (g)E 3 (κ ν Z) so that it makes sense to solve the problem with the following conditions,

The valley of Chamonix is considered in the same physical conditions as above, i.e. with a cloud between 7000m and 9000m. However we have not implemented the anisotropic scattering yet, so β = 0.


Nicolas Peyret
email: nicolas.peyret@supmeca.fr

Gaël Chevallier

Jean-Luc Dion

Gael Chevallier

Dynamic Damping in Joints: Multiscale Model Taking into Account Defects in a Nominally Plane Surface

Keywords: Jointed structures, damping, dry friction, hysteretic damping

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

Introduction

Estimating the damping remains a important challenge in the computation of the vibration levels of jointed-structures. Computer aided design and finite element analysis are used to predict respectively inertia and stiffness with fairly good accuracy, but damping is often badly estimated, meaning that vibration amplitudes are also wrongly predicted. Linear damping is well taken into account in simulation tools, but the Softwares to study friction-induced damping (FID) remains quite sparse, whereas joints, such as welded points, bolted joints or rivets, are widely used to link the parts of mechanisms and structures.

Essentially, there are two distinct dissipation mechanisms. First, damping can be caused by friction between sliding parts, see for instance [START_REF] Berthillier | Blades forced response analysis with friction dampers[END_REF], [START_REF] Whiteman | Displacement-dependent dry friction damping of a beam-like structure[END_REF], [START_REF] Poudou | Modeling and analysis of the dynamics of dry-friction-damped structural systems[END_REF], [START_REF] Korkmaz | Analyse dynamique expérimentale d'une poutre amortie par un amortisseur à frottement sec[END_REF], Dion et al. [2013] and [START_REF] Chevallier | Dynamic friction evolution during transient sliding[END_REF]. In this case, there is a macroscopic motion between the parts of the joint. Second, damping can be caused by the structure itself through a joint between the parts. In such case, there is no macroscopic motion between the parts. Damping is induced by micro sliding or partial sliding. This means that there are sticking zones and sliding zones between the parts. Many learned authors [START_REF] Goodman | Analysis of slip damping with reference to turbine blade vibration[END_REF][START_REF] Beards | The damping of structural vibration by rotational slip in joints[END_REF][START_REF] Pian | Structural damping of a simple built-up beam with riveted joints in bending[END_REF][START_REF] Ungar | The status of engineering knowledge concerning the damping of builtup structures[END_REF][START_REF] Ahmadian | Identification of bolted lap joints parameters in assembled structures[END_REF][START_REF] Metherell | Instantaneous energy dissipation rate in a lap joint -uniform clamping pressure[END_REF][START_REF] Rogers | Energy dissipation through joints: Theory and experiments[END_REF][START_REF] Goyder | Development of a method for measuring damping in bolted joints[END_REF][START_REF] Peyret | Micro-slip induced damping in planar contact under constant and uniform normal stress[END_REF][START_REF] Peyret | Energy dissipation by micro-slip in an assembly, analytic and experimental approach[END_REF] have designed experimental testing devices to highlight the fact that the damping value in assembled structures is much greater than the inherent damping value of the material. To measure and identify the dissipated energy during both dissipation mechanisms, total sliding and partial sliding, Dion et al. [2013] have proposed a specific experimental setup.

The present work starts from the idea of a previous paper, see [START_REF] Peyret | Micro-slip induced damping in planar contact under constant and uniform normal stress[END_REF]. In the latter, the testing device was a clamped-clamped bending beam built with three parts assembled together as shown in Fig. 1. The beam material is Steel Fig. 1. View of the specimen. On the left: the three parts of the beam. On the right: the testing device assembled. Faces Γ U are screwed on an anti-vibration marble. Strain gauges are glued on both sides (2) of the beam to measure the normal load. Piezoelectric patches (3) are glued on both sides of the beam to apply a torque T piezo to excite the first vibration mode. A quick clamp is used to apply the normal load before the tests on faces Γ N .

C35, with a Young modulus E = 220 GPa and a Poisson Coefficient ν = 0.3. Cut was obtained by electrical discharge machining. The joints have been geometrically placed in order to avoid coupling between the normal load and the tangential load inside the joint, see Appendix B. In [START_REF] Peyret | Nonlinear dynamic behavior modelling of a planar friction interface in a structure assembly[END_REF], the shape of the beam has been improved in order to obtain the greatest possible damping and to minimize the sensitivity of the results to manufacturing errors, see Fig. 1. The final design of this new testing device was presented in Dion et al. [2013]. The latter paper presents an original experimental procedure based on stop-sine excitation, see also [START_REF] Heller | Experimental identification of nonlinear dynamic properties of built-up structures[END_REF], obtained through piezoelectric transducers. The authors justify that their excitation method is more efficient that Dirac (shock) or Heaviside (release of a static load) excitations. Signal processing tools are proposed to post-process the instantaneous damping and frequency of the first mode. Finally, this work [Dion et al., 2013] experimentally demonstrates the presence of micro-slip damping, by comparing the vibration behavior of a monolithic part and a built-up structure with friction-joints. This comparison is the best way to determine the added damping, [START_REF] Segalman | A modal approach to modeling spatially distributed vibration energy dissipation[END_REF]. Figure 2 sums up the different results of these studies. These curves are repositioned at the intrinsic damping of the assembly (not evaluated in the models) with the goal to highlight only the damping due to the interfaces. In blue, the theoretical loss factor considering that the interface is modelled as a nominally flat surface, this result show no damping for the lower amplitudes. This result was obtained in [START_REF] Peyret | Micro-slip induced damping in planar contact under constant and uniform normal stress[END_REF], where the authors compute the energy losses assuming that the interfaces are perfectly plane. The modeling of the interface by a flat surface cannot show that there are energy losses for the smallest oscillations.

To our opinion, the surface defects generated during the manufacturing process are responsible for the gap between experimental and numerical results. These defects are challenging the assumptions of flat contact interfaces and uniform normal stress fields. The contact area between both parts is geometrically defined by the nominal contact area A 0 , i.e., the area where both solids seem to be in contact. In fact, because of the geometrical defects, the real contact area A R only consists of the top surfaces of the defects. A R is naturally smaller than the geometric surface, and the real contact pressure p R is much higher than the nominal pressure p 0 . To compute A R and p R , one of the most commonly used model is the Greenwood-Williamson model [START_REF] Greenwood | Energy dissipation through joints: Theory and experiments[END_REF][START_REF] Greenwood | Surface roughness and contact: An apology[END_REF]. The height and the shape of the roughness are the main parameters to compute A R and p R . In the Greenwood-Williamson model, the asperities are assumed to be uniformly distributed in the horizontal plane, without interaction between them, and their height is defined by a Gaussian distribution. This law seems to best match the experimental measurements and has been the subject of many mathematical developments. Moreover, the asperities are expected to have spherical ends, permitting each contact asperity to be approximated as an hertzian contact. Such models focus on real surfaces with geometrical defects of 3rd order standards [NF EN ISO 4287-4288, 1998] and [NF EN ISO 12085, 1998]. For more complexity, one can refer to: [START_REF] Greenwood | The contact of two nominally flat rough surfaces[END_REF], who took into account the roughness of different forms and other kinds of statistical distribution; Goryacheva and Dobychin [1991], who showed interactions between asperities; [START_REF] Tsukizoe | On the mechanism of contact between metal surfaces: The penetrating depth and the average clearance[END_REF] and [START_REF] Mccool | Predicting microfracture in ceramics via a microcontact model[END_REF], who improved the uniformly distribution on the surface; [START_REF] Pugliese | Rough contacts between actual engineering surfaces: Part II. Contact mechanics[END_REF], who took into account the plasticity of the contact. [START_REF] Bjorkland | A random model for micro-slip between nominally flat surfaces[END_REF] and [START_REF] Eriten | Surface roughness effects on energy dissipation in fretting contact of nominally flat surfaces[END_REF] combined the asperity model of [START_REF] Greenwood | Energy dissipation through joints: Theory and experiments[END_REF] with Mindlins treatment of the tangential contact problem to define a statistical estimate of the resulting dissipation, a similar approach was used by [START_REF] Bureau | Elasticity and onset of frictional dissipation at a non-sliding multi-contact interface[END_REF] to interpret the results of experimental dissipation measurements under microslip conditions. In previous works, the authors proposed, also similar approaches with models and experimental correlations to take into account the energy losses in the contact. A rheological contact model based on the extended Greenwood model with micro-contacts and statistical distributions was developed and studied. Experimental results and simulations are compared in order to assess and discuss the model, [Bouchaala et al., 2013[Bouchaala et al., , 2014]].

The geometrical uncertainties have an important effect on the loss factor and the stiffness of the structure. To minimize the uncertainties on simulations results, we can minimize the geometrical uncertainties by geometrical specifications of the roughness. But this solution strongly depends on the statistical distribution of the roughness, if robust results is expected, it could be possible to define an geometry with lesser influence on results. In this work, only the geometrical defects taken into account are of 1st and 2nd (order standards [NF EN ISO 4287-4288, 1998] and [NF EN ISO 12085, 1998]), form defects and waviness defects. The purpose of this study is to find a simple model with a quite good correlation to the experimentation. This, to be able to propose which kind of defect is possible to produce to reach a given level of damping. First, the aim is to compute the real contact zones and the real pressure field in the joints. Second, to calculate the energy losses, the tangential behavior is introduced. This is done thanks to Mindlin's model. Finally, the model is generalized to a multi-contact interface.

Model of the Interface

In order to model the contact of a real surface, a multiple spherical contact model, based on Hertz' and Mindlin's theories, is proposed.

Parameters of multiple spherical contact model

The geometrical defects can be described as the members of a statistical population. The experimental observations (Fig. 8) show that the zones in contact can be considered as spheres and their population can be divided into two families, see Fig. 3: • G-Spheres, the largest, that ensure the rigidity of the interface. They are the first to come into contact during the normal loading. They might not slide completely, thus they are only in partial sliding. • P -Spheres, the smallest. They are the last to come into contact. They might slide partially or completely.

The interface loading is characterized by the normal displacement δ n of the two normal planes constituting the interface. Before the loading, the difference between "G" spheres "P" spheres the heights of the two families of spheres is called h. There are m G-spheres of radius ρ G and n P -spheres of radius ρ P . Figure 3, on the right, describes the balance between the global normal load N , the tangential load T and the distribution of local normal and tangential forces N P , T P , N G , T G of both families of asperity.

N = mN G + nN P T = mT G + nT P (1)

Physical model for one asperity

The distribution of the normal load N is a function of the normal displacement δ n and is defined, for a single spherical contact, by Hertz' theory [START_REF] Hertz | On the contact of elastic solids[END_REF] as follows:

δ n = 1 ρ 3N ρ 4E * 2 3
(2)

with 1 E * = 1-ν 2 1 E1 + 1-ν 2 2 E2 and 1 ρ = 1 R1 + 1 R2
, where E i , ν i and R i are respectively the modulus of elasticity, the poisson's ratio and the radius of curvature of the two solids in contact.

Once the normal load N is applied, a tangential displacement δ is applied between the two surfaces in contact, see Fig. 4. The Mindlin theory [START_REF] Mindlin | Compliance of elastic bodies in contact[END_REF][START_REF] Mindlin | Effects of an oscillating tangential force on the contact surfaces of elastic spheres[END_REF], gives the expression of the displacement induced by partial sliding of the contact surface, as

δ = 3µN 16 G * a 2 -c 2 a 3 (3) with G * = 2-ν 2 1 G1 + 2-ν 2 2 G2 and G j = Ej 2(1+νj )
From the previous equations, Johnson [1987] offers the solution for a spherical contact, initially loaded with a constant normal load N , and then subjected to a tangential load T oscillating between ±T max with T max < µN, see Fig. 5. This loading is divided into three phases:

• Phase 0, curve OA: initial tangential loading, T ∈ [0 T max ] δ = µN 16a G * 1 -1 - T µN 2 3
(4)

• Phase 1, curve ABC: cycle tangential unloading, T ∈ [T max -T max ] δ = µN 16a G * 2 1 - T max -T 2µN 2 3 -1 - T max µN 2 3 -1 (5) • Phase 2, curve CDA: cycle tangential loading, T ∈ [-T max T max ] δ = µN 16a G * -2 1 - T max + T 2µN 2 3 + 1 - T max µN 2 3 + 1 (6)
The area defined by the cycle gives the dissipated energy for a cycle (5).

W diss = T dδ (7)
This energy can be expressed as a function of the tangential displacement δ:

W diss = 9µ 2 N 2 10a G * 1 -1 - δ δ lim 5 2 - 5 6 1 -1 - δ δ lim 3 2 2 - δ δ lim (8)
where δ lim = 3µN G * 16a defines a limit value of tangential displacement for each type of spherical contact for which there is the passage from the macro-slip phase to the total-slip one, see Mindlin's theory [START_REF] Mindlin | Compliance of elastic bodies in contact[END_REF][START_REF] Mindlin | Effects of an oscillating tangential force on the contact surfaces of elastic spheres[END_REF]. This limit is defined on a loading cycle including macro-slip, see Fig. 5.

Physical model for multiple contacts

Using the force relations given in Eq. ( 1) and the Hertz' model of Eq. ( 2), the normal load on the interface constituted with several contacts can be expressed as:

N = 4 3 E * mρ 1 2 G δ 3 2 n + nρ 1 2 P (δ n -h) 3 2 (9) if δ n -h > 0 otherwise N = 4 3 E * mρ 1 2 G δ 3 2 n ( 10
)
In the following, we are going to consider the case with both P-Spheres and G-Spheres in contact, see Eq. ( 9). Once the normal load is applied, a tangential displacement δ is then imposed. The latter is assumed to be the same for every spherical contact. According to the value of δ, three phases can be highlighted, see Fig. 6:

• "Generalized partial slip" (GPS): G-and P-spheres do not slip totally.

• "Pseudo-partial slip" (PPS) : G-spheres remain in partial slip and P-spheres slip totally.

• "Total slip" (TS) : G-and P-spheres slip totally. The distribution of the tangential load is function of the normal displacement δ n :

• For the phase of generalized-partial slip (GPS)

T = mµN G 1 -1 - δ δ G lim 3 2 + nµN P 1 -1 - δ δ P lim 3 2 (11)
• For the phase of pseudo-partial slip (PPS)

T = mµN G 1 -1 - δ δ G lim 3 2 + nµN P (12)
The loading cycle of the multiple spherical contacts is plotted in Fig. 6. The comparison of the later with the shape of the cycles given by a single spherical contact, of Fig. 5, shows a new behavior during the loading especially when the P-Spheres begin to slip completely. Equation ( 7), which defines the energy dissipated by a single spherical contact, with Eqs. ( 11) and ( 12), are used to define the energy dissipated by a multiple spherical contact:

W diss = 9µ 2 10 G * (w G + w P ) (13)
where w G is always in partial-slip

w G = m N 2 G a G 1 -1 - δ δ G lim 5 2 - 5 6 1 -1 - δ δ G lim 3 2 2 - δ δ G lim (14)
and w P is defined according to the phase: in generalized partial slip (GPS):

w P = n N 2 P a P 1 -1 - δ δ P lim 5 2 - 5 6 1 -1 - δ δ P lim 3 2 2 - δ δ P lim ( 15
)
in pseudo-partial slip (PPS):

w P = µN P δ (16)

Analysis of the model parameters influence

It is interesting to analyze the evolution of the friction-induced vibration-damping through the loss-factor η. The mathematical definition of the loss factor is

η = W diss 2πW ext (17)
where W diss is the disipated energy in the interface and W ext is the work of external forces. The loss-factor is a function of the multiple-contact model parameters.

The values of these parameters can be debatable: the number of P and G-spheres depends on the manufacturing process, the friction coefficient can change during the vibration motion according to wear, temperature, etc. In order to visualize the effects of each parameter, a study of the sensibility of the parameters has been carried out, see Table 1. The results have been summarized in Fig. 7 and discussed in the following section. As said in Sec. 2.1, the model contains only two kinds of spheres, under the assumption that the G-spheres provide the rigidity of the interface and remain in partial slip whereas the P-spheres slip totally when δ becomes large enough. The identification of the contacts can be tricky. For the studied specimen, Sec. 3 gives the characteristics of the surfaces : 5 G-spheres and 18 P-spheres can be considered from the measurements. The other parameters are chosen from our knowledge, see line 1 of Table 1. This set of parameters provides the reference case which will be shown on Fig. 7. On X-axis, the displacement υ at the center of the beam is plotted, while on Y -axis the loss factor η is shown.

The first design parameters of the parametric study are n and m, respectively the number of P-spheres and G-spheres. As n increases, the maximum damping also increases, however the vibration amplitude υ has to be greater, as it can be seen in Fig. 7(b). On the contrary, see Fig. 7(c), when m increases, the damping decreases: this is due to the fact that δ n is fixed, which implies that the global normal load, increases when m increases. The number of G-spheres for a given δ n directly influences the resultant normal load and the tangential stiffness of the interface. If δ n increases, as h is fixed, the compression of the P-spheres no longer allows the phase of macro-slip: this implies that the damping increases for greatest the values of υ, see Fig. 7(d).

In fact, the parameters are not independent, when the normal displacement increases, P-spheres are more compressed and therefore they have the same behavior than G-spheres: no phase of macro-slip. However others spheres, previously not in contact, find themselves slightly compressed and dissipate energy by macro-slip. So increasing the normal displacement leads to an increase of the number m of G-spheres. Figure 7(e) presents the evolution of the loss factor, increasing both the normal displacement δ n and the number m of G-spheres. While the number of P-spheres n remains constant, see Table 1. The damping function remains slightly dependent on the amplitude. To conclude this analysis, it is interesting to consider the effect of a variation of the friction coefficient µ. This leads to an evolution of the threshold of macro-slip for the P-spheres, see Fig. 7(f).

Experimental Characterization of the Joint Surfaces of the Test Bench

In order to obtain the real contact surfaces, 3D measurements were made, with a ZEISS OI 442 machine, see Fig. 8. On each surface, a million points were measured.

After filtering and re-sampling, a mesh constituted by 1200 × 800 pixels has been obtained. The contact takes place at the highest points over the imaginary perfect plane. To choose a cloud of points, defining the contact-area, only the "higher" points, 10% of the set of points, were retained. Figure 8 shows that the surfaces present irregularities on which contacts. To determine the influence of the defects on the vibration-damping, the surface irregularities are modeled as sphere caps, see Fig. 9, where df and ∅(m) are respectively the height and the diameter of the sphere cap. R, the radius of curvature is calculated from df and ∅ using the following expression:

R = ∅ 2 4 + df 2 2df ( 18
)
In accordance with the model, two types of spheres are identified on each surface, see Table 2. Following a classical assumption, see [START_REF] Greenwood | Energy dissipation through joints: Theory and experiments[END_REF] and [START_REF] Greenwood | Surface roughness and contact: An apology[END_REF], in a contact between two surfaces, one is considered to be flat, while the other one includes all the defects. The numbers m = 5 and n = 18 have been evaluated by counting the number of defects of both sides of the interface. There are 5 that belong to the big defects category and 18 that belong to the small defect category. Thus, during the normal loading of the interface, the G-spheres ensure the rigidity of the interface and are the first and most compressed, whereas the P-spheres are only slightly compressed.

Numerical Simulation of the Test Bench

The energy dissipated is plotted in Fig. 10. As function of the displacement δ, i.e., Eq. ( 13), the figure shows three behaviours: for small displacements W diss is nil; for intermediate displacements W diss is governed by Eq. ( 15) of the GPS and for high displacements W diss is governed by Eq. ( 16) of the PPS. Overall, W diss is governed by a non-continuous polynomial function with an exponent between 0 and 2. Taking into account the external force applied on the structure to obtain the dynamic behavior, it is possible to define the loss factor of the structure:

η = W diss 2πW ext ( 19
)
where the work of external forces W ext can be expressed as function of the displacement at the center of the beam υ:

W ext = 1 2 K m υ 2 (20)
W ext is governed by a quadratic function. Thus the loss-factor is first nil due to W diss , then it increases and it tends to zero when δ tends to infinite, see Fig. 7(a). The latter is plotted for a 460 N pre-load and for the data summarized in Table 1.

The transition from GPS to PPS leads to a drastic change in the evolution of the loss factor. To compare the experimental results with the results provided by the model, one needs to isolate the contribution of interfaces studied from the other sources of damping. The two structures, assembly and monolithic, have the same geometry and the same boundary conditions, it is thus possible to consider two distinct sources of damping. We distinguish the dissipation induced by the interfaces (friction joints) f d2 from the dissipation provided by all the others sources f d1 (clamps, intrinsic damping of the material, etc.):

W diss 0
mϋ + kυ = f d1 ( υ, υ) + f d2 ( υ, υ) (21)
Thus, the dissipated energy can be written as

W diss = (f d1 ( υ, υ) + f d2 ( υ, υ))dυ (22)
Both energies are separable by subtraction. It is considered that W d1 can be identified by a test on the monolithic beam and W d2 can be calculated by subtracting the damping rate of the monolithic beam (see Fig. 11). The results of the simulations, with plane surfaces in interfaces [START_REF] Peyret | Micro-slip induced damping in planar contact under constant and uniform normal stress[END_REF] and with two types of defects in interfaces, are positioned at the intrinsic damping of the assembly besause this contribution are not taken into account in these models. It is possible to compare models simulation to the experimental results [Dion et al., 2013]. This comparison of the experimental results with those given by the model highlights a better agreement than in [START_REF] Peyret | Micro-slip induced damping in planar contact under constant and uniform normal stress[END_REF]. Taking into account geometrical defaults in interfaces, even if the parameters of the retained geometry for model is quite simple, it allows a good agreement. The choosen model has a good agreement with experimental results, except for small amplitudes. This difference is probably due to the assumption considering only two types of spheres of different heights. In fact this model does not take into account the progressivity of the spheres height in the interfaces and the fact that the irregularities of the surface outcrop the contact and dissipate some energy at small vibration amplitudes.

Conclusion

Unlike models that take into account the interface as a geometrically perfect plane [START_REF] Peyret | Micro-slip induced damping in planar contact under constant and uniform normal stress[END_REF], taking into account form and waviness defects of the surface allows to evaluate damping at the interfaces also at the smallest vibration amplitudes. A model that takes into account form and waviness surface defects was developed in order to model the evolution of damping as a function of strain amplitude. This new approach is a multi-spherical contact modeling of the interface, based on the behavior of a single spherical contact developed by Hertz and Mindlin. This multi-scale model considers the presence of two types of contacts: contacts which ensure the rigidity of the interface and contacts which dissipate energy. This modeling allows to obtain a better correlation between the simulation results and the experimental measurements. Further work aims to link macro-models used for vibration analysis [START_REF] Quinn | Modal analysis of jointed structures[END_REF], with local parameters. According to the function of the studied assembly, this work allows to design contact surfaces. In order to design a stiff assembly, the model can be used to define the maximum admissible surface defect. On the contrary, in order to design a dissipative assembly, the model can be used to define controlled or desired surface defects.

Appendix B Design of the Experimental Device

The principle of the experimental device Fig. 1 was justified in previous works, see [Dion et al., 2013] and [START_REF] Peyret | Micro-slip induced damping in planar contact under constant and uniform normal stress[END_REF]: it is a clamped-clamped beam excited on its first bending mode. The beam is built up with three parts linked by two planar joints. A normal pre-load N is applied on both extremities of the beam before there are clamped to the ground. The three parts remain in contact even during the bending motion due to the normal pre-load N and the friction between the counter-parts. This design allows avoiding coupling between pre-load induced normal stresses in the planar joints and the vibration-induced shear stresses due to the zero bending moment at the location of the joints, see Fig. B.1. In the previous paper (see [START_REF] Peyret | Micro-slip induced damping in planar contact under constant and uniform normal stress[END_REF]), the beam was studied under a quasi-static assumption (QSA), i.e., the normal and shear stresses were computed using the definition of the static bending moment when the beam was loaded with a centered transverse load. A more recent work (see [Dion et al., 2013]) showed that QSA leads to errors in defining the zerobending-moment section and for the amplitude of the transverse load in this section. By using the definition of the first bending mode (MA), the bending moment is equal to zero in the cross sections located at x = 89.6% of L/4, see Fig. B.1. Moreover the transverse load is 20% higher for MA than for QSA. By improving the shape of the beam, it is possible to define the position of the cutting section which leads to a zero bending moment. This section has to be well placed even if the experimental mode shape is slightly different from the theoretical one due to machining errors. The constant section over the length of the beam does not allow obtaining this robustness. The difference between MA and QSA illustrates this lack of robustness in the design well.

The beam has to be as thick as possible in order to maximize the damping of the structure. Machining errors can lead to badly positioned cutting-sections and therefore to coupling between the normal load and the vibration motion. In order to ensure robustness in the face of machining errors, the authors propose a design that allows obtaining a nil bending moment over a non-zero length Fig. B.2. In order to obtain uniformly distributed normal stress over the cutting plane, the thickness of the beam changes smoothly over its length. Finally, particular attention was given to obtaining the highest frequency for the second bending mode in order to avoid coupling between the first two modes. All these design constraints led the authors to define the geometry of the

Fig. 2 .

 2 Fig. 2. Experimental results: damping versus the normal load and versus the displacement magnitude. Both the monolithic and the cut beam were tested. The latter was tested under four normal load levels in the range 100-1000 N. Simulation results: damping versus the displacement magnitude for a normal load value of 500 N.

Fig. 3 .

 3 Fig. 3. Left: Simplified description of the population of geometrical defects. Right: Balance between local and global forces.

Fig. 4 .

 4 Fig. 4. Schematic elastic spheres pressed by a constant normal force N and subjected to an oscillating tangential force T .

Fig. 6 .

 6 Fig. 6. Loading cycle for multiple contacts model. Tangential load T (N ) versus δ (m) displacement in the contact.

Fig. 7 .

 7 Fig. 7. Results of the parametric study.

Fig. 8 .Fig. 9 .

 89 Fig. 8. 3D measurements of test-bench interfaces (in mm) from ZEISS OI442.

Fig. 10 .

 10 Fig. 10. Energy dissipated for a multi-contact interface. W diss (J) versus displacement δ(m).

  Fig. B.1. Bending moment (left) and transverse force (right) belonging to the cross sections of the beam studied. The dashed lines indicate the cutting planes located at x = 100 and 300 mm under the quasi static assumption and at x = 89.6 mm and 310.4 mm under the modal assumption. These results were obtained considering a analytical beam model under a Timoshenko assumption.

Fig

  Fig. B.2. Bending moment plotted over the beam length. Particular attention should be given to the width of the zero moment zone which makes the device robust to machining errors.

Fig

  Fig. B.3. Left: New design, 3D view of half of the beam. The clamp is located on Vol. 1. Volume 2 is cut in order to make the joint. Volume 3 has been designed to ensure the normal stresses are uniformly distributed. Right: Normal stress distribution under the first mode shaped prescribed displacement. The normal stresses are zero in the neighborhood of the cutting plane.

  Fig. B.3: the plane of the beam is symmetrical; it has two clamped blocks and two interface areas . The blocks around the two interface zones are designed to maximize the width of the zone where the bending moment is nil, see Fig. B.2.

Table 1 .

 1 Parametric study: Values of the parameters.

	m	n	N(N )	δn(10 -7 m)	h(10 -7 m)	µ	Figure
	5	18	460	2.5	2.4	0.2	7a (Ref)
	5	4, 9	450, 454	2.5	2.4	0.2	7b
		36, 72	474, 501
	1, 3	18	103, 282	2.5	2.4	0.2	7c
	8, 11		728, 997
	5	18	572, 710	2.7, 2.9	2.4	0.2	7d
			869	3.1
	3, 5	18	250, 460	2.3, 2.5	2.2, 2.4	0.2	7e
	7		716	2.7	2.6
	5	18	460	2.5	2.4	0.15, 0.18	7f
						0.25, 0.3

Table 2 .

 2 Model parameters.

	Spheres	Number	df (m)	∅(m)	R(m)
	G	n = 5	1.7 10 -6	1.7 10 -2	23.8
	P	m = 18	1.2 10 -6	1.6 10 -2	26.7
	h		2 .4 10 -7 m

New simulation results Experimental results from Dion & al. [2013] Simulation results from Peyret & al. [2010] Corrected Loss Factor %

  Fig. 11. Comparison between experimental results, simulation results obtained under perfectly flat surfaces assumption and simulation results obtained taking under account surface defects.

		1.4					0.8	Cut Beam, N=500N
		1.2						Cut Beam, N=500N
	η [%]						0.6
	Measured Loss Factor	0,6 0,8 1	0.2	0.4	0.6	0.8	0.4 0.2	Cut Beam, N=500N
		0,4		Intrinsic damping of the
				structure coming from the
				material,
				clamps, etc
		0,2
		0	0.2	0.4	0.6	0.8	1	1.2	1.4
					Amplitude µ m( )


Mouna Fradi
email: mounafradi@yahoo.com

Raoudha Gaha
email: raoudha.gaha@gmail.com

Abdelfattah Mlika
email: abdelfattah.mlika@gmail.com

Faïda Mhenni
email: faida.mhenni@supmeca.fr

Jean Yves Choley
email: jean-yves.choley@supmeca.fr

Design of an Electronic Throttle Body Based on a New Knowledge Sharing Engineering Methodology

Keywords: Collaborative design, Mechatronic system, Methodology, ETB, KARREN

To minimize the number of iterations and correction returns while designing a system, sharing crucial parameters and data between different actors is needed. Since mechatronic systems are considered complex because of their multi-disciplines, their design requires collaborative work to ensure the sharing of parameters between contributors from different domains. This paper proposes a new methodology based on the collaborative design to choose the architecture of a mechatronic system. This methodology is structured around three main phases: the pre-collaboration phase, the collaboration phase and the postcollaboration phase. The proposed methodology has been validated by applying it on a mechatronic system called Electronic Throttle Body (ETB). In order to share the different established activities and capitalize knowledge, KARREN (Knowledge Acquisition and reuse for Robust Engineering) a collaborative tool from DPS (Digital Product Simulation) a French company, is used in this project to help us to choose the appropriate architecture of the Electronic Throttle Body.

Introduction

Mechatronics is defined as a synergic integration of mechanical, electronics, information systems and computer systems in the design of a product in order to increase and optimize the functionality. This alliance requires a simple and fast data exchange between different disciplines to obtain a product of high quality, low price and a short time. For this reason, the concept of collaborative engineering has replaced sequential engineering to improve communication between project actors and to ensure robust and efficient sharing of knowledge. In this context, knowledge Based-Engineering (KBE) is a technology able to capture and reuse knowledge to reduce costs and time of development. This technology has large application such as knowledge reuse and collaboration support [START_REF] Mcharek | Knowledge sharing for mechatronic systems design and optimization[END_REF]. In this direction, this paper proposes a new methodology based on a collaborative design to choose the most appropriate architecture of a mechatronic system using KBE application to support collaboration between actors.

The paper is organized in several sections. Section 2 presents a review of related works. In Sect. 3, our methodology is presented in details. In Sect. 4, the proposed methodology is illustrated on a case study. Section 5 gives some results and Sect. 6 summarizes and concludes the paper.

2 State of the Art

Mechatronic Systems Development Methodologies

As previously mentioned, mechatronic systems are known with their multidisciplinary nature, which increases the complexity of the design. In order to decrease this complexity, various groups of researchers to develop mechatronic systems with high quality and low cost have done significant works. Four types of development methodologies are discussed in this section.

First, the methodologies based on the optimization with integration of control were presented as methodologies to develop mechatronic systems. In this context, Zhang et al. [START_REF] Zhang | Integrated design of mechanical structure and control algorithm for a programmable four-bar linkage[END_REF] proposed an integrated approach for mechatronic design to facilitate the control system. To improve this work [START_REF] Li | Design for control-a concurrent engineering approach for mechatronic systems design[END_REF] developed a design framework called Control For Design (DFC). This approach emphasizes the importance of designing the control parameters simultaneously with the structural parameters even though it is possible to change the controller parameters after the system is built. Another methodology has been proposed by Villarreal-Cervantes et al. [START_REF] Villarreal-Cervantes | Structure-control mechatronic design of the planar 5r 2dof parallel robot[END_REF] to formulate the mechatronic design problem of the 5R 2Dof parallel robot and its Proportional-Integral-Derivative (PID) controller. They integrated a nonlinear dynamic optimization problem. This work concerns parallel planar robots. Nevertheless, some applications require high performance positioning systems with a spatial parallel mechanism. In this framework, the contribution of Lara-Molina et al. [START_REF] Dumur | Structure-control optimal design of 6-DOF fully parallel robot[END_REF] takes into account the structure and control design variables simultaneously. They proposed an optimal design methodology for Stewart-Gough robot. This robot is a three-dimensional mechanism where its movable platform is connected to the fixed base by six legs. This methodology aims at maximizing the accuracy of positioning to optimize the overall performance of the robot. In these previous works, effective approaches have been introduced. However, the improvement of the system performance based on the possibility of modifying controller parameters was overlooked, which represents an important criteria in design evaluation.

Methodologies based on the design evaluation are another type, which aims at developing mechatronic systems. Moulianitis et al. [START_REF] Moulianitis | A model for concept evaluation in design-an application to mechatronics design of robot grippers[END_REF] proposed an evaluation index that includes the complexity, flexibility and intelligence of mechatronic systems. In addition, Hammadi et al. [START_REF] Hammadi | A new multicriteria indicator for mechatronic system performance evaluation in preliminary design level[END_REF] proposed a multi-criteria performance indicator called Mechatronic Design Indicator (MDI) as a neuronal network of radial basis function. This indicator provides accurate information for decision making while decreasing design time. However, this methodology overlooks other mechatronic design requirements such as risk assessment, cost, etc. As an improvement to this work, Mohebbi et al. [START_REF] Mohebbi | Mechatronic multicriteria profile (MMP) for conceptual design of a robotic visual servoing system[END_REF] developed a new multi-criteria profile. This profile includes five criteria, which are cost, flexibility, complexity, intelligence and reliability of mechatronic systems. These criteria can be determined using a specific mathematical function. This methodology is applied to the conceptual design of a 6-Dof robotic manipulator. Following this work, Moulianitis et al. [START_REF] Moulianitis | A new index based on mechatronics abilities for the conceptual design evaluation[END_REF] proposed a new index includes seven criteria, which are configurability, adaptability, dependability, decisional autonomy, motion ability, perception ability and interaction ability. In formulating this index, the cost and the complexity are taken into account. The proposed index is applied to the design of firefighting robot in order to evaluate its design solutions. However, all these approaches are not extended for use in the detailed design phase.

Due to their flexibility in design representation, bond graphs and evolutionary algorithms are considered as an efficient method for system representation in detailed phases as well as preliminary phases. Seo et al. [START_REF] Seo | Toward a unified and automated design methodology for multi-domain dynamic systems using bond graphs and genetic programming[END_REF] developed a unified and automated methodology for the design of mechatronic systems. The proposed approach consists in combining Bond-Graphs (BG) and Genetic programming (GP). However, this work is oriented to the structural part in mechatronic systems and lacks the controller part that is always present in these systems. In order to achieve the best topology and parameters of a mechatronic system, Behbahani and De Silvia [START_REF] Behbahani | Niching genetic scheme with bond graphs for topology and parameter optimization of a mechatronic system[END_REF] proposed a new evolutionary algorithm. This algorithm makes it possible to explore efficiently the design space. In this work, the control part has been taken into account. Another methodology combining (BG) and (GP) is presented by Samarakoon et al. [START_REF] Samarakoon | Design evolution of engineering systems using bond graphs and genetic programming[END_REF]. This approach was implemented in an industrial fish cutting machine. Nevertheless, the nonlinearity of the machine has been overlooked in this work.

Finally, although the structure of mechatronic systems has been treated during these previous works, the control part has been overlooked in some of them. This absence can influence the reality of these systems since the control is an essential part in mechatronic products. In this context, several works to model mechatronic systems have been done using the object-oriented modeling language Modelica and Dymola software. This tool makes it possible to control mechatronic systems while modeling their structures in the same environment. Modelica has been used by Ferretti et al. [START_REF] Ferretti | Virtual prototyping of mechatronic systems[END_REF]to model and simulate the behavior of a mechatronic machining center. Another work was done by Mcharek et al. [START_REF] Mcharek | Modeling and multiobjective optimization of an Electronic Throttle in open-loop[END_REF]. In this work, a mechatronic system called Electronic Throttle Body (ETB) was modelled with Modelica language. The contribution of this work is represented by taking into account the integration of the fail-safe technique in the ETB model.

Methodologies Based on Collaborative Design and Knowledge Sharing

Over many decades, the complexity of products is in increase. This increase implies also an increase in the number of people needed to develop these systems. Actors need to interact and exchange crucial data to obtain an ideal product. This can be achieved by following collaborative design. Such a strategy enable companies to reduce costs, improve quality and reduce time to market while innovating their products. In this context, researchers developed several methodologies in order to exchange and share knowledge. First, methodologies based on Product Lifecycle Management (PLM) have been used to exchange knowledge between actors. Chen et al. [START_REF] Chen | A framework for integrated design of mechatronic systems[END_REF] developed an approach of dividing the mechatronic system into components and classifying them into two parts an electrical and a mechanical one. After this classification, the constraints between components are defined and a table summarizing the components and the connection between them is established. This methodology reduces the complexity of mechatronic systems but the classification of components can cause problems while developing these systems. Alexopoulos et al. [START_REF] Alexopoulos | A web-services oriented workflow management system for integrated digital production engineering[END_REF] proposed a methodology for integrating data mechatronic system into the PLM system. In this approach, Automation Markup Language (AutomationML) has been used as an intermediate format to exchange data between different tasks. This language facilitates the integration and exchange of data. However, there are not many translation tools available because AutomationML is a new format. Another interesting study is proposed by Abid et al. [START_REF] Abid | SysML approach for the integration of mechatronics system within PLM systems[END_REF] based on the use of behavioral and structural diagrams of system Modeling Language (SysML) to model structure of a complex product in a PLM system. The PLM in this case makes it possible to coordinate the exchange between designers. This methodology is validated with a 3D robot. All these methodologies have shown their efficiency in the exchange of data between actors. However, they do not allow an easy access to data embedded in different models. Second, in the field of simulation and modeling, the goal is to provide an integrated environment for engineers in order to achieve optimal multi-domain design. This can be satisfied using technology called Process Integration and Design Optimization (PIDO). There are several tools to apply this technology such as Model Center [START_REF]Phoenix Integration[END_REF] and Isight. In this context, Azzouzi et al. [START_REF] Azzouzi | Process integration and design optimization technologies for modelling improvement[END_REF] proposed a methodology aims at improving system modeling by applying PIDO technology. This approach has been validated by applying it to a hybrid vehicle using an optimization tool (Model Center). This tool is used to identify unknown parameters of a system. An improvement of this methodology can be made by minimizing the error between model results and experimental data. However, PIDO technology is used in the optimization of a single architecture not to share knowledges in order to choose the most appropriate architecture among several ones of a system. Finally, Knowledge Based-Engineering (KBE) tools are used in collaborative design since they allow knowledge capitalization and reuse. Several methodologies were done using KBE tools. Badin [START_REF] Badin | Ingénierie hautement productive et collaborative à base de connaissances métier: vers une méthodologie et un méta-modèle de gestion des connaissances en configurations[END_REF], during his PhD research work, developed a new methodology based on configuration management of crucial knowledge for product design and numerical simulation. This approach allows engineers to capitalize, reuse and maintain coherently knowledges between different domains. As an improvement of this methodology, decision support approaches to guide designers can be integrated. Mcharek et al. [12], also proposed a new methodology based on KBE concept. This approach is characterized by its ability to support designers during design phases as well as its ability to reuse design results. The contribution of Mcharek et al. [START_REF] Mcharek | Knowledge sharing for mechatronic systems design and optimization[END_REF] is illustrated by its application to the design of a mechatronic system called Electronic Throttle Body (ETB). We will consider the same system during the application of the methodology. Therefore, our contribution in this paper is illustrated by the development of a new methodology, which will help us not only to collaborate between different actors of one architecture of a mechatronic system but also to compare between several architectures and to choose the most appropriate one.

Proposed Methodology

The interaction between several companies working at the same time for the same project has become a standard for manufacturers and more practically in the aeronautical and automotive fields. These fields are known by their complexity because of the integration of different disciplines. Therefore, each actor must have access to the right information, at the right time and with the right information format. Mechatronic systems are one of complex systems, which are characterized by the combination of mechanical, computer systems and electronics aided by control strategies. These systems can be represented by various architectures or solutions, which make their design very tedious. The design process is composed of seven steps as shown in Fig. 1. In this context, our contribution aims at developing a new methodology that helps engineers to choose the best architecture of mechatronic system among several ones in the conceptual phase. This methodology is divided into three main phases. The first phase of our methodology starts with an idea of designing a system in order to satisfy one or more specific needs. Once the function of the system is defined, a search of different design solutions must be established. If the different solutions have been selected, we can create the models for each solution. These models have different forms such as 3D, fluid, procurement and dynamic models. During the second phase, all project actors will collaborate with each other. This phase represents the core of our methodology since at this stage the actors can share their work to detect the inconsistencies between them and reduce correction returns. This phase starts with the extraction of crucial data from different established models and from the requirements imposed on the system. Once important data are extracted, they can be shared in the collaboration support. Indeed, only the data characterizing the interdependence between disciplines need to be shared. The last step of the second phase consists in visualizing and managing conflicts either in one model or between several models. In the third phase, a comparison between the different solutions or architectures of the system must be done and the best solution can be chosen based on this comparison. Figure 2 illustrates the global process of our methodology.

Case Study

The proposed approach is used for the choice of the best architecture among three of an electronic Throttle Body (ETB). The ETB adjusts the airflow entering to the combustion engine. Following the methodology steps, we apply it to the ETB. This application starts with finding the different solutions.

For the ETB we choose three different architectures one with a DC motor, one with a brushless motor and another one with a stepper motor. After choosing these Fig. 1 Product development process [START_REF] Masoudi | A guidance system for selecting an appropriate ecodesign checklist in the early stages of product development[END_REF] architectures, the different models for each architecture must be established. For this application, three models will collaborate with each other: a dynamic model using the environment Dymola and the object-oriented language Modelica, a 3D model using CATIA and procurement model, which consists in finding the most appropriate motors conforming within simulation results. The ETB model with a DC motor in Dymola environment is illustrated in Fig. 3. The DC motor is replaced by the brushless motor in the second alternative and by the stepper motor in the third one. After creating these dynamic models with Dymola software, the performances of each motor such as torque, power, voltage, etc. must be visualized. In the procurement models, three motors have be chosen conforming within visualized performances. Indeed, a single 3D model must be created since the difference between the three alternatives concerns motors type not the geometry.

Fig. 2 The global process of the proposed methodology Fig. 3 ETB model with a DC motor in Dymola environment Once all models have been created, crucial data can be extracted and shared between actors. Sharing data consists in organizing it in Information Core Entities (ICEs), which are generic entities aims to organize and capitalize crucial data, as presented in Table 1. In the end of this application, the existing conflicts are visualized and managed and the best architecture is chosen.

Results and Discussion

Our methodology was tested in KARREN platform. Figure 4 summarizes the obtained results. Only the chosen DC motor in procurement model can be implemented in the 3D model. The other chosen motors have incompatible dimensions with the 3D model. Therefore, we choose the first alternative of the ETB with a DC motor.

Conclusions

This paper has presented a new methodology based on collaborative design and knowledge sharing to choose the best architecture among several ones of a mechatronic system. This methodology was applied to an Electronic Throttle Body. In this application, we considered three alternatives each one with a different type of motors (with a DC motor, brushless motor and stepper motor) and we created three models for each alternative which are 3D model, procurement model and dynamic model. In our future work, we will add another domain such as safety system to improve our methodology.

Fig. 4

 4 Fig. 4 Obtained results in KARREN platform

Table 1

 1 Organization of crucial data

	Parameters Unit Description	Constraints
	ICE_performance motor
	P	W	Power of motor -
	C	N.m Torque of motor <1.5 N.m
	ICE_dimension motor
	L	mm Motor length	-
	D	mm Motor diameter -

Acknowledgements. The authors would like to express their acknowledgments to the mechanical laboratory of Sousse (Tunisia) and QUARTZ laboratory in SUPMECA-Paris.