Hivert Hélène
email: helene.hivert@ens-lyon.fr

NUMERICAL SCHEMES FOR KINETIC EQUATION WITH DIFFUSION LIMIT AND ANOMALOUS TIME SCALE

Keywords: 2010 Mathematics Subject Classification. Primary: 35B40, 65L04, 41A60 ; Secondary: 65M22 BGK equation, anomalous diffusion limit, multiscale scheme, asymptotic preserving scheme, asymptotic analysis

In this work, we propose numerical schemes for linear kinetic equation, which are able to deal with a diffusion limit and an anomalous time scale of the form ε 2 (1 + |ln(ε)|). When the equilibrium distribution function is a heavy-tailed function, it is known that for an appropriate time scale, the mean-free-path limit leads either to diffusion or fractional diffusion equation, depending on the tail of the equilibrium. The bifurcation between these two limits is the classical diffusion limit with anomalous time scale treated in this work. Our aim is to develop numerical schemes which work for the different regimes, with no restriction on the numerical parameters. Indeed, the degeneracy ε → 0 makes the kinetic equation stiff. From a numerical point of view, it is necessary to construct schemes able to undertake this stiffness to avoid the increase of computational cost. In this case, it is crucial to capture numerically the effects of the large velocities of the heavy-tailed equilibrium. Since the degeneracy towards the diffusion limit is very slow, it is also essential to respect the asymptotic behavior of the solution, and not only the limit. Various numerical tests are performed to illustrate the efficiency of our methods in this context.

1. Introduction. The modeling of a large amount of particles, from a theoretical or numerical point of view, is a very active field of research. The direct application of Newton's laws leads to a large system of coupled equations, one for each particle of the system. Since such a huge number of unknowns is beyond the reach of numerical computations, the so-called microscopic scale is not adapted to the numerical analysis of these systems. Instead, an approach based on statistical physics is preferred, where the distribution function of particles f ε (t, x, v) depending on the time t ≥ 0, the position x ∈ R d , and the velocity v ∈ R d is considered. The parameter ε ∈ (0, 1] denotes here the Knudsen number, which is proportional to the mean free path of the particle.

Provided an initial condition f (0, x, v) = f in (x, v), we consider the following kinetic equation for f ε

Θ(ε)∂ t f ε + εv • ∇ x f ε = L(f ε ). (1)
The quantity Θ(ε) = ε 2 (1 + ln (ε)), is a suitable scaling parameter to be chosen according to the nature of L, in order to capture a non-trivial dynamic when ε goes to 0. The linear operator L describes the collisions of the particles. We will consider the particular case of the BGK operator

L(f ε ) = ρ ε M -f ε , (2)
with

ρ ε (t, x) = R d f ε (t, x, v)dv =: f ε (t, x).
Here, the equilibrium function M is even, positive and normalized to 1

M (-v) = M (v) > 0 for all v ∈ R d and M = 1.
In what follows, the integral in v of f on a domain D will always be denoted f D .

The integral is on the whole space R d when no domain is specified.

The asymptotic analysis of (1) when ε → 0 can be found in [START_REF] Mellet | Fractional diffusion limit for collisional kinetic equations[END_REF]. Our aim in this paper is to construct numerical schemes for (1) that fit with this asymptotic analysis. The properties of the equilibrium function M , and the scaling Θ(ε) of (1) make various asymptotic behaviors arise. For instance, considering the scaling Θ(ε) = ε with a local Maxwellian equilibrium in [START_REF] Bardos | Fluid dynamic limits of kinetic equations. i. formal derivations[END_REF] and the Boltzmann operator instead of the linear collision operator L, leads when ε goes to 0 to a fluid limit for [START_REF] Bardos | Fluid dynamic limits of kinetic equations. i. formal derivations[END_REF] (see [START_REF] Bardos | Fluid dynamic limits of kinetic equations. i. formal derivations[END_REF][START_REF] Bardos | Fluid dynamic limits of kinetic equations. ii: Convergence proofs for the boltzmann equations[END_REF]). With the BGK collision operator, the scaling Θ(ε) = ε 2 , and a global Maxwellian equilibrium M (v) = (2π) -d/2 e -v 2 /2 , the solution f ε of (1) degenerates when ε → 0 to a distribution at equilibrium f = ρM , where ρ solves a diffusion equation ∂ t ρ -∇ x • (D∇ x ρ) = 0, with initial condition ρ in = f in . A rigorous derivation of diffusion-type equations in this last case was first investigated in [START_REF] Wigner | Nuclear reactor theory[END_REF][START_REF] Bardos | Diffusion approximation and computation of the critical size[END_REF][START_REF] Bensoussan | Boundary layers and homogenization of transport processes[END_REF][START_REF] Degond | Diffusion limit for non homogeneous and non-microreversible processes[END_REF]. When M is a Maxwellian and Θ(ε) = ε 2 , the diffusion coefficient D is finite because of the exponential decay of M for large v, and is given by

D = v ⊗ vM . (3)
However, for non Maxwellian equilibrium, it may happen in many situations that the coefficient D becomes infinite. An important example is the case of a so-called heavy-tailed function M (v) ∼ 1/|v| d+α for large v, α ∈ (0, 2). In this case, the scaling Θ(ε) = ε 2 is not the suitable choice and does not lead to a non trivial dynamics when ε goes to 0. It is well-known that, in this case, the time scale Θ(ε) should be modified according to α in order to capture a non-trivial dynamics, which turns out to be a fractional diffusion model [START_REF] Mellet | Fractional diffusion limit for collisional kinetic equations[END_REF][START_REF] Ben Abdallah | Fractional diffusion limit for collisional kinetic equations: A Hilbert expansion approach[END_REF]. Heavy-tailed equilibria arise in the study of granular media (see [START_REF] Ernst | Scaling solutions of inelastic Boltzmann equations with over-populated high energy tails[END_REF][START_REF] Bobylev | Boltzmann equations for mixtures of Maxwell gases: Exact solutions and power like tails[END_REF][START_REF] Bobylev | On some properties of kinetic and hydrodynamic equations for inelastic interactions[END_REF]), astrophysical plasmas (see [START_REF] Mendis | Cosmic dusty plasma[END_REF][START_REF] Summers | The modified plasma dispersion function[END_REF]), tokamaks (see [START_REF] Del Castillo-Negrete | Non-diffusive transport in plasma turbulence: A fractional diffusion approach[END_REF]) or in economy (see [START_REF] Kleiber | Statistical size distributions in economics and actuarial sciences[END_REF]). Usually, the fractional diffusion equation describes the motion of the particles driven by a Levy process (see [START_REF] Moor | Limites diffusives pour des équations cinétiques stochastiques[END_REF][START_REF] Moor | Fractional diffusion limit for a stochastic kinetic equation[END_REF]), while the classical diffusion is governed by a Brownian process.

From a numerical point of view, the case α ∈ (0, 2) has been treated in [START_REF] Crouseilles | Numerical schemes for kinetic equations in the anomalous diffusion limit. Part I: The case of heavy-tailed equilibrium[END_REF][START_REF] Crouseilles | Multiscale numerical schemes for kinetic equations in the anomalous diffusion limit[END_REF][START_REF] Wang | An asymptotic-preserving scheme for linear kinetic equation with fractional diffusion limit[END_REF]. In this paper, we consider the so-called critical case α = 2, of the above heavy-tailed equilibrium, which induces a different asymptotic behavior when ε goes to 0. For a sake of simplicity we will consider here the following particular case

M (v) = m |v| < 1 m/|v| d+2 |v| ≥ 1, ( 4
)
where m is such that M = 1. With this assumption, it has been shown in [START_REF] Mellet | Fractional diffusion limit for collisional kinetic equations[END_REF] that if Θ(ε) ∼ ε 2 | ln(ε)| in [START_REF] Bardos | Fluid dynamic limits of kinetic equations. i. formal derivations[END_REF] the solution f ε of (1) converges weakly when ε → 0 to ρ 0 M , where ρ 0 solves a classical diffusion equation

∂ t ρ 0 -c d m∆ρ 0 = 0, ρ 0 (0, x) = ρ in (x) = f in (x). (5)
Here c d is a coefficient that only depends on the dimension d (c 1 = 2, c 2 = π, and c 3 = 8π/3). To have a non-vanishing of Θ(ε) for ε = 1, we shall consider

Θ(ε) = ε 2 (1 + |ln(ε)|) . (6)
In this work, we prove that when the initial data f in is at equilibrium, the convergence when ε → 0 of the solution f ε of (1) to the solution of ( 5) is very slow. However, we will show that f ε approaches much faster an intermediate equilibrium ρε M , where ρε is the solution of the following equation (in Fourier variable)

∂ t ρε + a ε (k) ρε = 0, t > 0, k ∈ R d ρε (0, k) = ρin (k), k ∈ R d , (7)
and

a ε (k) = 1 Θ(ε) ε 2 (k • v) 2 1 + ε 2 (k • v) 2 M , (8)
where ρε is the space Fourier transform of ρε .

The main goal of this work is to construct numerical schemes for (1) which do not need to be refined in order to capture the right asymptotic behavior, when ε goes to 0. There are, in fact, some stiffness in the problem [START_REF] Bardos | Fluid dynamic limits of kinetic equations. i. formal derivations[END_REF]. First, the smallness of ε imposes severe conditions on the numerical parameters, if a naive approach is used. We will follow an Asymptotic Preserving (AP) strategy [START_REF] Jin | Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations[END_REF][START_REF] Klar | An asymptotic-induced scheme for nonstationary transport equations in the diffusive limit[END_REF][START_REF] Klar | An asymptotic preserving numerical scheme for kinetic equations in the low Mach number limit[END_REF] to overcome this difficulty. Namely, if we consider a problem P ε which degenerates into a problem P 0 when ε tends to 0, then an AP scheme S h ε should enjoy the following properties: 1. at fixed ε, S h ε must be consistent with P ε when the discretization parameter h tends to 0. 2. S h ε must degenerate into a scheme S h 0 consistent with P 0 when ε goes to 0. An AP scheme can even enjoy the stronger property of being Uniformly Accurate (UA), meaning that its accuracy does not depend on ε.

The construction of AP schemes for the case of classical diffusion with classical time scale Θ(ε) = ε 2 has been widely investigated (see [START_REF] Jin | Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations[END_REF][START_REF] Jin | Discretization of the multiscale semiconductor Boltzmann equation by diffusive relaxation schemes[END_REF][START_REF] Jin | Discretization of the multiscale semiconductor Boltzmann equation by diffusive relaxation schemes[END_REF][START_REF] Jin | Uniformly accurate diffusive relaxation schemes for multiscale transport equations[END_REF][START_REF] Klar | An asymptotic-induced scheme for nonstationary transport equations in the diffusive limit[END_REF][START_REF] Buet | Asymptotic preserving scheme and numerical methods for radiative hydrodynamic models[END_REF][START_REF] Carrillo | Numerical schemes of diffusion asymptotics and moment closures for kinetic equations[END_REF][START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF][START_REF] Lemou | Micro-macro schemes for kinetic equations including boundary layers[END_REF]), but the heavy-tail of the equilibrium brings an additional stiffness, which is not usual in the classical cases. Indeed, the anomalous time scale [START_REF] Bobylev | On some properties of kinetic and hydrodynamic equations for inelastic interactions[END_REF] in ( 1) is designed to capture the effect of the high velocities in the asymptotic analysis. From a numerical point of view, taking into account these high velocities is then crucial to ensure the AP property of the scheme. In previous works (see [START_REF] Crouseilles | Multiscale numerical schemes for kinetic equations in the anomalous diffusion limit[END_REF][START_REF] Crouseilles | Numerical schemes for kinetic equations in the anomalous diffusion limit. Part I: The case of heavy-tailed equilibrium[END_REF][START_REF] Crouseilles | Numerical schemes for kinetic equations in the anomalous diffusion limit. Part II: Degenerate collision frequency[END_REF]), we investigated AP schemes for kinetic equations in the fractional diffusion limit. In the critical case, the same methodology can not be applied since the resulting schemes do not respect the dynamics of the convergence towards the diffusion equation [START_REF] Bardos | Fluid dynamic limits of kinetic equations. i. formal derivations[END_REF]. This is due to the slow convergence in this critical case, and a specific study will be performed to capture the right asymptotic behavior.

Three numerical schemes will be investigated in this paper. The first one is based on a fully implicit in time scheme for (1) using a Fourier variable in space. A suitable modification is introduced in this case, and its AP behavior is proved. Then, a scheme based on a micro-macro decomposition of the distribution function f ε is presented. It allows to avoid, if needed, the use of the Fourier transform and of the time implicit character of the previous scheme. Moreover, it can be adapted to deal with more general collision operators in [START_REF] Bardos | Fluid dynamic limits of kinetic equations. i. formal derivations[END_REF]. Once again, the convergence towards [START_REF] Bobylev | Boltzmann equations for mixtures of Maxwell gases: Exact solutions and power like tails[END_REF] has to be treated with care to ensure the good asymptotic behavior of the scheme. Eventually, we propose an approach based on an integral formulation of (1) in Fourier variable, which enjoys the stronger UA property.

To summarize, we will show that (1) approaches [START_REF] Bobylev | Boltzmann equations for mixtures of Maxwell gases: Exact solutions and power like tails[END_REF] with a rate ε 1 + |ln(ε)|, while ( 7) approaches the limit model ( 5) with the slower rate 1/ (1 + |ln(ε)|), and we will construct numerical schemes which respect this asymptotic behavior. This can be illustrated by the following diagram

P ε : (1) ------------→ ε √ 1+|ln(ε)| Pε : (7) ------------→ 1/(1+|ln(ε)|) P 0 : (5) h→0 -----→ -----→ h→0 -----→ h→0 S h ε ------------→ ε √ 1+|ln(ε)| Sh ε ------------→ 1/(1+|ln(ε)|) S h 0 (9)
The schemes S h ε and Sh ε are respectively consistent with the problems P ε and Pε when ε > 0 is fixed, and S h 0 is consistent with P 0 . The scheme S h ε approaches Sh

ε with rate ε 1 + |ln(ε)| while Sh ε approaches S 0 with rate 1/ (1 + |ln(ε)|)
when the discretization parameter h is fixed and when ε goes to 0. The asymptotic behavior P ε → Pε → P 0 will be proved in suitable functional spaces and illustrated by several numerical tests.

The paper is organized as follows. In the next section, we will start by establishing the convergence rates of (1) towards [START_REF] Bensoussan | Boundary layers and homogenization of transport processes[END_REF], and of (7) to [START_REF] Bensoussan | Boundary layers and homogenization of transport processes[END_REF]. Then, in Section 3, we present an asymptotic preserving scheme for [START_REF] Bobylev | Boltzmann equations for mixtures of Maxwell gases: Exact solutions and power like tails[END_REF], and the three asymptotic preserving schemes for [START_REF] Bardos | Fluid dynamic limits of kinetic equations. i. formal derivations[END_REF], which are tested in Section 4.

2. Degeneracy to the diffusion limit. In this section, we show that the convergence of (1) towards the diffusion limit (5) can be quantified by two steps as described by the diagram (9). This will be the basis of the numerical methods proposed in Section 3. In what follows, the space Fourier transform will be often used. We will denote f (t, k, v) (resp. ρ(t, k)) the space Fourier transform of f (t, x, v) (resp. ρ(t, x))

f (t, k, v) = R d e -ik•x f (t, x, v)dx.
We will use the following weighted

L 2 norm f (x, v) 2 L 2 M -1 = R d R d |f (x, v)| 2 1 M (v) dvdx, ( 10
)
where M is given by (4). We have the following propositions:

Proposition 1. Let f ε ∈ L ∞ 0, T ; L 2 M -1 be the solution of (1) with initial con- dition f in (x, v) = ρ in (x)M (v), with ρ in ∈ H N d R d (N 1 = 3, N 2 = N 3 = 6)
, where H N d R d denotes the usual Sobolev space, and let T > 0. Then, there exists a constant C such that

f ε -ρε M L ∞ (0,T ;L 2 M -1 ) ≤ CT ε 1 + |ln(ε)| ρ in H N d , (11)
where ρε ∈ L ∞ 0, T ; L 2 solves (7) with initial condition ρ in , and with M defined in (4).

Proposition 2. Assume that the previous assumptions hold, and let ρε ∈ L ∞ 0, T ; L 2 (R d ) be the solution of [START_REF] Bobylev | Boltzmann equations for mixtures of Maxwell gases: Exact solutions and power like tails[END_REF] with

ρ in ∈ H M d R d (M 1 = 2, M 2 = M 3 = 4)
, and let T > 0. Then, there exists a constant C such that

ρε -ρ 0 L ∞ (0,T ;L 2 (R d )) ≤ CT 1 + |ln(ε)| ρ in H M d , (12)
where ρ 0 ∈ L ∞ 0, T ; L 2 (R d ) solves (5), with initial condition ρ in .

These two propositions establish that the convergence ε → 0 of (1) makes appear two dynamics with very different rates in ε. Indeed, in a first step, the solution f ε approaches (relatively) quickly a distribution function at equilibrium ρε (t, x)M (v), with ρε solution of [START_REF] Bobylev | Boltzmann equations for mixtures of Maxwell gases: Exact solutions and power like tails[END_REF]. Then, in a second step, this density ρε goes much more slowly to the solution ρ 0 of the limit diffusion equation [START_REF] Bensoussan | Boundary layers and homogenization of transport processes[END_REF]. Hence, the approximation of f ε by its diffusion limit ( 5) is valid only for very small ε. At the numerical level, it implies that a large range of ε may stay out of reach if no appropriate strategy is used in the numerical schemes. The next parts present numerical methods for (1) degenerating when ε to 0 into numerical methods solving [START_REF] Bobylev | Boltzmann equations for mixtures of Maxwell gases: Exact solutions and power like tails[END_REF]. In particular, the schemes enjoy the AP property that is • when ε is fixed, the scheme is consistent with (1),

• the difference between the numerical scheme S h ε for (1) and the numerical scheme Sh ε for (7) converges to 0 when ε goes to 0. Afterwards, letting ε go to zero in the scheme Sh ε for (7) makes it degenerate in a scheme S h 0 solving the limit diffusion equation [START_REF] Bensoussan | Boundary layers and homogenization of transport processes[END_REF]. This approach enables to capture numerically the two scales of convergence of the solution of (1) to its diffusion limit, as described by [START_REF] Carrillo | Numerical schemes of diffusion asymptotics and moment closures for kinetic equations[END_REF].

The proofs of Prop. 1 and Prop. 2 follows the ideas developed in [START_REF] Ben Abdallah | Fractional diffusion limit for collisional kinetic equations: A Hilbert expansion approach[END_REF] in the case of a fractional diffusion limit for kinetic equation, coming from a heavy-tailed equilibrium function. We start by proving the following lemma, which gives the limit of a ε when ε goes to 0. Lemma 2.1. Considering a ε defined in [START_REF] Buet | Asymptotic preserving scheme and numerical methods for radiative hydrodynamic models[END_REF], there exists a constant C such that,

∀k ∈ R d 1. If d = 1, a ε (k) -2mk 2 ≤ C 1 + |ln(ε)| 1 + k 2 . ( 13
)
2. If d = 2, 3, |a ε (k) -c d m|k| 2 | ≤ C 1 + |ln(ε)| 1 + |k| 4 (14)
with c 2 = π and c 3 = 8π 3 . Proof of Lemma 2.1. In the 1-dimensional case, with the change of variables w = εv, the coefficient a ε reads

a ε (k) = m 1 + |ln(ε)| k 2 v 2 1 + ε 2 k 2 v 2 |v|<1 + m 1 + |ln(ε)| k 2 w 2 1 + k 2 w 2 1 |w| 3 |w|≥ε , (15)
where the last term in (15) can be computed for nonzero k with the change of variables u = 1/w

m 1 + |ln(ε)| k 2 w 2 1 + k 2 w 2 1 |w| 3 |w|≥ε = 2mk 2 + 2mk 2 1 + |ln(ε)| 1 2 ln ε 2 + 1 k 2 -1 ,
and the first term is bounded by Ck 2 /(1 + |ln(ε)|). Eventually the inequality ln(x)/2 ≤ 1 + x for x > 0, gives [START_REF] Moor | Fractional diffusion limit for a stochastic kinetic equation[END_REF].

When d = 2, 3, with the change of variables w = ε|k|v for nonzero k, (8) reads

a ε (k) = m 1 + |ln(ε)| (k • v) 2 1 + ε 2 (k • v) 2 |v|≤1 + m|k| 2 1 + |ln(ε)| (w • e) 2 1 + (w • e) 2 1 |w| d+2 |w|≥ε|k| , (16)
where e denotes any unitary vector. The first term is lower than C|k| 2 , it remains to consider the second one. Depending on the dimension, a change of variables in polar or spherical coordinates can be applied in it. Since no additional difficulty arise when d = 3, we treat here the case of polar coordinates, when d = 2. Choosing e such that w • e = |w| cos(θ), it can be rewritten as

(w • e) 2 1 + (w • e) 2 1 |w| d+2 |w|≥ε|k| = I 1 + I 2 + I 3 , with I 1 = 2π θ=0 1 r=ε|k| r 2 cos 2 (θ) 1 r 3 drdθ, (17)
I 2 = 2π θ=0 1 r=ε|k| r 2 cos 2 (θ) 1 1 + r 2 cos 2 (θ) -1 1 r 3 drdθ ( 18
)
I 3 = 2π θ=0 ∞ r=1 r 2 cos 2 (θ) 1 + r 2 cos 2 (θ) 1 r 3 drdθ. ( 19
)
At this stage, the last two integrals can be bounded

|I 2 | ≤ 2π θ=0 1 r=ε|k| r 4 cos 4 (θ) 1 + r 2 cos 2 (θ) 1 r 3 drdθ ≤ 2π θ=0 1 r=ε|k| r cos 4 (θ)drdθ ≤ C(1 + |k| 2 ) |I 3 | ≤ 2π 0 ∞ r=1 1 r 3 drdθ ≤ C,
and the first one can be computed explicitly

I 1 = -π ln(|k|) -π ln(ε).
Eventually, the inequality ln(|k|) ≤ |k| yields [START_REF] Moor | Limites diffusives pour des équations cinétiques stochastiques[END_REF].

Remark 1. The proof of Lemma 2.1 provides, for nonzero k, a rewriting of a ε in which the diffusion limit clearly appears. In dimension 1, it is

a d=1 ε (k) = 2mk 2 + m 1 + |ln(ε)| k 2 v 2 1 + ε 2 k 2 v 2 |v|<1 + 2mk 2 1 + |ln(ε)| 1 2 ln ε 2 + 1 k 2 -1 , (20)
and in dimension 2 it reads

a d=2 ε (k) = πm|k| 2 + m 1 + |ln(ε)| (k • v) 2 1 + ε 2 (k • v) 2 |v|≥1 - m|k| 2 1 + |ln(ε)| π (1 + ln(|k|)) - m|k| 2 1 + |ln(ε)| 2π θ=0 1 r=ε|k| r cos 4 (θ) 1 + r 2 cos 2 (θ) drdθ + m|k| 2 1 + |ln(ε)| 2π θ=0 ∞ r=1 cos 2 (θ) 1 + r 2 cos 2 (θ) 1 r drdθ. (21)
Eventually, a change of variables in spherical coordinates provides the case of dimension 3

a d=3 ε (k) = 8π 3 m|k| 2 + m 1 + |ln(ε)| (k • v) 2 1 + ε 2 (k • v) 2 |v|≥1 - m|k| 2 1 + |ln(ε)| 8π 3 (ln(|k|) + 1) - 2πm|k| 2 1 + |ln(ε)| π θ=0 1 r=ε|k| r sin 5 (θ) 1 + r 2 sin 2 (θ) drdθ + 2πm|k| 2 1 + |ln(ε)| π θ=0 ∞ r=1 sin 2 (θ) 1 + r 2 sin 2 (θ) 1 r drdθ. (22)
Proof of Prop. 2. We denote by ρ0 and ρε the solutions of ( 5) and ( 7) in space Fourier variable, with initial condition ρin , and a 0 (k) = c d mk 2 . They satisfy

ρε (t, k) = e -taε(k) ρin (k), ρ0 (t, k) = e -ta0(k) ρin (k).
Since a ε (k) ≥ 0 and a 0 (k) ≥ 0, the inequality |e

-x -e -y | ≤ |x -y| for x ≥ 0, y ≥ 0 yields        |ρ 0 -ρε | ≤ Ct 1 + |ln(ε)| 1 + k 2 |ρ in | if d = 1, |ρ 0 -ρε | ≤ Ct 1 + |ln(ε)| 1 + |k| 2 2 |ρ in | if d = 2, 3,
and gives [START_REF] Crouseilles | Numerical schemes for kinetic equations in the anomalous diffusion limit. Part II: Degenerate collision frequency[END_REF].

The convergence of the solution of ( 7) to the solution of ( 5) is very slow in ε, but Prop. 1 states that the degeneracy of ( 1) to ( 7) happens much faster. The proof of Prop. 1 is based on the proof of the similar result stated in [START_REF] Ben Abdallah | Fractional diffusion limit for collisional kinetic equations: A Hilbert expansion approach[END_REF] for the case of a kinetic equation with heavy-tailed equilibrium degenerating to a fractional diffusion equation. It uses the following properties of the collision operator, established in [START_REF] Degond | Diffusion limit for non homogeneous and non-microreversible processes[END_REF] and [START_REF] Mellet | Fractional diffusion limit for collisional kinetic equations[END_REF] in a more general case: Proposition 3. The collision operator L given by (2) has the following properties 1. The operator L is bounded in

L 2 M -1 (R d × R d ). 2. There exists a constant C > 0 such that ∀f ∈ L 2 M -1 , R d L(f ) f M dv ≤ - 1 2 R d |f -f M | 2 1 M dv. 3. For h ∈ L 2 M -1 , the equation L(g) = h has a solution if and only if h = 0. If g = 0, this solution is unique. Proof of Prop. 1.
As in [START_REF] Ben Abdallah | Fractional diffusion limit for collisional kinetic equations: A Hilbert expansion approach[END_REF], the proof relies on the following Hilbert expansion of the solution f ε of ( 1)

f ε = ρε M + g 1 + g 2 + r
, where ρ, g 1 , g 2 and r solve

(1 + εv • ∇ x ) g 1 = -εv • ∇ x ρε M (23) g 2 M -g 2 = -g 1 M + Θ(ε)∂ t ρε M (24)
∂ t r + ε Θ(ε) v • ∇ x r = 1 Θ(ε) ( r M -r) -∂ t g 1 -∂ t g 2 + ε Θ(ε) v • ∇ x g 2 (25)
with initial condition ρε (0,

x) = ρ in (x), g 1 (0, x, v) = g 2 (0, x, v) = r(0, x, v) = 0.
Considering [START_REF] Kleiber | Statistical size distributions in economics and actuarial sciences[END_REF] in Fourier variable leads to

ĝ1 = - iεk • v 1 + iεk • v ρε M. (26)
Furthermore, the left-hand side of ( 24) is such that g 2 M -g 2 = 0. Still in Fourier variable, ( 24) then reads -g 1 + Θ(ε)∂ t ρε = 0, that is using ( 26) and the symmetry of [START_REF] Buet | Asymptotic preserving scheme and numerical methods for radiative hydrodynamic models[END_REF]. Note that solving this equation leads to

M ∂ t ρε (t, k) + a ε (k) ρε (t, k) = 0, with a ε defined in
ρε (t, k) = e -aε(k)t ρin (k), ∂ t ρε (t, k) = -a ε (k)e -aε(k)t ρin (k). (27)
As ρ solves [START_REF] Bobylev | Boltzmann equations for mixtures of Maxwell gases: Exact solutions and power like tails[END_REF], the proof of the proposition relies on establishing estimates for

g 1 L ∞ (0,T ;L 2 M -1 ) , g 2 L ∞ (0,T ;L 2 M -1 ) and r L ∞ (0,T ;L 2 M -1 )
. Since a ε ≥ 0, ( 27) and the expression of ĝ1 in [START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF] give

R d |ĝ 1 | 2 dv M ≤ Θ(ε)a ε (k)|ρ in | 2 ,
with Θ(ε) defined in [START_REF] Bobylev | On some properties of kinetic and hydrodynamic equations for inelastic interactions[END_REF]. The estimation of a ε given by Lemma 2.1 yields

g 1 L ∞ (0,T ;L 2 M -1 ) ≤ Cε 1 + |ln(ε)| ρ in H N d , (28)
with N 1 = 1, and

N 2 = N 3 = 2.
The solution g 2 of ( 24) is unique if g 2 = 0. It can be rewritten as

-g 2 = -g 1 M + Θ(ε)∂ t ρε M,
and ( 26) and ( 7) imply g 2 = 0. The estimation of r can be done with (25) integrated on the characteristic curves, multiplied by r/M and integrated in velocity. Since g 2 = 0 and the term in r M -r gives a negative integral with the second point of Prop. 3, it reads

1 2 d dt R d r 2 (t, X(t), v) M dv ≤ - R d ∂ t g 1 (t, X(t), v) r(t, X(t), v) M dv,
where X(t) = x + ε Θ(ε) tv. After an integration in time and space, and using a Cauchy-Schwarz inequality, it eventually writes

1 2 r L ∞ (0,T ;L 2 M -1 ) ≤ ∂ t g 1 L 1 (0,T ;L 2 M -1 )
, so that an upper bound for the right-hand side is needed to complete the proof of the proposition. A differentiation in time of ( 26), together with [START_REF] Mellet | Fractional diffusion limit for collisional kinetic equations[END_REF], gives

∂ t g 1 2 L 2 M -1 ≤ Θ(ε) R d a ε (k) 3 e -aε(k)t |ρ in | 2 dk, hence, using Lemma 2.1 ∂ t g 1 L 1 (0,T ;L 2 M -1 ) ≤ CT ε 1 + |ln(ε)| ρ in H N d , which completes the proof.
3. Numerical schemes. This section presents three numerical schemes designed to approximate the solution of (1), and to respect the asymptotic behavior of Prop. 1 and Prop. 2 when ε becomes small. These schemes are based respectively on an implicit formulation of (1) in Fourier variable, a micro-macro decomposition, and an integral formulation of (1). These three strategies have been investigated in [START_REF] Crouseilles | Numerical schemes for kinetic equations in the anomalous diffusion limit. Part I: The case of heavy-tailed equilibrium[END_REF][START_REF] Crouseilles | Numerical schemes for kinetic equations in the anomalous diffusion limit. Part II: Degenerate collision frequency[END_REF] in the case of the fractional diffusion limit of the kinetic equation, α ∈ (0, 2). Here we extend them to the case α = 2. In this critical case, we required that not only the limit equation ( 5), but the asymptotic behavior (7) of the numerical solution is prescribed.

In what follows, we will denote by t n = n∆t, 0 ≤ n ≤ N , such that N ∆t = T the time discretization, and we will set f n ≈ f (t n ). All the numerical tests will be performed at time T = 0.1, the time step ∆t being specified in each case. The schemes require a velocity discretization, we will denote

f Nv,D = ∆v Nv j=1 v j ∈D f (v j ), (29)
the quadrature formula for the integration in velocity on the domain D. When no domain is specified, the integration is for all the index 1 ≤ j ≤ N v . In the tests, we will consider a uniform discretization of the domain |v| ≤ 10, with N v points symmetrically distributed on both sides of the origin to ensure that vM Nv = 0. Namely, we will consider N v = 2N v with N v = 100 and use the grid

v j = -10 + ∆v 2 + (j -1)∆v, 1 ≤ j ≤ 2N v , (30)
with ∆v = 10/N v . The normalization constant m of the equilibrium M is chosen such that M Nv = 1. The space domain is bounded, and we consider periodic conditions, allowing the use of the Fourier variable. At the discrete level, the discrete Fourier variable will be used. In the tests, we will consider x ∈ [-1, 1], discretized with N x points such that

x i = -1 + (i -1)∆x, 1 ≤ i ≤ N x , (31)
where ∆x = 2/N x . This space discretization is linked to the Fourier modes we use for the computation of the discrete Fourier transform, that are the integers k such that -N x /2 ≤ k ≤ N x /2. At the discrete level, the usual L 2 norm in space, and the norm L 2 M -1 defined in (10) will be computed using the space and velocity discretization

ρ 2 L 2 Nx = ∆x Nx i=1 |ρ(x i )| 2 , f 2 L 2 M -1 ,Nx ,Nv = ∆x∆v Nx i=1 Nv j=1 |f (x i , v j )| 2 M (v j ) . ( 32
)
Eventually, in the following numerical tests, we will take N x = 32 and consider the initial condition

f in (x, v) = (1 + sin(πx))M (v). ( 33
)
We want our schemes to enjoy a stronger property than the AP one. Indeed, they respect the two steps of the convergence of (1) towards ( 5) following the diagram [START_REF] Carrillo | Numerical schemes of diffusion asymptotics and moment closures for kinetic equations[END_REF].

The peculiar asymptotic behavior of (1) comes from the effect of the high velocities in the equilibrium function M . Indeed, its heavy-tailed character makes the usual diffusion coefficient (3) be infinite. The anomalous scaling Θ(ε) is chosen to counterbalance the weight of M for large velocities. It is then necessary to take them into account in the numerical computations. One idea would be to consider an adaptive grid for large velocities, to use as much large velocities as needed. However, linking the stiffness ε and the discretization parameters is not in the spirit of AP strategies. We propose a method based on analytical modifications of the terms degenerating into (8) in the schemes. These modifications are done consistently in the schemes, and ( 8) is carefully discretized to ensure the degeneracy to the limit equation [START_REF] Bensoussan | Boundary layers and homogenization of transport processes[END_REF]. Once it is done, we are able to prove that the schemes respect the asymptotic behavior of (1).

3.1. Discretization of a ε . Since the degeneracy of (1) to the diffusion limit [START_REF] Bensoussan | Boundary layers and homogenization of transport processes[END_REF] goes through the equation [START_REF] Bobylev | Boltzmann equations for mixtures of Maxwell gases: Exact solutions and power like tails[END_REF], it is necessary to handle the degeneracy of ( 7) to [START_REF] Bensoussan | Boundary layers and homogenization of transport processes[END_REF] numerically to design AP schemes for [START_REF] Bardos | Fluid dynamic limits of kinetic equations. i. formal derivations[END_REF]. In Fourier variable, both ( 5) and ( 7) are differential equations, which can be solved analytically, and that can also be solved with classical schemes. For instance, a semidiscrete-in-time implicit scheme for [START_REF] Bobylev | Boltzmann equations for mixtures of Maxwell gases: Exact solutions and power like tails[END_REF]

reads ρn+1 ε -ρn ε ∆t + a ε (k) ρn+1 ε = 0, ( 34
)
and it degenerates into a scheme solving (5) if a ε (k) → a 0 (k) = mc d |k| 2 when ε goes to 0. If the velocity integrations are done continuously, this is ensured thanks to [START_REF] Moor | Fractional diffusion limit for a stochastic kinetic equation[END_REF]. Unfortunately, the discrete velocity integration breaks this property. Indeed, with [START_REF] Summers | The modified plasma dispersion function[END_REF] the discrete version of a ε (k) writes

1 Θ(ε) ε 2 k 2 v 2 1 + ε 2 k 2 v 2 M Nv = ∆v 1 + |ln(ε)| Nv j=1 k 2 v 2 j 1 + ε 2 k 2 v 2 j M (v j ) (35) ∼ ε→0 ∆v 1 + |ln(ε)| Nv j=1 k 2 v 2 j M (v j ).
Since the quadrature uses a finite number of points, the sum above is finite whilst its continuous version is infinite. Then, the naive numerical version of a ε (k) tends to 0 when ε → 0, which is not satisfactory. Consequently, it is necessary to modify analytically a ε before applying the discretization [START_REF] Summers | The modified plasma dispersion function[END_REF], to ensure the correct asymptotic behavior of a ε . In Remark 1, the expressions a ε (k) rewritten in ( 20)-( 21)-( 22) make the limit behavior of a ε clearly appear. Hence, at the discrete level, we rather use a discrete version of these terms, where the integrals are computed with an usual quadrature [START_REF] Summers | The modified plasma dispersion function[END_REF]. It ensures the discrete version a Nv ε of a ε to degenerate to the expected diffusion coefficient when ε tends to 0. To do so, in dimension 1, we use the following discrete version of a ε , from (20)

a Nv ε (k) =              m 1 + |ln(ε)| k 2 v 2 1+ε 2 k 2 v 2
Nv,|v|≤1

+2m|k| 2 1 + 1 2 ln ε 2 + 1 |k| 2 -1 1+|ln(ε)| , k = 0 0, k = 0, (36)
which respects the estimates ( 13) at the discrete level. Eventually, with a Nv ε , the Euler scheme for [START_REF] Bobylev | Boltzmann equations for mixtures of Maxwell gases: Exact solutions and power like tails[END_REF] degenerates when ε goes to 0, to a Euler scheme for [START_REF] Bensoussan | Boundary layers and homogenization of transport processes[END_REF].

The numerical tests highlight this behavior. Indeed, Fig. 1 displays the solutions given by the scheme (34), with a ε computed with the discretizations (35) and (36). The solutions are computed with ∆t = 10 -2 , for a range of ε. We also plot the numerical solution of the limit model approximated by ρn+1 0 -ρn

0 ∆t + 2m|k| 2 ρn+1 0 = 0. ( 37
)
When a ε is computed directly with (35), we observe that the solution of (34) does not tend to the limit diffusion when ε goes to 0, but to the initial condition ρ in = f in Nv , with f in defined in (33). Conversely, the solution of (34) with a ε computed with (36), comes close to the solution of (37) when ε goes to 0. Moreover, Fig. 2 suggests that this convergence in ε happens with the speed 1/(1 + |ln(ε)|), proved in Prop. 2. To obtain this figure, we computed the solution ρ ε,∆t of (34) with (36), and the solution ρ ∆t of (37), with ∆t = 10 -2 . Then, we plotted the relative error

Error(ε) = ρ ε,∆t -ρ ∆t L 2 Nx ρ ∆t L 2 Nx , with • L 2
Nx defined in (32), as a function of ε ∈ (0, 1). Figure 1. For ∆t = 10 -2 , the solutions of (37) and (34) at time T = 0.1 for different values of ε. Left: when a ε is computed with (35). Right: when a ε is computed with (36).

Implicit scheme (IS).

The first idea to write a scheme for (1) which enjoys the AP property is to use a fully implicit scheme for (1) written in Fourier variable. Indeed, it is well-known that in the case of the classical diffusion limit, such a scheme degenerates when ε → 0 into a scheme solving the limit equation. In our case, such an approach does not work, due to the effects of large velocities. Then, it is necessary to suitably write the scheme to make this degeneracy appear. We For ∆t = 10 -2 , the relative error between the solution of the scheme (34) and the limit scheme (37) at time T = 0.1, in function of ε (log scale).

start with (1) written in space Fourier variable, and we consider a fully implicit time discretization

f n+1 -f n ∆t + 1 Θ(ε) (1 + iεk • v) f n+1 = 1 Θ(ε) ρn+1 M, (38)
where the dependence in ε of f n+1 , f n and ρ n+1 has been omitted to simplify the notations. This scheme degenerates when ε → 0 into

f n+1 = ρn+1 M.
Therefore, provided that the density ρn+1 is known, the scheme (38) enjoys the AP property, and f n+1 reads

f n+1 = Θ(ε) f n + ∆t ρn+1 M Θ(ε) + ∆t + iε∆t k • v . ( 39
)
An expression for ρn+1 is given by a fully implicit time discretization of (1) integrated in velocity

ρn+1 -ρn ∆t + 1 Θ(ε) iεk • v f n+1 = 0,
which writes, using the expression of f n+1 in (39)

ρn+1 -ρn ∆t + iεk • v f n Θ(ε) + ∆t + i∆tεk • v + ∆t Θ(ε) iεk • vM Θ(ε) + ∆t + iε∆tk • v ρn+1 = 0. ( 40
)
Using the symmetry of M , the last term of the previous expression can be rewritten as ∆t Θ(ε)

iεk • vM Θ(ε) + ∆t + iε∆tk • v ρn+1 = ∆t 2 Θ(ε) ε 2 (k • v) 2 (Θ(ε) + ∆t) 2 + ε 2 ∆t 2 (k • v) 2 M ρn+1 (41) = 1 Θ(ε) ε 2 (k • v) 2 Θ(ε)+∆t ∆t 2 + ε 2 (k • v) 2 M ρn+1 . ( 42
)
The expression (41) indicates that the term is of order ∆t 2 for fixed ε, and thus can be removed or modified with no incidence on the accuracy of the scheme. The expression (42) shows that it degenerates into the coefficient a ε (k) when ε decreases.

Since the second term of (40) vanishes for small ε, if the integrations in velocity are done continuously, the scheme degenerates when ε to 0 into a scheme solving [START_REF] Bobylev | Boltzmann equations for mixtures of Maxwell gases: Exact solutions and power like tails[END_REF].

To ensure the AP property at the fully discrete level, (42) can be approximated by

∆t Θ(ε) + ∆t 2 1 Θ(ε) ε 2 (k • v) 2 1 + ε 2 (k • v) 2 M ρn+1/2 = ∆t Θ(ε) + ∆t 2 a ε (k)ρ n+1/2 ,
where ρn+1/2 can be taken equal to ρn or ρn+1 , depending on the desired limit scheme (implicit or explicit). Then the coefficient a ε is numerically approximated as in Section 3.1 to ensure the AP property. Eventually, we have the following proposition Proposition 4. We consider the following scheme defined for all k and all time indices

0 ≤ n ≤ N, N ∆t = T by            f n+1 -f n ∆t + 1 Θ(ε) (1 + iεk • v) f n+1 = 1 Θ(ε) ρn+1 M ρn+1 -ρn ∆t + iεk • v f n Θ(ε) + ∆t + i∆tεk • v Nv + ∆t ∆t + Θ(ε) 2 a Nv ε (k)ρ n+1/2 = 0, (43)
with Θ(ε) and a Nv ε defined in (6) and (36), and with initial condition

f 0 (k, v) = fin (k, v), ρ0 (k) = f 0 (k, •)

Nv

. The quantity ρn+1/2 can be chosen equal to ρn or ρn+1 depending on the desired asymptotic scheme (implicit or explicit in time). This scheme has the following properties:

1. The scheme is of order 1 for fixed ε > 0 and preserves the total mass. 2. The scheme is AP: for a fixed ∆t, the scheme solves the diffusion equation (5) when ε goes to 0

ρn+1 -ρn ∆t + mc d |k| 2 ρn+1/2 = 0, (44)
with c 1 = 2, c 2 = π, c 3 = 8π 3 and where m has been defined in (4). Remark 2. The degeneracy of (43) to a scheme solving (5) respects the two dynamics discussed in the diagram (9). Indeed, the numerical tests show that (43) degenerates with speed

ε 1 + |ln(ε)| to ρn+1 -ρn ∆t + a Nv ε (k)ρ n+1/2 = 0, (45)
which solves [START_REF] Bobylev | Boltzmann equations for mixtures of Maxwell gases: Exact solutions and power like tails[END_REF]. Then letting ε going to 0 in the previous scheme makes it degenerate to (44) with speed 1/(1 + |ln(ε)|).

Proof. The first point of Prop. 4 is a direct consequence of the fact that we used an Euler scheme for [START_REF] Bardos | Fluid dynamic limits of kinetic equations. i. formal derivations[END_REF], and that the modifications we applied to it were done consistently. Taking k = 0 in the second line of (43) yields the conservation of the total mass. Eventually, with the modification of the term (42), the AP character is straightforward. Indeed, letting ε becomes small with fixed ∆t in (43) leads to

   f n+1 = ρn+1 M ρn+1 -ρn ∆t + a Nv ε (k)ρ n+1/2 = 0,
which degenerates into a scheme solving (5) when ε → 0.

Micro-macro scheme (MMS).

In the previous part, we proposed an implicit scheme for (1), which has the same behavior as the solution of (1) when ε goes to 0. However, the use of the Fourier variable may be restrictive in some cases, and its extension to more general collision operators seems difficult. Moreover, the implicit treatment of the transport operator induces high computational cost, especially in multi-dimensional cases. We propose here a scheme based on a micromacro decomposition of the distribution function f ε , which respects the behavior of f ε when ε → 0. It treats the transport terms explicitly, and can be extended to the case of general collision operators of [START_REF] Mellet | Fractional diffusion limit for collisional kinetic equations[END_REF] with a strategy similar to [START_REF] Lemou | Relaxed micro-macro schemes for kinetic equations[END_REF]. Denoting ρ ε = f ε , we consider the decomposition f ε (t, x, v) = ρ ε (t, x)M (v) + g ε (t, x, v), such that g ε = 0. Inserting it in (1) and integrating in velocity gives

∂ t ρ ε + ε Θ(ε) v • ∇ x g ε = 0. ( 46
)
Then, we multiply (46) by M and we subtract it from (1) to get an equation for g ε

∂ t g ε + ε Θ(ε) v • ∇ x ρ ε M + ε Θ(ε) (v • ∇ x g ε -v • ∇ x g ε M ) = - 1 Θ(ε) g ε . (47)
A semi-discrete-in-time numerical scheme can be designed following [START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF], in which we implicit the stiffest terms

       g n+1 -g n ∆t + ε Θ(ε) v • ∇ x ρ n M + ε Θ(ε) (v • ∇ x g n -v • ∇ x g n M ) = - 1 Θ(ε) g n+1 ρ n+1 -ρ n ∆t + ε Θ(ε) v • ∇ x g n+1 = 0.
(48) Here, the dependence in ε of ρ n and g n is omitted to simplify the notations. For fixed ε, it is a consistent scheme, it remains to check wether it preserves the correct asymptotic. With the first line of (48), g n+1 = O(ε) when ∆t is fixed, and then g n+1 = -εv • ∇ x ρ n M + o(ε). When reported in the second line of (48), this gives the following limit scheme for ρ

ρ n+1 -ρ n ∆t - 1 1 + |ln(ε)| v • ∇ x (v • ∇ x ρ n M ) = 0, (49)
which does not corresponds to the expected diffusion limit. Indeed, the term into brackets is infinite when the integration in velocity is done continuously, but it is finite when the discretization ( 29) is applied. Hence, in this latter case, when ε → 0, the limit scheme is ρ n+1 = ρ n , which is not consistent with [START_REF] Bensoussan | Boundary layers and homogenization of transport processes[END_REF]. It is then necessary to modify (48) to ensure the correct asymptotic behavior. To do so, we modify the macro equation by making appear the inverse of the transport operator. We first remark that v • ∇ x g n+1 = v • ∇ x f n+1 , we then use a semi discrete implicit formulation of f as in the previous part (see (38) for instance)

f n+1 = Θ(ε) Θ(ε) + ∆t I + ∆t Θ(ε) + ∆t εv • ∇ x -1 f n + ∆t Θ(ε) + ∆t I + ∆t Θ(ε) + ∆t εv • ∇ x -1 ρ n+1 M.
This expression is injected in the flux of the macro equation of (48) to get

ρ n+1 -ρ n ∆t + ε Θ(ε) + ∆t v • ∇ x I + ∆t Θ(ε) + ∆t εv • ∇ x -1 f n (50) + ε Θ(ε) ∆t Θ(ε) + ∆t v • ∇ x I + ∆t Θ(ε) + ∆t εv • ∇ x -1 ρ n+1 M = 0. ( 51
)
Since the scheme is first order in time, the bracket in (50) can be simplified consistently in v • ∇ x f n to avoid the inversion of a differential operator. For the same reason, we could remove (51) with no incidence on the consistency of the scheme; however, it degenerates when ε → 0 to the diffusion of ( 5). Therefore, this term is kept to ensure the AP property of the scheme, but it must be treated with care to take correctly the effects of the high velocities into account. In Fourier variable, and using the evenness of M , (51) reads

1 Θ(ε) ∆t Θ(ε) + ∆t 2 ε 2 (k • v) 2 1 + ∆t Θ(ε)+∆t 2 ε 2 (k • v) 2 M ρn+1 .
This term is of order ∆t 2 when ε is fixed, and it degenerates to a ε (k)ρ n+1 (with a ε defined in ( 8)) when ε becomes small with ∆t fixed. Thus, (51) can be modified consistently in ∆t 2 /(Θ(ε) + ∆t) 2 a ε (k)ρ n+1/2 , which ensures both the AP character of the scheme, and the fact that the numerical solution solves [START_REF] Bobylev | Boltzmann equations for mixtures of Maxwell gases: Exact solutions and power like tails[END_REF] when ε belongs to the intermediate regime in which the solution f ε of (1) solves [START_REF] Bobylev | Boltzmann equations for mixtures of Maxwell gases: Exact solutions and power like tails[END_REF]. Here, ρn+1/2 can be chosen equal to ρn or ρn+1 .

The passage from the semi-discrete-in-time scheme to the fully discretized scheme can then be done with the quadrature [START_REF] Summers | The modified plasma dispersion function[END_REF], and the suitable discretization of a ε defined in [START_REF] Buet | Asymptotic preserving scheme and numerical methods for radiative hydrodynamic models[END_REF], as in section 3.1. The derivatives in space are treated with a classical first order upwind scheme; for a sake of simplicity in the notations, we still write here the spatial derivatives as continuous. Eventually, we have the following proposition. Proposition 5. We consider the following scheme (discretized in time and velocity) defined for all x ∈ R d , v j , 1 ≤ j ≤ N v and all time indices 0 ≤ n ≤ N, N ∆t = T > 0 by

         g n+1 -g n ∆t + ε Θ(ε) v • ∇ x ρ n M + ε Θ(ε) v • ∇ x g n -v • ∇ x g n Nv M = - g n+1 Θ(ε) ρ n+1 -ρ n ∆t + ε Θ(ε) + ∆t v • ∇ x g n Nv + ∆t Θ(ε) + ∆t 2 F -1 a Nv ε ρn+1/2 = 0, ( 52
)
with Θ(ε) and a Nv ε defined in [START_REF] Bobylev | On some properties of kinetic and hydrodynamic equations for inelastic interactions[END_REF] and in (36), where F -1 denotes the inverse of the Fourier transform, and with initial condition ρ 0 = f in Nv , and g 0 = f in -ρ 0 M . The quantity ρn+1/2 can be chosen equal to ρn or ρn+1 depending on the desired asymptotic scheme (implicit or explicit in time). This scheme has the following properties:

1. The scheme is of first order in time for any fixed ε > 0 and preserves the total mass. 2. The scheme is AP: for a fixed ∆t, the scheme solves the diffusion equation [START_REF] Bensoussan | Boundary layers and homogenization of transport processes[END_REF] when ε goes to 0

ρn+1 -ρn ∆t + mc d |k| 2 ρn+1/2 = 0, ( 53
)
with c 1 = 2, c 2 = π, c 3 = 8π 3
, and where m has been defined in (4). Remark 3. As in the case of the implicit scheme, this scheme enjoys a property stronger than being AP, since it respects the diagram [START_REF] Carrillo | Numerical schemes of diffusion asymptotics and moment closures for kinetic equations[END_REF]. Hence, it respects the dynamics of the convergence to the solution of the diffusion equation ( 5) when ε → 0. Indeed, the numerical tests show that (52) degenerates with speed ε 1 + |ln(ε)| to the scheme (45), which solves [START_REF] Bobylev | Boltzmann equations for mixtures of Maxwell gases: Exact solutions and power like tails[END_REF]. The degeneracy to the limit scheme (53) happens much more slowly, with speed 1/(1 + |ln(ε)|).

Remark 4. If we want to avoid the use of Fourier variable, we can replace

F -1 a Nv ε ρn+1/2 by -mc d ∆ x ρ n+1/2
. This modification is consistent with (1), and still enjoys the AP property, but the convergence to the limit scheme does not make a scheme for [START_REF] Bobylev | Boltzmann equations for mixtures of Maxwell gases: Exact solutions and power like tails[END_REF] arise in the intermediate regime in ε.

Proof. Since the first point comes immediately from the derivation of the scheme, let us prove the second point. The first line of (52) gives that g n+1 = o(1) when ε → 0. Once injected in the second line, it yields

ρ n+1 -ρ n ∆t + F -1 a Nv ε ρn+1/2 = o(1),
which is a consistent approximation of [START_REF] Bobylev | Boltzmann equations for mixtures of Maxwell gases: Exact solutions and power like tails[END_REF]. Finally, a Nv ε tends to mc d |k| 2 when ε → 0, which concludes the proof.

Integral formulation based scheme (DS).

In the previous parts, we proposed two numerical schemes solving [START_REF] Bardos | Fluid dynamic limits of kinetic equations. i. formal derivations[END_REF], which in addition respect the asymptotic behavior of the continuous solution when ε goes to 0. They enjoy the AP property, but their precision is not uniform in ε. In this section, we propose a scheme based on an integral formulation of (1) in Fourier variable. This approach was investigated in [START_REF] Crouseilles | Numerical schemes for kinetic equations in the anomalous diffusion limit. Part I: The case of heavy-tailed equilibrium[END_REF]- [START_REF] Crouseilles | Numerical schemes for kinetic equations in the anomalous diffusion limit. Part II: Degenerate collision frequency[END_REF] in the case of the fractional diffusion limit. This scheme was shown to be uniformly accurate (UA) : its precision does not depend on the stiffness parameter ε. A similar strategy is performed here, to deal with the diffusion case with an anomalous scaling.

We start from the Duhamel form for fε , solution of (1) in space Fourier variable

fε (t, k, v) = A 0 (t, k, v) + t Θ(ε) 0 e -s(1+iεk•v) ρε (t -Θ(ε)s, k)M (v)ds, where A 0 (t, k, v) = e -t(1+iεk•v) fin (k, v). ( 54
)
Evaluating at time t n+1 and integrating in velocity, we get

ρε (t n+1 , k) = A 0 (t n+1 , k, •) + n j=0 t j+1 Θ(ε) t j Θ(ε) e -s(1+iεk•v) M ρε (t n+1 -Θ(ε)s, k)ds.
(55) To write a scheme on ρε , we perform a quadrature of order 2 in the integral in s;

we have ∀j ∈ [[0, n]], ∀s ∈ tj Θ(ε) , tj+1 Θ(ε) , ρε (t n+1 -Θ(ε)s, k) = 1 - Θ(ε)s -t j ∆t ρε (t n+1 -t j ) + Θ(ε)s -t j ∆t ρε (t n+1 -t j+1 ) + O(∆t 2 ). ( 56
)
Assuming that the time derivatives of ρε are uniformly bounded in ε, the remainder is uniformly bounded by a term of magnitude O ∆t 2 independently of ε. Inserting (56) in (55) yields the following scheme for ρε

ρn+1 = A 0 (t n+1 , k, •) + n j=0 b j ρn-j + c j ρn+1-j , (57)
where the dependence in ε of ρε has been omitted to simplify the notations. The coefficients b j , c j , 0 ≤ j ≤ n are given by

b j = t j+1 Θ(ε) t j Θ(ε) Θ(ε)s -t j ∆t e -s(1+iεk•v) M ds c j = t j+1 Θ(ε) t j Θ(ε) 1 - Θ(ε)s -t j ∆t e -s(1+iεk•v) M ds.
The consistency of this scheme comes directly from its derivation, but the discretization of the integrals in velocity appearing in b j and c j is crucial to ensure the AP property. Indeed, as in the previous parts, their continuous limit when ε goes to 0 makes the diffusion (5) appear. But the effects of the high velocities are not captured if ( 29) is used. Therefore, it is necessary to modify analytically b j and c j before applying the velocity discretization. The velocity integration of b j and c j can be rewritten as

e -s(1+iεk•v) M = e -s + e -s(1+iεk•v) -e -s M |v|<1 + e -s(1+iεk•v) -e -s M |v|≥1 , because M = 1.
After that, the change of variables w = εv is performed in the last term to capture the high velocity effects, and the integrals in s of b j and c j are computed exactly to get

b j = -e -t j+1 Θ(ε) + Θ(ε) ∆t e -t j Θ(ε) 1 -e -∆t Θ(ε) + m e -t j+1 Θ(ε) - e -t j+1 Θ(ε) (1+iεk•v) 1 + iεk • v |v|<1 (58) + Θ(ε) ∆t m e -t j Θ(ε) (1+iεk•v) 1 -e -∆t Θ(ε) (1+iεk•v) (1 + iεk • v) 2 -e -t j Θ(ε) 1 -e -∆t Θ(ε) |v|<1 + ε 2 m 1 |w| d+2   e -t j+1 Θ(ε) - e -t j+1 Θ(ε) (1+ik•w) 1 + ik • w   |w|≥ε + Θ(ε) ∆t ε 2 m 1 |w| d+2 e t j Θ(ε) (1+ik•w) 1 -e -∆t Θ(ε) (1+ik•w) (1 + ik • w) 2 |w|≥ε - Θ(ε) ∆t ε 2 m 1 |w| d+2 e -t j Θ(ε) (1 -e -∆t Θ(ε) ) |w|≥ε , c j = e -t j Θ(ε) - Θ(ε) ∆t e -t j Θ(ε) 1 -e -∆t Θ(ε) (59)
- Θ(ε) ∆t m e -t j Θ(ε) (1+iεk•v) 1 -e -∆t Θ(ε) (1+iεk•v) (1 + iεk • v) 2 -e -t j Θ(ε) 1 -e -∆t Θ(ε) |v|<1 + m e -t j Θ(ε) (1+iεk•v) 1 + iεk • v -e -t j Θ(ε) |v|<1 + ε 2 m 1 |w| d+2   e -t j Θ(ε) (1+ik•w) 1 + ik • w -e -t j Θ(ε)   |w|≥ε (60)
- Θ(ε) ∆t ε 2 m 1 |w| d+2 e -t j Θ(ε) (1+ik•w) 1 -e -∆t Θ(ε) (1+ik•w) (1 + ik • w) 2 |w|≥ε + Θ(ε) ∆t ε 2 m 1 |w| d+2 e -t j Θ(ε) 1 -e -∆t Θ(ε) |w|≥ε .
However, it is still not sufficient to ensure the AP character of the scheme when the integrals are discretized since, when j = 0, the term (60) reads

Θ(ε) m 1 + |ln(ε)| (k • v) 2 1 + ε 2 (k • v) 2 |v|<1 + m 1 + |ln(ε)| 1 |w| d+2 (k • w) 2 1 + (k • w) 2 |w|≥ε ,
which is equal to Θ(ε)a ε , with a ε defined in [START_REF] Buet | Asymptotic preserving scheme and numerical methods for radiative hydrodynamic models[END_REF] and rewritten in [START_REF] Degond | Diffusion limit for non homogeneous and non-microreversible processes[END_REF] (note that the equality ( 15) is valid for any dimension). In Section 3.1, we established that the discretization of this term is crucial to ensure the AP property of our schemes. Following the strategy presented previously, we then replace (60) by an appropriate discretization Θ(ε)a Nv ε defined in (36) when j = 0. Eventually, c 0 rewrites

c 0 = Θ(ε)a Nv ε (k) + 1 (61) - Θ(ε) ∆t 1 -e -∆t Θ(ε) - Θ(ε) ∆t m 1 -e -∆t Θ(ε) (1+iεk•v) (1 + iεk • v) 2 -1 + e -∆t Θ(ε) |v|<1 - Θ(ε) ∆t mε 2 1 |w| d+2 1 -e -∆t Θ(ε) (1+ik•w) (1 + ik • w) 2 -1 + e -∆t Θ(ε) |w|≥ε .
The expressions b j (j ≥ 0), c j (j ≥ 1) and c 0 , defined in ( 58)-( 59)-( 61), computed with the discretization (29) make the scheme (57) AP. Moreover, let us remark that the scheme only uses ρ and not f . It implies that the computation of the whole distribution f is not necessary, which represents a gain of computational cost, since the problem for f is of higher dimension than the problem on ρ. Let us note that f can even so be computed in Fourier variable with a scheme similar to (57). Indeed, Duhamel formula of (1) between the times t n and t n+1 leads to

fε (t n+1 , k, v) = e -∆t Θ(ε) (1+iεk•v) fε (t n , k, v) + ∆t Θ(ε) 0 e -s(1+iεk•v) ρε (t n+1 -Θ(ε)s, k)M (v)ds. (62)
As in the scheme for ρ, we apply the quadrature (56) in the integral in s, and the scheme for f reads

f n+1 = e -∆t Θ(ε) (1+iεk•v) f n + β ρn M + γ ρn+1 M,
where the dependence in ε of f has been omitted, and with

β = ∆t Θ(ε) 0 Θ(ε)s ∆t e -s(1+iεk•v) ds = - e -∆t Θ(ε) (1+iεk•v) 1 + iεk • v + Θ(ε) ∆t 1 -e -∆t Θ(ε) (1+iεk•v) (1 + iεk • v) 2 (63)
γ = ∆t Θ(ε) 0 1 - Θ(ε)s ∆t e -s(1+iεk•v) ds = 1 1 + iεk • v - Θ(ε) ∆t 1 -e -∆t Θ(ε) (1+iεk•v) (1 + iεk • v) 2 . (64)
Eventually, we have the following proposition.

Proposition 6. Provided the initial condition f 0 (k, v j ) j∈[[1,Nv]] = fin (k, v j ) j∈[[1,Nv]] , ρ0 (k) = f 0 (k, •)
Nv , we consider the following scheme, defined for all k and all time indices 0 ≤ n ≤ N, N ∆t = T , by

         ρn+1 (k) = A 0 (t n+1 , k, v) Nv + b 0 ρn (k) + n j=1 b j ρn-j (k) + c j ρn+1-j (k) 1 -c 0 f n+1 = e -∆t Θ(ε) (1+iεk•v) f n+1 + β ρn M + γ ρn+1 M, (65)
with A 0 , β, γ defined in (54)-( 63)-(64), and b j (j ≥ 0), c j (j ≥ 1) and c 0 defined in (58)-( 59)-(61) and computed with [START_REF] Summers | The modified plasma dispersion function[END_REF]. This scheme has the following properties:

1. The scheme is first order in time and preserves the total mass.

2. The scheme is AP: for a fixed ∆t, it solves the diffusion equation (5) when ε goes to zero ρn+1 -ρn ∆t + mc d |k| 2 ρn+1 = 0, with c 1 = 2, c 2 = π, c 3 = 3π 8 , and where m has been defined in (4). 3. Moreover, the semi-discrete-in-time scheme enjoys the UA property: it is first order uniformly in ε

∃C > 0, sup ε∈(0,1] f N (•, •) -f (T, •, •) L ∞ k,v ≤ C∆t.
Remark 5. The UA property of the scheme ensures that it respects the dynamics of the convergence towards the diffusion equation [START_REF] Bensoussan | Boundary layers and homogenization of transport processes[END_REF].

Proof. The conservation of the mass can be proved by induction on the expression of ρn+1 for k = 0. Moreover, since the proof of the third point of Prop. 6 is very similar to the one presented in [START_REF] Crouseilles | Numerical schemes for kinetic equations in the anomalous diffusion limit. Part I: The case of heavy-tailed equilibrium[END_REF]- [START_REF] Crouseilles | Numerical schemes for kinetic equations in the anomalous diffusion limit. Part II: Degenerate collision frequency[END_REF] in the case of the fractional diffusion limit, we just prove the AP character of the scheme (65). First of all, let us remark that A 0 ,b j and c j (j ≥ 1) defined in ( 54)-( 58)-( 59) are exponentially small when ε goes to 0. The expressions (58) and (61) then give the following equivalents for b 0 and

c 0 b 0 = Θ(ε) ∆t + Θ(ε)o (1) , c 0 = Θ(ε)a Nv ε + Θ(ε) ∆t + Θ(ε)o(1).
Hence the scheme for ρ degenerates when ε goes to 0 into ρn+1 -ρn ∆t + a Nv ε ρn+1 = 0.

The asymptotic behavior of a Nv ε has already been studied, it remains to consider the limit of the expression of f n+1 when ε goes to 0. From (63) and (64), we have

β = O(Θ(ε)), γ = 1 + O(Θ(ε)),
and eventually the scheme for f n+1 degenerates for small ε into f n+1 = ρn+1 M.

4. Numerical tests. In this section, we propose numerical tests in dimension 1 to validate the schemes of the previous parts. For simplicity, we will denote IS the implicit scheme of Prop. 4, MMS the micro-macro scheme of Prop. 5, and DS the integral formulation based scheme of Prop. 6. The implicit Euler scheme (45) for the equation [START_REF] Bobylev | Boltzmann equations for mixtures of Maxwell gases: Exact solutions and power like tails[END_REF], with a ε computed with (36), will be denoted quasi-diff. Eventually, the scheme (44) for the limit diffusion equation ( 5) will be denoted diff.

To highlight the consistency of the schemes, we will compare their solutions to the solution of a reference scheme which is an explicit Euler scheme in Fourier variable for (1)

         f n+1 -f n ∆t + 1 Θ(ε) (1 + iεkv) f n = 1 Θ(ε) ρn M ρn+1 = f n+1 Nv f 0 (k, v) = fin (k, v). (66)
As detailed in the beginning of Section 3, all the tests are performed at time T = 0.1, the time step being specified in each case. We will consider N x points in space on the interval [-1, 1] with periodic boundary conditions, and the velocities are discretized with N v = 199 points symmetrically distributed on the interval [-10, 10]. The initial condition is (33). In the DS scheme, the velocity integrations with the variable w in (58)-( 59)-( 61) are computed with the same grid as the integrals with the v variable.

4.1. Implicit scheme (IS). In this section, we test the properties of the IS scheme. First of all, its consistency is tested. The left part of Fig. 3 shows that, for ε = 1, its solution coincides with the solution of the explicit scheme (66). For this figure, the two solutions are computed with the time step ∆t = 10 -2 . Then, the solution ρ ∆t of the IS scheme for a range of ∆t is compared to the solution ρ ∆t ref =10 -6 of the IS scheme for ∆t ref = 10 -6 at time T = 0.1. The relative error between these densities

Error consistency (∆t, ε) = ρ ref erence -ρ ∆t L 2 Nx ρ ref erence L 2 Nx , (67)
with

• L 2
Nx defined in (32), is displayed in the right-hand side of Fig. 3 in function of ∆t. As it is a line with slope 0.98, the numerical order of the method is 1. Then, the AP character of the IS scheme is tested. The dynamics of the convergence towards the diffusion limit ( 5) is highlighted in the left part of Fig. 4, where the densities obtained with the IS scheme are presented for a range of ε, and compared to the densities given by quasi-diff and diff schemes. These solutions are computed with ∆t = 10 -2 . For intermediate ε, the solution of the IS scheme sticks to the solution of quasi-diff, and the two densities goes together to the solution of diff when ε tends to 0. The convergence towards the solution of quasi-diff happens with speed ε 1 + |ln(ε)|, as suggested by the error study presented in the right-hand side of Fig. 4. Denoting by f the distribution function obtained with the IS scheme, and ρ quasi-dif f the density obtained with the quasi-diff scheme for ∆t = 10 -2 , the quantity

Error AP (ε) = f -ρ quasi-dif f M L 2 M -1 ,Nx ,Nv , (68)
with 4.2. Micro-macro scheme (MMS). In this section, we focus on the MMS scheme and we highlight its properties with the same tests as in the previous section. As the transport operator is treated with a classical upwind scheme in the MMS scheme, the CFL condition imposes to take time steps smaller than in the previous section. The left-hand side of Fig. 5 shows that the solutions of the MMS and explicit schemes with ∆t = 10 -4 and ε = 1 coincide. The numerical order of the MMS scheme is studied in the right-hand side of Fig. 5, where the consistency error (67) for ε = 1 is displayed in function of ∆t. For this figure, the density ρ ref erence is the density obtained with the MMS scheme for ∆t = 10 -6 and ρ ∆t are densities obtained with the MMS scheme for a range of ∆t. As expected, the numerical order of the method is 1. Concerning the AP character of the MMS scheme, the left-hand side of Fig. 6 shows, for a range of ε, the densities obtained with the MMS, quasi-diff and diff schemes with ∆t = 10 -4 . Once again, the solution of the MMS scheme first joins the solution of the quasi-diff scheme when ε becomes small and the two reach the solution of diff scheme together when ε tends to 0. The speed of the convergence towards the solution of quasi-diff scheme is ε 1 + |ln(ε)| stated in Prop. 1. Indeed, the right-hand side of Fig. 6 displays the error (68), where ρ quasi-dif f , and the solution f given by the MMS scheme are computed with ∆t = 10 -4 .

• L 2 M -1 ,Nx ,Nv

Integral formulation based scheme (DS).

In this section, we test the integral formulation scheme DS. For ε = 1, we check in the left-hand side of Fig. 7, that the solutions of the explicit and DS schemes are close. This figure was obtained with ∆t = 10 -2 . The order of accuracy of the DS scheme is studied in the right-hand side of Fig. 7. For ε = 1, it displays the error (67) between ρ ref erence defined as the solution of the DS scheme for ∆t = 5 • 10 -5 , and the solution of the DS scheme ρ ∆t for a range of ∆t. In a logarithmic scale, the slope of the line obtained is close to 1, as expected. Then, we can study the AP character of the DS scheme. The dynamics of the convergence towards the diffusion limit ( 5) is highlighted in the left-hand side of Fig. 8. It presents the densities obtained with the DS, quasi-diff and diff schemes for ∆t = 10 -2 . The solution of the DS scheme has the right behavior, since it is very close to the solution of quasi-diff for intermediate ε, and remains stuck to it as it goes to the solution of diff when ε goes to 0. The right-hand side of Fig. 8 displays the error study of the convergence in ε of the solution of the DS scheme to the solution of the quasi-diff scheme. The error (68) is displayed in a logarithmic scale in this figure, where both ρ quasi-dif f and the solution f of the DS scheme are computed with ∆t = 10 -2 . The convergence happens with the expected ε 1 + |ln(ε)| rate.

To highlight the fact that the DS scheme is of order 1 uniformly in ε, we compare the results given by the DS scheme for ∆t ref = 5 • 10 -5 to the results given by the same DS scheme for different values of ∆t and ε. With these densities, the error (67) is displayed in Fig. 9 as a function of ε. We observe that the error curves are stratified with respect to ε, showing the uniformity of the scheme with respect to ε.

Conclusion.

In this paper, we presented three AP schemes for a kinetic equation with a diffusion limit and an anomalous time scale We first proved that the convergence of the solution of the kinetic equation towards the solution of the asymptotic model happens with a slow rate in ε, but that this convergence goes through an intermediate model with a quicker convergence rate. The degeneracy ε → 0 makes the problem stiff, which complicates the numerical computations, but moreover the slow rate of the convergence towards the asymptotic model makes the limit model unattainable in the numerical computations. Hence, we proposed an adaptation of the AP approach for numerical schemes, by requiring in addition the schemes to respect the dynamics of the degeneracy towards the limit equation. A second stiffness arise in the design of numerical schemes for kinetic equations with the anomalous diffusion limit we studied, since the high velocities play a crucial role in the asymptotic analysis of the model. From a numerical point of view, it is then necessary to take these high velocities correctly into account to ensure the AP property of the schemes. To deal with this stiffness, we proposed a method based on analytical transformation of the terms leading to the asymptotic model in the schemes. A similar approach was proposed in [START_REF] Crouseilles | Numerical schemes for kinetic equations in the anomalous diffusion limit. Part I: The case of heavy-tailed equilibrium[END_REF][START_REF] Crouseilles | Numerical schemes for kinetic equations in the anomalous diffusion limit. Part II: Degenerate collision frequency[END_REF] in the case of a fractional diffusion limit for the kinetic equation.

The first scheme we presented is based on a fully implicit discretization of the kinetic equation in Fourier variable, and the second one uses a micro-macro decomposition of the distribution function and treats the transport explicitely. Both of them enjoy the AP property and respect the dynamics of the convergence towards the limit model. The last scheme we proposed, which is based on an integral formulation of the kinetic equation in space Fourier variable, even enjoys the UA property.

In a near future, we aim at extending this work to more general context, considering an integral collision operator. In this context, the study can not be simplified by the use of the space Fourier variable. However, schemes based on a micro-macro decomposition of the distribution can be written. Moreover, the case of fractional diffusion limit of kinetic equations with an integral collision operator can also be treated with such an approach.

Figure 2 .

 2 Figure2. For ∆t = 10 -2 , the relative error between the solution of the scheme (34) and the limit scheme (37) at time T = 0.1, in function of ε (log scale).

Figure 3 .

 3 Figure 3. Left: for ∆t = 10 -2 and ε = 1, the solutions of the IS scheme and of the explicit scheme (66). Right: The relative consistency error (67) for the IS scheme (log scale).

Figure 4 .

 4 Figure 4. Left: for ∆t = 10 -2 , the solutions of IS, quasi-diff and diff schemes for different values of ε. Right: for ∆t = 10 -2 , the error between the solution of IS and quasi-diff scheme in function of ε (log scale).

15 Figure 5 .

 155 Figure 5. Left: for ∆t = 10 -4 and ε = 1, the solutions of the MMS scheme and of the explicit scheme (66). Right: The relative consistency error (67) for the MMS scheme (log scale).

Figure 6 .

 6 Figure 6. Left: for ∆t = 10 -4 , the solutions of MMS, quasi-diff and diff schemes for different values of ε. Right: for ∆t = 10 -4 , the error (68) between the solution of MMS and quasi-diff scheme in function of ε (log scale).

Figure 7 .

 7 Figure 7. Left: for ∆t = 10 -2 and ε = 1, the solutions of the DS scheme and of the explicit scheme (66). Right: The relative consistency error (67) for the DS scheme (log scale).

Figure 8 .

 8 Figure 8. Left: for ∆t = 10 -2 , the solutions of DS, quasi-diff and diff schemes for different values of ε. Right: for ∆t = 10 -2 , the error (68) between the solution of DS and quasi-diff scheme in function of ε (log scale).

Figure 9 .

 9 Figure 9. The error (67) as a function of ε. The density ρ ref erence is the density given by the DS scheme for ∆t ref = 5 • 10 -5 , and ρ ∆t are the densities given by the DS scheme for different values of ∆t (log scale).

H ÉL ÈNE HIVERT *

The author is supported by ERC starting grant MESOPROBIO.


A Le V An

B Peseux'

BOUNDARY ELEMENT ANALYSIS OF AN INTEGRAL EQUATION FOR THREE-DMENSIONAL CRACK PROBLEMS

In another paper, the authors proposed an integral equation for arbitrary shaped three-dimensional cracks. In the present paper, a discretization of this equation using a tensor formalism is formulated. This approach has the advantage of providing the displacement discontinuity vector in the local basis which varies as a function of the point of the crack surface. This also facilitates the computation of the stress intensity factors along the crack edge. Numerical examples reported for a circular crack and a semi-elliptical surface crack in a cylindrical bar show that one can obtain good resuhs, using few Gaussian poilnts and no singular elements.

INTRODUCTION

The integral equation method has been thoroughly experienced in structure analysis. Integral equations derived from the Somigliana representation were successfully applied to threedimensional problems, such as in References 1 and 2. In the latter reference, a particular symmetrical problem of cracked bodies was also investigated. Integral equations with kernels containing singular solutions, usually referred to as Kupradze elastic potentials, [START_REF] Kupradze | Three-Dimensional Problems of the Mathematical Theory of Elasticity and Thermoelasticity[END_REF] are particularly well suited to crack analysis. The use of these potentials in linear fracture mechanics is advantageous as it permits one to obtain directly stress intensity factors from the computed vector density. Embedded cracks in an infinite medium were studied in References 4-7. In particular, Bui [START_REF] Bui | [END_REF] showed that for plane cracks the mode I is entirely uncoupled from modes II and III. Later on, integral equations for three-dimensional cracks were proposed; 8 • 9 in the last reference use was made of a Kupradze double layer potential, and the unknown was the vector density directly related to the displacement discontinuity through the crack surface. In this paper, we carry out the discretization of the integral equation proposed in Reference 9, with a view to studying imbedded or surface crack problems.

For numerical purposes, the main difficulty is the reckoning of singular integrals defined in the sense of the principal value. Cruse 1 succeeded in giving the closed form of the principal value integral, using planar boundary elements. A more general way to evaluate two-dimensional singular integrals was investigated by K.azantzakis and Theocaris, 10 who also reduced these to one-dimensional finite-part integrals. 11 Among the most recent works, one can find References 12 and 13 where exact expressions for some specific two-dimensional finite!-part integrals are derived.

As for Kupradze potentials, Bui [START_REF] Bui | [END_REF] took a small polygon instead of a circle to perform principal value integrals, it seems, however, that this procedure cannot be used with ease. Lastly, one can find in Reference 14 an improvement of the evaluatjon of integrals of order 1 /r, 1 fr 2 , 1 fr 3 , by means of a minimum number of Gaussian points. As the considered distance r, small though it may be, remains always finite, these integrals are actually not singular ones. In any case, difficulties must certainly be expected when dealing with principal value integrals defined on curved surfaces.

In our analysis, we limit consideration to classical Gaussian quadrature, with no actual principal value computations. Nevertheless, two numerical examples given in the last section show that one can obtain fair results, using few Gaussian points, no singular integral computations and in particular, no singular elements around the crack edge. First, we deal with the embedded pennyshaped crack. The numerical results obtained proved to fit well to the analytical solution of Sneddon. 15 Then we treat the problem of a semi-elliptical crack in a cylindrical bar under combined loads. In the opening mode, the numerical results can be compared to the existing results 16 -18 which involved analogous crack geometries.

THEORETICAL BACKGROUND

Let us consider the problem of an elastic solid ~ containing an arbitrary shaped, imbedded or surface crack S, with the mixed boundary conditions:

y 0 eSuD.: lim t(x, Dx) = tg(y 0 , Dyo) !i)\S ~ x-. y~, n .. = n 70 (Ia) or: y 0 eSuD: lim u(x) = u 8 (Yo) (lb) !i)\S 3x -+ y~
where the superscript g stands for given, and the double sign + is related to the surface orientation, locally defined by the normal nyo at each point y 0 .

Displacement field expression

Let the displacement field be expressed by means of the Kupradze double-layer potential of the first kind

3 u(x)= LTT(y-x, n,)<p(y)d,S (2)
where the density cp is a vector function defined in the crack S, and T the Kupradze tensor The displacement discontinuity on the crack is obtained from

[u(y 0 )] =u(yri)-u(yi))='P<yo) ( 3
)
Using (3), one can calculate the stress intensity factors which are now directly related to the density cp.

Integral equations

It was shown in Reference 9 that, under the assumptions SeC

A3(u, v)~ y =F(u, v)ES cJ»(u, v)=q>(F(u, v))=<p(y)
t: the unit vector tangent to the crack edge aS, oriented according to the orientation of S. One should notice that the symbol pv before the two-dimensional integral relates to the principal value on S (not on A), and that the expression still holds in the general case when F.u is not

::lo---- . . e, i
Figure 1. General crack configuration with respectively a particular and an ordinary point, Yo and y perpendicular to F. v• The kernel of the line integral has been expressed for brevity as the product of a linear operator~ and the density q». Such line integrals appear when surface crack problems are involved.

In the next section, we shall concern ourselves with the discretization of the two-dimensional integral. As previously mentioned, this discretization will not include the study of the principal value.

DISCRETIZATION OF THE INTEGRAL EQUATION

Shape functions

The co-ordinates transformation is given by -r: e --+(u, v) = ( (N (e)){ u }e, (N(e) ){ v} e) (5) where e is a point of the reference element and the superscript e relates to an ordinary element. In this paper, we are interested to the most commonly encountered co-ordinate systems: the Cartesian, the polar and the cylindrical ones. In Cartesian co-ordinates (u, v)=(y 1 , y 2 ), in polar co-ordinates (u, v) = (p, 8) (Figure 2) and in cylindrical co-ordinates, with p constant, (u, v) = (8, y 3 ) (Figure 3).

Interpolation of~

In practice, it is advantageous to obtain the components of the density~ in the local basis which varies with the pointy of S. For this reason, we attach much importance to the choice of the coordinates system and we will consider different ways to interpolate the components of q, in the local basis. Use will be made of isoparametric elements.

In Cartesian co-ordinates, for instance, the local basis b is invariant; it coincides with the global, --fixed basis e of the space E 3 : The density q, is considered in the fixed basis:

~=cfJ1e1 +cfJ2e2+t/>3e3
Each of the components of ~ is then interpolated classically by

cPi = (N (~)){ cPi} e, iE { 1, 2, 3}
In cylindrical co-ordinates, the local basis varies with the point y:

The density will be considered in the variable basis:

cl » = 4> PeP+ l/Joe6 + ¢3 e3
and its components are now interpolated by

cP i = ( N (e)) { </J i} e, iE{p, e, 3} Relations (6) in turn yield Cn b { cfJ( e) } = b [ .Af ( ~) J • C 11 { cP} e (6a) (6b) (7)
where c, denotes the canonical basis of IR". The letters b, c" added in small characters recall the bases under consideration. In the case of four-node elements, the matrix [ JV] can be explicited as fol1ows:

Nl 0 0 l N 2 0 0 I N3 0 0 I N4 0 0 Cl2 o 1 o o I o o I b[JV (~)] = 0 Nl Nz N3 0 N4 0 I I I 0 0 N 1 I 0 0 Nz: 0 0 N3 I 0 0 N4 I I
and the density t1l in the case of cylindrical co-ordinates can be expressed by ( tP )e = ((t/Jp, t/Jo, tP3)\ (t/Jp, t/Jo, ¢3) 2 , (t/Jp, t/Jo, ¢3) 3 , (t/Jp, tPo, tP3) 4 )e

One can then express the derivatives of t1l with respect to the parameters u and v:

where

.d.. = .4{ . .d.. e 'JI,u u 'JI .d.. = vi( . .d.. e '+', v v 'f' if u =I(} if u = 0' if v =I(} if v = (}
The matrical representation of .# 1 in the bases b and c 12 (case of four-node elements) is

0 -N~ ol 0 -N2 ol 0 -NJ ol 0 -N4 0 I I I ctz I I I b[Ai' 1 (~)] = Nt 0 0 1 N 2 0 0 I N3 0 01 N 4 0 0 0 0 o I o 0 o I o 0 ol 0 0 0 I l I

Discretization of the integral equation

(8a) (8b) (9)
We study next how to transform each of the terms of the right-hand side of (4). First, let us consider the terms containing «Jl. u• We have (cp,u, e, F,Jny 0 =[cp,u•(er 1\ F,v)]ny 0 Similarly As for the term «11. u 1\ e, we have

= [(Aucpe)•(er 1\ F.J]nyo (using (8a)) = [ny 0 ® A!J (er 1\ F.v)J«Jle (lOa) (lOb) (11)
where an implicit sum is implied on the repeated subscript iE{l, 2, 3} or {p, (}, 3}, b 0 =(b 0 J is the basis related to y 0 . It should be noted that the components of the given stress t(y 0 , ny 0 ) in the lefthand side of(4) are obtained in the basis b 0 . Furthermore, in interacting or surface crack problems, b 0 may differ or not from b, according to the relative position of the particular point y 0 and the ordinary one y. Now in virtue of (11):

(F,v •ny 0 )«Jl,u 1\ er = (F, v •ny 0 )[b 0 i ® AJ(er 1\ b 0 J]cpe (1 Oc)
In the same manner, the fourh terms of ( 4) can be expressed by

It should be noted that, in (12d) for instance, as the given stress in the left-hand side is expressed in the basis b 0 , the vector e,. between braces { } (column-matrix, see Appendix) must be expressed in the basis b 0 , whereas the same e,. between angle brackets ( ) (row-matrix) must be expressed in the basis b.

The terms containing «P. v are discretized in the same manner; they yield relations analogous to (12) where «P.v takes the place of <l»,u, and F,u of F,v• Thus, for each particular point y 0 of the crack, we have discretized the integral equation using a tensor formalism, equations (10). This equation being discretized as a vectorial equation, and not as three separate scalar equations, this formalism allows us a compact programmation for numerical purposes. One has to solve an algebraic system in the form

[A] {cf>}={B}
where [A] is a full, non-symmetrical matrix, {4>} the unknown vector and {B} contains the components of the given stresses. In interacting or surface crack problems, { cf>} contains different vector densities corresponding respectively to different cracks or boundaries of the solid !?}.

N.B.: From the dimensional point of view, the matrices appearing in equations (12) (terms between brackets [ ]) are expressed in the inverse of a length or are dimensionless according as Cartesian, polar or cylindrical co-ordinates are involved.

NUMERICAL PARTS

I. Embedded circular crack

First, the program has been tested on the simple case of a circular crack centred at 0, with radius a, under a uniform pressure p, the analytical solution of which is well known: 15 8

(1-v 2 )• pa -------~ cp(p)= xE Jt-{p/a) 2 •e 3 , p= IIOyll
As for any plane crack imbedded in the infinite medium, the three modes are uncoupled. The stress intensity factor in the opening mode is given by K

=2p~

The natural choice of parameters for the penny-shaped crack is (u, v)=(p, 0). One should notice that this choice yields a co-ordinates transformation which is well suited to the crack geometry: the crack edge is especially better shaped than with the usual grid pattern, as shown in Figure 4.

Moreover, this choice of parameters facilitates the computation of the stress intensity factors which require to be taken as a limit when moving along a direction perpendicular to the crack edge. The given stress vector is expressed at each point Yo in a local basis b 0 =(ep 0 ; e 80 ; e 3 ), whereas the vector density 4> at a pointy in b=(ep; e 0 ; e 3 ); this situation is shown in Figure 5.

Figure 6 represents a quarter of a mesh with 181 interior nodal points (7 nodes along the radius and 36 along the polar angle (J, each sector being equal to 1 0°). The numerical integration was performed with 2 x 2 Gaussian points. As expected, the numerical values of the radial and tangential displacement discontinuities q, P and ¢ 9 are zero; moreover, all the nodes of the same radius have identical ¢ 3 values. It is interesting to note that, in mixed modes problems, the formulation herein proposed permits us to obtain directly the radial and tangential components of cf» (l/Jp and c/> 9 ) without additional calculations, because cf» is expressed in the local polar basis. Table I and Figure 7 give respectively the values of c/> 3 normalized by the theoretical maximum value (c/> 3 max)theor=8(1-v 2 )pa/nE, and its curve compared to Sneddon's solution. It should be noticed that one obtains in fact

1 0•95393 0•8 0•6 0•49889 0•35901 0 0•91244 0•87363 0•72500 0•56512 0•47032 0•35041 0 E/(16n(l-v 2
))•c/> 3 directly from the program, so that the Young's modulus E datum is not necessary. The relative error on </> 3 is maximum at the crack centre and equals 8•75 per cent; that of the computed stress intensity factor K 1 is 3•97 per cent. Figure 8 represents a finer mesh with 529 interior nodes (13 along the radius and 48 along the polar angel 8, each sector being equal to 7•5°). Use is always made of 2 x 2 Gaussian points. Table II and Figure 7 give the values of ¢ 3 normalized by the theoretical maximum ¢rvalue. The relative error on ¢ 3 at the centre is 2•97 per cent; that of the computed stress intensity factor K 1 is 0• 25 per cent.

Semi-elliptical crack in a cylindrical bar under combined loads

Let us consider a cylindrical bar with radius R, containing a semi-elliptical crack with semi-axes a and b. The crack occurs perpendicularly in the lateral surface of the cylinder, along a line which will be referred to as the surface-line of the crack (Figure 9). Three types of loads will be applied uniquely on the crack, the external surface of the bar being left free: a uniform pressure, a linear pressure along the y 1 -axis, and a shear distribution corresponding to the torsion of the bar (Figure 1 0).

According to the discussion in Reference 9, the surface crack problem can be investigated by considering a two crack problem: the first crack is the studied one, the second is a cylindershaped crack high enough to simulate a free bar. We introduce the notations: respectively the unknown density functions over Scr and Scyl

• cl-cr ='Per ° F cr' cl»cyl = 'Pcyl ° F cyl
The edge of Scyt is the union of the upper and the lower crack edges of the cylinder, plus a vanishing, counterclockwise contour surrounding L. The set of equations of the problem can be written as follows:

t(y 0 , ny 0 )=vpf (kernel containing q,cr)dudv .lcr + f (kernel containing q,cyl )du dv

where the kernels in (l3a) and (13b), given by (4), are not made explicit because of their length; S~~f and S~~~ denote respectively the upper and lower part of the cylindrical crack, separated by the section containing Scr• The relation (13e) holds if we choose the normal to Scr upward with respect to the crack and the normal to Scyl outward with respect to the cylinder. With another choice of surfaces orientations, however, this relation still holds but the signs before each term may change. Furthermore, one can easily verify that, in various surface crack problems, one or more relations of type (13e) may be obtained. Relations of type (l3e) allow one to verify that, though every contour integral is different from zero, there always appears a set of thesewhich cancel one another.

For numerical purposes, a crack depth-radius ratio b/R=0•4 and a semi-axes ratio b/a= 5/1 ~0•714 are considered. The cylinder lateral surface is divided into 748 elements with 731 interior nodes, while the semi-elliptical crack contains but 58 elements and 60 interior nodes (Figure 11). So large a node number on the free surface is due to the use of four-node elements with their sides along the parallel of the cylinder or along its axis. It should be wiser to use curved cylindrical elements (i.e. elements with distinct cylindrical nodal values) so as to realize a relatively refined mesh around the surface-line together with a coarse mesh far away, and to reduce notably the node number on the free surface. We have, however, not chosen such elements because of the difficulties caused to automatic mesh generations. In fact, the sole aim of this paper is to show the feasibility and the validity of the proposed approach. Anyway, the interesting feature is the little node number on the crack itself. Like the previous numerical example, use is made of four-node elements and the numerical quadrature is performed with 2 x 2 Gaussian points. Tables III, IV and V give the 4>rvalues, ie{l, 2, 3}, resulting from three loading types. From Tables III andIV, one can note that, under uniform or linear pressure loadings, the crack is seen to be in the opening mode. Table V shows on the contrary that, under torsion loadings, the crack deforms in a mixed mode II+ III. With the same maximum stress values (i.e. u=umax• see notation, Figure 10), ¢ 3 -values due to uniform pressures are found to be greater than those due to linear pressure.

It is important to note that one obtains Cartesian components of q, over the crack, and cylindrical components on the cylinder. Some interesting results can be deduced from the <Pivalues, ie { p, e, 3}, over the free surface. In the opening mode for instance, the point belongs to the upper crack face and located at the centre of the ellipse has, besides a positive ¢ 3 , a negative ¢P. This means that it moves upward with respect to the crack and outward with respect to the cylinder, as to some extent predicted by physical considerations.

Figures 12, 13 and 14 represent the crack in non-deformed and deformed states, where for a better representation deformations are emphasized by a multiplier equal to 5.

It is rather attractive to represent the stress intensity factor variations along the crack edge, but in surface crack problems, the stress intensity factor definition is almost difficult to establish. Table III. Components of E/(16n(l-v 2 ))•cPcr resulting from a uniform pressure: ¢ 1 =¢ 2 =0, ¢ 3 #:0; •cl»cr resulting from a linear pressure with umax = 1: c/J 1 =c/> 2 =0, c/J3 #:0 Generally speaking, the elastic state varies from plane strain in the interior to plane stress at the free surface. In the following, we shall limit our consideration to the stress intensity factors at the deepest point of the crack. There exists at this point a plane strain state; in the case of v = 0 One can find in Reference 1 7 a list of stress intensity factors for various crack and solid geometries. The so-called surface crack in solid cylinder is studied under tensile loading and pure bending (Figure 15). These actually are almond-shaped cracks which occur perpendicularly to the lateral surface of the cylinder. For lack of matter for comparisons, we consider, however, the crack of this type with b/ R =0•4; at the deepest point of the crack, the computed values are expected to be close to the above ones. The stress intensity factor is then expressed in the additive form 1 7 which yields in tensile loading:

(v =0•3, b/R =0•4, bja = 5/7, b= 5, u= 1, the unit of a, b, R is [L], that of u is [F• L-2 ], that of the above components is [F•L -1 ]) E E E E Node 16n(l-v 2 ) •¢ 3 Node 16n(l-v 2 ) •¢ 3 Node 16tt(l-v 2 )•¢ 3 Node 16n(l-v 2 )•c/> 3
E E E E Node 16n(1-v 2 ). c/> 3 Node 16n(1-v 2 ) •c/> 3 Node 16tt(l-v 2 ). c/> 3 Node 16n(l-v 2 ) •c/> 3 1 0•3994 2 0•3927 3 0•3927 4 0•3722 5 0•
E t6n(l-v 2 ) •c/> 1 E 161t(l-v 2 ) •l/J 2 Node E 16n(l-v 2 ) •l/> 1 E 16n(1 -v 2 ) •l/> 2
K, = [a• Fo +amu' F1 JJ'i{h F 0 = G[0•752 + 1•286/J + 0•37 Y 3 ] F 1 = G[0•923 + 0•199 Y 4 ] G = 0•92(2/n)• 1/cos fJ• [tan fJ/fJ] 1 1 2 Y = 1-sin/3 fJ = n:/2 • b/2R K.f(a Jib)= 0•800
and in pure bending:

(

K 1 /( (1 max Jib)= 0•606.
Recently, in Reference 18, use is made of the finite element method to study semi-elliptical cracks in cylindrical bars under tensile loadings. When computing the stress intensity factor K 1 , the assumption that a plane strain state holds in the neighbourhood of the crack is made in order to apply the virtual crack extension method. Charts for K 1 at the centre of the crack are given for b/a=O•Ol, 0•2, 0•5, 1•0, 2•0 and for various bfR ratios. An extrapolation from these K 1 -charts then gives us approximately the interval containing the K 1 /((1J'i{h) ratio: for b/ R =0•4, bfa= 5/7, with a uniform pressure (1, K.f((f J'i{h)e[0•93, 1•00].

The discrepancy between the above results is mainly due to the crack geometry differences; nevertheless, the values obtained are of the same order. For precise comparisons to be made, investigations are being pursued to obtain more complete abacuses of stress intensity factors, reported for various crack configurations. CONCLUDING REMARKS Equation [START_REF] Bui | [END_REF] has been discretized and written in a tensor form, equation (10). Whereas lengthy expressions should be expected if [START_REF] Bui | [END_REF] had been dealt with as three separate scalar equations, this approach yields a concise form of the discretized equation and thus allows a compact way of programming. Moreover, a crack parametrization being chosen; this enables one to obtain directly the components ofthe displacement discontinuity vector at each point, as well as the stress intensity factors along the crack edge, in a local basis. Unlike the usual finite element method, no singular element is required near the crack edge. The method herein proposed can also easily be applied to non-symmetrical crack configuration problems for which the boundary integral equation method becomes inefficient. Moreover, since no a priori assumptions of a plain strain or plain stress are to be made, the integral equation method can be used to investigate the stress singularity in the neighbourhood of the crack.

The example on the circular crack shows that one can obtain good results, using only four-node elements and 2 x 2 Gaussian points. The problem of the semi-elliptical crack in a cylindrical bar was also investigated with combined loadings. The numerical results obtained in opening mode were found to be in the range of different existing results for analogous crack geometries. The proposed method also allows one to obtain results for more complex loadings, in particular for torsion loading to which few works in the literature are devoted. APPENDIX Notation b=(b 1 ; b 2 ; b 3 ) a local basis at a pointy in E 3 c,. the canonical basis of ~~~ f} the elastic body characterized by constants (E, v) of} the boundary of q) e=(e 1 ; e 2 ; e 3 ) the fixed basis of E 3 E 3 the three-dimensional Euclidean space F a parametrization of S, defined on a domain A in IR 2 : A3(u, v)-+ y = F(u, v)ES. In general, F,u may not be perpendicular to F.v I unit tensor ny normal to S at yeS S the crack surface as the crack edge t(y, n,) stress vector at yeS, with respect to normal.n, T Kupradze tensor y, y 0 respectively an ordinary, a particular point of S pvJ., surface integral on S, defined in the principal value sense {a} column-matrix of components of a vector a in a basis. One may encounter b{ a} which emphasizes the basis bin question (a) row-matrix, = {a} T, transpose of {a}

2p, 6 ( 1

 61 + p,)e,-ny T(y -X, Dy) 8n(A.+ 2Jl)•r2 [ny®er-er®ny-(e,:ny)I]-8n(A.+ 2p,)•r2 er®er 1 is the unit tensor, r= l l y-x !l, e,.=(y-x)/ll y-xl!.

  Vy0 ES, t(y 0 , O,.ol = 16 7<(~ _ v 2 ) pv 1 (2(<1>. •• e, F. ,)n,0 -(1-2 v)H+ .•• e, n,.o) F., +(F,v•ny 0 )<P,u A e,)+3(e,-<P,u)((F,v, ny 0 , e,)e,+(n 10 •e,)e, A F,v) -2(cl»,v, e,, F,u)Dyo + (1-2 v )( (clJ,v, e,, Dyo) F,u + (F,u •nyo)cl»,v A e,) -3(e, •<P,v)((F,u, Dy 0 , e,)e,+(ny 0 •e,)e, 1\ F,u)) dudv -16n(~-v2) f .. 'C(n,o, y-Yo• ~)op(y)d,l (4) with now (Figure l): r=lly-y 0 ll,e,=(y-yo)/IIY-Yoll F: a parametrization of S, defined on a domain A in IR 2 :

Figure 2 .Figure 3 .

 23 Figure 2. Interpolation with polar co-ordinates: (a) local co-ordinates ~ =(~ 1 , ~2 ); (b) global coordinates (p, 0)

  (e,:cl-,u)[(F,v, Dy 0 , e,.)e,.+(ny 0 •e,.)e,. 1\ F,vJ =[((F,v, ny 0 , e,)e,+(ny 0 •e,.)e,. 1\ F,v)®.AJ e,]cpe Eventually, relations (10) are expressed in matrical form by (lOd) (cp,u, e,, F,v) Dyo = [ { Dyo} ( (e, A F,v) [.A uJ)J { cp }e (12a) (<J»,u, e,., n,o) F,v = [ {F,v} ((e,. A Dyo) (.A uJ)] { cp }e (12b) (F,v 'Dyo)cl-,., A e,=(F,v•nyo) [ {bo;} ((e,. 1\ bOi}) [.A.,])] {cp }e (12c) (e,.•ci>, 11 )((F,v, n, 0 , e;.)e,+(ny 0 •e,.)e,. A F,v)

Figure 4 .Figure 5 .Figure 6 .

 456 Figure 4. Circular crack: (a) use of polar co-ordinates (u. v) = (p. 0); (b) usual grid pattern. use of Cartesian co-ordinates

Figure 7 .Figure 8 .

 78 Figure 7. Circular crack: </J 3 normalized by the theoretical maximum value: analytical results; analytical results; -x -x -computed values with the first mesh; -0 -0 -computed values with the second mesh

'Figure 9 .Figure 10 .

 910 Figure 9. Semi-elliptical surface crack in a cylindrical bar. L the surface-line of the crack

  ny 0 )=0=J (kernel containing cl-cr)dudv Acr + vp f (kernel containing «Pcy 1

1 Figure 11 .

 111 Figure 11. Surface crack mesh: 58 elements, 60 interior nodes. The elements are shaped into the elliptic co-ordinates defined on the ellipse

• 3 ,

 3 b/R=0•4, bja=5f7, the stress intensity factors are given by: for a uniform pressure, =0•820 for linear pressure, K.J(amaxfo) = 1/((Jmaxfo) X 2nj2n X lim { E/(16n(l-v 2 ))•c/J3/Jr} r-oo =0•595

Figure 12 .Figure 13 .

 1213 Figure 12. Surface crack under a uniform pressure: c/J 1 =c/J 2 =0, c/J 3 #;0. Crack shapes before and after deformation: (a) front view; (b) spatial view

Figure 14 .

 14 Figure 14. Surface crack under a torsion shear distribution: upper view

  [%] matrical representation of a linear transformation .AI from E 3 into IRn, also Cn denoted more precisely by e [AI'] =(a " b)•c (a, b, c) a®b tensor product defined by Vc E IR 3 , (a® b)c = a(b•c)

Table V .

 V Components of E/(16n(l-v 2 ))•_,cr resulting from a torsion shear distribution with tmax=l:c/>3=0, l/J1, c/>2:-#:0

		Node
		3722	6	0•3371	7	0•3371	8	0•2921
	9	0•3018	10	0•3193	11	0•3287	12	0•3317
	13	0•3287	14	0•3193	15	0•3018	16	0•2921
	17	0•2341	18	0•2394	19	0•2500	20	0•2592
	21	0•2643	22	0•2661	23	0•2643	24	0•2592
	25	Q-2500	26	0•2394	27	0•2341	28	0•1829
	29	Q-1967	30	0•2061	31	0•2129	32	0•2163
	33	0•2175	34	Q-2163	35	Q-2129	36	Q-2061
	37	0•1967	38	Q-1829	39	0•1300	40

•~