Guillaume Dujardin

Asymptotics of linear initial boundary value problems with periodic boundary data on the half-line and finite intervals

This paper deals with the asymptotic behaviour of the solutions of linear initial boundary value problems with constant coefficients on the half-line and on finite intervals. We assume that the boundary data are periodic in time and we investigate whether the solution becomes time-periodic after sufficiently long time. Using the Fokas' transformation method, we show that for the linear Schrödinger equation, the linear heat equation and the linearised KdV equation on the half-line, the solutions indeed become periodic for large time. However, for the same linear Schrödinger equation on a finite interval, we show that the solution, in general, is not asymptotically periodic; actually, the asymptotic behaviour of the solution depends on the commensurability of the time period T of the boundary data with the square of the length of the interval over π.

Introduction

This paper deals with the asymptotic behaviour of the solutions of linear initial boundary value problems (IBVP) with constant coefficients on the half-line and on finite intervals when the boundary data are periodic. Following the recent investigation of Bona and Fokas [START_REF] Bona | Initial-boundary-value-problems for linear and integrable nonlinear dispersive partial differential equations[END_REF], our main concern is the existence of a periodic asymptotic profile for the solutions of such problems. As we shall see, the situation is rather different on the half-line and on finite intervals, because in the case of a finite interval, the "waves" can be "reflected" from one boundary to the other.

We first look at linear problems on the half-line written in the form      ∂ t q + ω(-i∂ x )q = 0 t > 0, x > 0 q(0, x) = q 0 (x)

x > 0 ∂ j x q(t, 0) = f j (t) t > 0, 0 ≤ j ≤ N -1,

(1.1)
where ω is a polynomial of degree n ∈ N ⋆ , q 0 is a smooth initial datum, N ∈ {1, . . . , n} is the number of boundary data such that the problem is well-posed (see [START_REF] Fokas | Initial-boundary value problems for linear dispersive evolution equations on the half-line[END_REF]), and (f j ) 0≤j≤N -1 are N given smooth functions that are compatible with q 0 . We assume that f 0 , . . . , f N -1 are periodic functions of time with the same period and we look for a smooth asymptotic profile q p such that

• for all x ≥ 0, t → q p (t, x) is a periodic function of time,

• for all x ≥ 0, |q(t, x) -q p (t, x)| -→ t→+∞ 0.

In [START_REF] Bona | Initial-boundary-value-problems for linear and integrable nonlinear dispersive partial differential equations[END_REF], the existence of a periodic profile for the solution of the linear KdV equation on the half-line (corresponding to ω(-i∂ x ) = -ik 3 ) is stated, provided that the initial datum q 0 is homogeneous (q 0 ≡ 0) and that the boundary datum f 0 is periodic (N = 1 in that case). The proof is based on the new general transform method developped by A.S. Fokas [START_REF] Fokas | A unified transform method for solving linear and certain nonlinear PDEs[END_REF]. In this paper, we use the Fokas' method to derive the asymptotic behaviour of the solutions of such equations with periodic boundary data. More precisely, for the solutions of the linear Schrödinger equation, of the linear heat equation and of the linear KdV equation on the halfline, we prove the existence of an asymptotic periodic profile when the periodic boundary data are periodic and we provide explicit formulae for these profiles involving the Fourier coefficients of the boundary data.

We also investigate the asymptotic behaviour of the solutions of linear PDEs with constant coefficients on bounded intervals with periodic boundary data of the form:

        
∂ t q + ω(-i∂ x )q = 0 t > 0, x ∈ (0, L) q(0, x) = q 0 (x)

x ∈ (0, L) ∂ j

x q(t, 0) = f j (t) t > 0, 0 ≤ j ≤ N 1 -1 ∂ j

x q(t, L) = g j (t) t > 0, 0

≤ j ≤ N 2 -1, (1.2)
where L > 0 is given, ω is a polynomial of degree n ∈ N ⋆ , q 0 is a smooth initial datum, N 1 , N 2 ∈ {1, . . . , n} are the numbers of data at x = 0 and x = L respectively such that the problem is well-posed (see [START_REF] Fokas | A transform method for linear evolution PDEs on a finite interval[END_REF], Appendix A, and [START_REF] Pelloni | Well-posed boundary value problems for linear evolution equations on a finite interval[END_REF]), and (f j ) 0≤j≤N 1 -1 and (g j ) 0≤j≤N 2 -1 are N 1 + N 2 given smooth functions that are compatible with q 0 . We assume that f 0 , . . . , f N 1 -1 and g 0 , . . . , g N 2 -1 are periodic functions of time with the same period and we look for a smooth asymptotic profile q p such that

• for all x ∈ [0, L], t → q p (t, x) is a periodic function of time,

• for all x ∈ [0, L], |q(t, x) -q p (t, x)| -→ t→+∞ 0.

In the case of the linear Schrödinger equation, we use the formula derived in [START_REF] Fokas | A transform method for linear evolution PDEs on a finite interval[END_REF] to obtain results on the long-time behaviour of the solution. In particular, this formula allows us to give sufficient conditions depending on the link between the length L of the interval and the period T of the boundary data to obtain periodic solutions, as well as solutions that do not have any asymptotic profile.

The outline of the paper is the following. In Section 2, we consider the linear Schrödinger equation on the half-line (i.e. a problem of the form (1.1) with ω(k) = ik 2 ) with homogeneous initial datum (q 0 ≡ 0) and smooth periodic boundary datum f 0 . We recall the Fokas' method (see [START_REF] Fokas | A unified approach to boundary value problems[END_REF]), we prove the existence of a periodic profile q p and we provide an explicit formula for that profile involving the Fourier coefficients of f 0 (see Theorem 2.4, Theorem 2.7 and Theorem 2.14). Section 3 is devoted to another illustration of the efficiency of the integral representation method for the analysis of linear IBVP on the half-line with periodic boundary data through two examples : the linear heat equation (ω(k) = k 2 ) and the linearised KdV equation (ω(k) = -ik 3 ). In both cases, we follow the method we used in Section 2: we use the integral representation method, we prove the existence of an asymptotic profile for the exact solution of the problem and we provide explicit formula for these profiles (see Theorem 3.1, Theorem 3.5 and Theorem 3.6). In Section 4, we consider problems of the form (1.2) on a finite interval (0, L) with ω(k) = ik 2 . If f 0 and g 0 are T -periodic smooth functions, we show that, if L 2 /π and T are linearly dependent on Q, the solution of these linear Schrödinger equations can have, for example, unbounded L 2 -norm ( q(t) L 2 (0,L) -→ t→+∞ +∞) and hence they cannot be asymptotically periodic in general. If T and L 2 /π are linearly independent on Q, we show that the solution is not asymptotically periodic in general as well. Finally, we provide numerical experiments in Section 5 for illustration.

The linear Schrödinger equation on the half-line

In this section, we consider the linear Schrödinger equation on the half-line with a periodic boundary datum and we apply the Fokas' method (see [START_REF] Fokas | A unified approach to boundary value problems[END_REF]). We first look at the asymptotic behaviour of the solution when the boundary datum is a sine function, we prove the existence of a periodic profile in that case and give an explicit formula for that profile using contour deformations of the involved integrals and Cauchy's residue theorem. Then, we apply the same method to the case of a general periodic boundary datum. We prove the existence of a periodic profile in that case and provide an explicit formula for that profile involving the Fourier coefficients of the boundary datum. Recall that we provide numerical experiments in Section 5 for illustration.

The problem

Following [START_REF] Bona | Initial-boundary-value-problems for linear and integrable nonlinear dispersive partial differential equations[END_REF], we consider the linear Schrödinger equation on the half-line written in the form of the inital-boundary-value problem:

     i∂ t q + ∂ xx q = 0 t > 0, x > 0 q(0, x) = 0 x > 0 q(t, 0) = f 0 (t) t > 0. (2.1)
This problem is of the form (1.1) with ω(k) = ik 2 . Note that N = 1 in that case.

Since the problem is linear, we restrict ourselves to the case q 0 ≡ 0. Moreover, we assume that f 0 is a smooth periodic function of time.

The Fokas' method

We apply the method described in [START_REF] Bona | Initial-boundary-value-problems for linear and integrable nonlinear dispersive partial differential equations[END_REF] to Problem (2.1). Our goal is to prove the integral representation formula (2.8) of the solution q in terms of the boundary datum f 0 .

Step 1: Write (2.1) into a divergence form

The first step of the method consists in rewriting the evolution partial differential equation in (2.1) as a family of linear partial differential equations in divergence form. In order to do that, we define for t > 0, x > 0 and k ∈ C

X(t, x, k) := -kq(t, x) + i∂ x q(t, x),
where q denotes the exact solution of Problem (2.1). Then we derive that for all t > 0, x > 0 and k ∈ C,

∂ t e -ikx+ω(k)t q -∂ x e -ikx+ω(k)t X(t, x, k) = 0. (2.2)

Step 2: Obtain a global relation from (2.2)

The next step of the method consists in deriving a global relation involving transforms of the boundary datum f 0 (t) = q(t, 0) and the Neumann derivative of the solution at the boundary ∂ x q(t, 0). In order to do that, we consider the Fourier transform of the solution of Problem (2.1): for all t > 0 and all k ∈ C such that Im (k) ≤ 0, we define q(t, k) = +∞ 0 q(t, x)e -ikx dx, and for all t > 0 and all κ ∈ C, f0 (t, κ) = t 0 e κs q(s, 0)ds and f1 (t, κ) = t 0 e κs ∂ x q(s, 0)ds.

(2.3)

Using the divergence form (2.2) of Problem (2.1), we obtain the following global relation: for all t > 0 and all k ∈ C such that Im (k) ≤ 0,

e ω(k)t q(t, k) = k f0 (t, ω(k)) -i f1 (t, ω(k)).
(2.4)

Step 3: Inverse the global relation (2.4)

Assuming sufficient smoothness and decay on the exact solution q of the Problem (2.1), we derive from (2.4) using the inverse Fourier transform that

∀t, x > 0, q(t, x) = 1 2π +∞ -∞ e ikx-ω(k)t k f0 (t, ω(k)) -i f1 (t, ω(k)) dk. (2.5)

Step 4: Eliminate unknown boundary conditions

Following [START_REF] Bona | Initial-boundary-value-problems for linear and integrable nonlinear dispersive partial differential equations[END_REF], we denote

D = {k ∈ C | Re (ω(k)) < 0},
and we note that for all z ∈ C, we have

z ∈ D ⇐⇒ -z ∈ D.
Therefore, we use the transformation z → -z in (2.4) to get that

∀k ∈ R, ∀t > 0, e ω(k)t q(t, -k) = -k f0 (t, ω(k)) -i f1 (t, ω(k)). (2.6)
This yields

∀k ∈ R, ∀t > 0, f1 (t, ω(k)) = ie ω(k)t q(t, k) + q(t, -k) 2 .
(2.7)

Hence, using Jordan's lemma, we have

∀t, x > 0, 1 2iπ +∞ -∞ e ikx-ω(k)t f1 (t, ω(k))dk = 1 2 q(t, x).
Eventually, plugging this result into the integral representation of the solution formula (2.5) yields

∀t, x > 0, q(t, x) = 1 2π +∞ -∞ 2ke ikx-ik 2 t f0 (t, ω(k))dk.
(2.8)

2.2.5

The integral representation of the solution of (2.1)

In order to sum up our previous calculations, we state the following:

Theorem 2.1 Assume f 0 is a given smooth periodic function of time. Then the solution q of the corresponding problem (2.1) is given by formula (2.8).

The remaining parts of this section are devoted to the use of the integral representation formula (2.8) for the description of the asymptotic behaviour of q. In Section 2.3, we look at the asymptotics when the boundary datum is a sine function. Then we look at the long time behaviour when the function f 0 is any smooth periodic function in Section 2.4.

Sine functions as boundary data for Problem (2.1)

Inspired by [START_REF] Flyer | A hybrid analytical-numerical method for solving evolution partial differential equations. I. The half-line[END_REF], we look at the exact solution q of Problem (2.1) with the simple periodic boundary datum:

f 0 (t) = sin(bt), (2.9)
for some b > 0. We are interested in the long time behaviour of the function q. Our goal is to prove Theorem 2.4 that provides an explicit formula for the asymptotic profile of q. Using the definition of the time-transform of the boundary datum f 0 (see (2.3)), we derive that f0 (t, k) = b + ke kt sin(bt) -be kt cos(bt)

k 2 + b 2 .
We split the corresponding integrand in the integral representation formula (2.8) in the following way. We define for all t, x > 0 and all k ∈ C such that b 2 -k4 = 0,

F 1 (t, x, k) = 2ke ikx be -ik 2 t b 2 -k 4 ,
and

F 2 (t, x, k) = 2ke ikx ik 2 sin(bt) -b cos(bt) b 2 -k 4 .
This way, we have

2ke ikx e -ω(k)t f0 (t, ω(k)) = F 1 (t, x, k) + F 2 (t, x, k). Note that that k → F 1 (t, x, k) + F 2 (t, x, k
) is a holomorphic function on C (singularities in ±i |b| and ± |b| are removable), but k → F 1 (t, x, k) and k → F 1 (t, x, k) have singularities in ±i |b| and ± |b|. Now, we perform two contour deformations on the integral (2.8) in order to derive the asymptotic behaviour of q. Firstly, we define the following complex paths:

∀r < 0, Γ 1 (r) = re -i π 4 , ∀θ ∈ (-π, 0), Γ 2 (θ) = √ b(1 + e iθ ), ∀r > 0, Γ 3 (r) = 2 √ b + r.
This way, the exact solution of Problem (2.1) reads, after contour deformation in formula (2.8)

q(t, x) = 1 2π 0 -∞ (F 1 (t, x, Γ 1 (r)) + F 2 (t, x, Γ 1 (r)))e -i π 4 dr + 1 2π 0 -π (F 1 (t, x, Γ 2 (θ)) + F 2 (t, x, Γ 2 (θ)))i √ be iθ dθ + 1 2π +∞ 0 (F 1 (t, x, Γ 3 (r)) + F 2 (t, x, Γ 3 (r)))dr. (2.10)
Secondly, we use the two following complex paths

∀r > 0, Γ + to derive that 1 2π +∞ 0 F 1 (t, x, Γ 3 (r))dr = 1 2π +∞ 0 F 1 (t, x, Γ - 3 (r))e -i π 4 dr, 1 2π +∞ 0 F 2 (t, x, Γ 3 (r))dr = 1 2π +∞ 0 F 2 (t, x, Γ + 3 (r))e i π 4 dr,
and we can plug these formulae in the last line of (2.10) to get another expression for the exact solution q of Problem (2.1). This way of writing the exact solution, though it is maybe not the simplest one, has two important features. First of all, it decouples the quantities that have different behaviours (some of them tend to 0 when t tends to +∞ (See Lemma 2.2) whereas the others carry the information about the profile of the solution when t tends to +∞ (See Lemma 2.3)). Second of all, the integrals involved in this expression of the exact solution are suitable for numerical computations since they are either integrals of smooth functions over a finite interval, either integrals of smooth exponentially decaying functions over a half-line.

A precise statement on the asymptotic behaviour of the solution q(t, x) of Problem (2.1) with boundary datum (2.9) when x > 0 is given and t tends to +∞ will be derived from the two following lemmas. Lemma 2.2 For all x > 0, the following limits hold:

0 -∞ F 1 (t, x, Γ 1 (r))e -i π 4 dr -→ t→+∞ 0, 0 -π F 1 (t, x, Γ 2 (θ))i √ be iθ dθ -→ t→+∞ 0, +∞ 0 F 1 (t, x, Γ - 3 (r))e -i π 4 dr -→ t→+∞ 0.
Lemma 2.3 The residue theorem ensures that

1 2iπ Γ 1 +Γ 2 +Γ + 3 F 2 (t, x, k)dk = e ix √ b e -ibt -e -x √ b e ibt 2 .
Proof. For all t, x > 0, we have

(k - √ b)F 2 (t, x, k) = 2ke ikx b cos(bt) -ik 2 sin(bt) (k 2 + b)(k + √ b) -→ k→ √ b e ix √ b e -ibt 2 ,
and similarly,

(k -i √ b)F 2 (t, x, k) = 2ke ikx b cos(bt) -ik 2 sin(bt) (k 2 -b)(k + i √ b) -→ k→i √ b - e -x √ b e ibt 2 .
Finally, we derive the following result on the asymtpotic behaviour of the solution q of problem (2.1) with boundary datum (2.9): Theorem 2.4 For all b > 0, let q be the solution of Problem (2.1) with boundary datum (2.9). The following estimate holds for all x > 0:

q(t, x) ∼ t→+∞ ie ix √ b e -ibt -ie -x √ b e ibt 2 .
(2.11)

Proof. Use the results of Lemma 2.2 and Lemma 2.3 in expression (2.10).

General periodic boundary data

In this section, we look at the asymptotic behaviour of the solution q of Problem (2.1) when the boundary datum f 0 is a smooth periodic function of time. Our goal is to prove Theorem 2.7 and Theorem 2.14 that prove the existence of a periodic profile for q and provide explicit formulae for that profile involving the Fourier coefficients of f 0 .

Zero-mean periodic boundary data

Since Equation (2.1) is linear and posed on the half-line, we can assume without loss of generality that the periodic boundary datum f 0 has period 2π. Since

f 0 is smooth, we write ∀t ∈ R, f 0 (t) = n∈Z fn e int , (2.12)
where the Fourier coefficients ( fn ) n∈Z tend to 0 rapidly as |n| tends to +∞. We assume that f 0 is compatible with the homogeneous initial datum at t = 0.

In other words, we assume that

f 0 (0) = n∈Z fn = 0. (2.13)
Moreover, we assume that f 0 has zero mean-value:

2π 0 f 0 (t)dt = 2π f0 = 0. (2.14)
With these hypotheses on the boundary datum g 0 and following the method introduced in Section 2.3, we rewrite the exact solution q of Problem (2.1) (see Lemma 2.5). We first define the three complex paths γ 1 ,γ 2 and γ 3 as follows:

∀r ∈ (-∞, √ 2), γ 1 (r) = re -i π 4 (2.15) ∀r ∈ (1, +∞), γ 2 (r) = -i + r ∀r ∈ ( √ 2, +∞), γ 3 (r) = re -i π 4 .
(2.16) Lemma 2.5 Let q be the exact solution of Problem (2.1) with 2π-periodic smooth boundary datum f 0 written in the form (2.12) and satisfying (2.13) and (2.14).

The following identity holds for all t, x > 0:

q(t, x) = 1 2iπ n∈Z fn γ 1 +γ 2 2ke ikx e int -e -ik 2 t (n + k 2 ) dk.
(2.17)

Remark 2.6 Note that for all t, x > 0, every integrand is a smooth function of

k since the singularities (in k = ± √ -n if n ≤ -1 or in k = ±i √ n if n ≥ 1) are removable.
Theorem 2.7 We have the following asymptotic periodic profile for the exact solution of Problem (2.1) with smooth 2π-periodic boundary data f 0 written in the form (2.12) and satisfying (2.13) and (2.14): for all x > 0,

q(t, x) ∼ t→+∞ ∞ n=1 fn e -x √ n e int + f-n e ix √ n e -int .
(2.18)

Proof. Using identity (2.17) of Lemma (2.5), we write that for all n ∈ Z \ {0},

γ 1 +γ 2 2ke ikx e int -e -ik 2 t (n + k 2 ) dk = γ 1 +γ 2 2ke ikx e int (n + k 2 ) dk - γ 1 +γ 2 2ke ikx e -ik 2 t (n + k 2 ) dk.
Note that, by contour deformation,

γ 1 +γ 2 2ke ikx e -ik 2 t (n + k 2 ) dk = γ 1 +γ 3 2ke ikx e -ik 2 t (n + k 2 ) dk,
and for all x > 0, this quantity tends to 0 when t tends to +∞, uniformly in n. Moreover, using Jordan's lemma, we get

1 2iπ γ 1 +γ 2 2ke ikx e int (n + k 2 ) dk = e ix √ -n e int if n ≤ -1 e -x √ n e int if n ≥ 1.
Remark 2.8 Note that, setting b = 1 in (2.9), Theorem 2.7 implies

∀x > 0, q(t, x) ∼ t→+∞ 1 2i e -x e it - 1 2i e ix e -it .
This result is exactly the corresponding relation (2.11) of Theorem 2.4.

Remark 2.9 Formula (2.18) of Theorem 2.7 ensures that the Sobolev regularity (in time) of the asymptotic profile of the exact solution of Problem (2.1) is the same as that of the periodic datum f 0 (recall that f 0 is assumed sufficiently smooth anyway).

Remark 2.10 Further analysis of the convergence of γ 1 +γ 3 2ke ikx e -ik 2 t (n+k 2 ) dk to 0 shows that the convergence in (2.18) is uniform in x on bounded sets.

General periodic boundary data

In order to get rid of the somehow artificial condition (2.14), we compute the limit as t tends to +∞ of the formal exact solution of Problem (2.1) obtained using formula (2.8) with boundary datum f 0 (t) ≡ 1 in Lemma 2.12. Before that, we define the following complex paths (note that they are a shifted version of the family (γ 1 , γ 2 , γ 3 ) defined in section 2.4.1):

∀r ∈ (-∞, 0), Γ 1 (r) = -i + re -i π 4 ∀r ∈ (0, +∞), Γ 2 (r) = -i + r ∀r ∈ (0, +∞), Γ 3 (r) = -i + re -i π 4 ,
and provide us with another lemma:

Lemma 2.11
The following equality holds:

1 2iπ Γ 1 +Γ 3 e -iz 2 z dz = 1 2 .
Proof. Using contour deformation, the left hand side of the identity reads

1 2iπ -1 -∞ e -ik 2 k dk + 1 2iπ 0 -π ie -ie 2iθ dθ + 1 2iπ +∞ 1 e -ik 2 k dk.
The left and the right term cancel out and the double of the last term is

1 2iπ 0 -π ie -ie 2iθ dθ + 1 2iπ π 0 ie -ie 2iθ dθ = 1 2iπ C e -ik 2 k dk,
where C stands for the positive unit circle. Hence the result is proved using Cauchy's residue theorem:

2 1 2iπ Γ 1 +Γ 3 e -iz 2 z dz = 1 2iπ C e -ik 2 k dk = 1.
We are now able to prove the result about the asymptotic behaviour of the formal exact solution of Problem (2.1) obtained using formula (2.8) with boundary datum f 0 (t) ≡ 1:

Lemma 2.12 For all x > 0, 1 2π +∞ -∞ 2ke ikx e -ik 2 t t 0 e ik 2 s dsdk -→ t→+∞ 1.
Proof. Note that for all t, x > 0, one has

1 2π +∞ -∞ 2ke ikx e -ik 2 t t 0 e ik 2 s dsdk = 1 2iπ +∞ -∞ 2e ikx 1 -e -ik 2 t k dk.
Moreover, since the function k → e ikx 1-e -ik 2 t k is analytic, we have

+∞ -∞ 2e ikx 1 -e -ik 2 t k dk = Γ 1 +Γ 2 2e ikx 1 -e -ik 2 t k dk.
Using the linearity of the integral, we derive that

+∞ -∞ 2e ikx 1 -e -ik 2 t k dk = Γ 1 +Γ 2 2e ikx k dk - Γ 1 +Γ 2 2e ikx-ik 2 t k dk.
The quantity we want to compute is now the sum of two contributions. Using Jordan's lemma, we compute directly the first contribution:

1 2iπ Γ 1 +Γ 2 2e ikx k dk = 2.
Then, we compute the second contribution using another contour deformation:

Γ 1 +Γ 2 2e ikx-ik 2 t k dk = Γ 1 +Γ 3 2e ikx-ik 2 t k dk.
We have, setting z = k √ t,

Γ 1 +Γ 3 e -ik 2 t k e ikx dk = Γ 1 +Γ 3 e -iz 2 z e iz x √ t dz.
Note that e

iz x √ t -→ t→+∞ 1,
and for all t ≥ 1 and r ∈ R,

e -i(re i π 4 -i) 2 re -i π 4 -i e i(re i π 4 -i) x √ t ≤ e -2 √ 2 2 r 1+r √ 2 2 1 2 + r + √ 2 2 2 1/2 e max(0,1+ √ 2 2 r)x .
Since this bound is integrable on R and independent of t ≥ 1, we derive that

Γ 1 +Γ 3 e -iz 2 z e iz x √ t dz -→ t→+∞ Γ 1 +Γ 3 e -iz 2 z dz.
This last integral is equal to iπ by Lemma 2.11. Finally,

1 2iπ Γ 1 +Γ 2 2e ikx-ik 2 t k dk -→ t→+∞ 1,
and the result is proved.

Remark 2.13 This lemma even allows us to consider boundary data f 0 that do not satisfy the compatibility condition (2.13).

Finally, we are able to use the linearity of the problem to combine Theorem 2.7 with Lemma 2.12 to obtain the following theorem: Theorem 2.14 Assume f 0 is a given sufficiently smooth 2π-periodic complexvalued function written in the form (2.12) and satisfying (2.13). We have the following result for the asymptotic behaviour of the exact solution q of Problem (2.1): for all x > 0,

q(t, x) ∼ t→+∞ f0 + ∞ n=1 fn e -x √ n e int + f-n e ix √ n e -int .
Remark 2.15 This result gives an explicit description of the asymptotic profile of the exact solution of the linear Schrödinger equation on the half-line (2.1) for homogeneous initial datum and periodic boundary datum.

The linear heat equation and the linear KdV equation on the half-line

This section is devoted to another illustration of the efficiency of the Fokas' method for the analysis of the asymptotic behaviour of the solutions of linear initial boundary value problems on the half-line. As we did in Section 2, we apply this integral representation method to the solutions of linear evolution partial differential equations with constant coefficients of the form (1.1) on the halfline with periodic boundary data and we use contour deformation and Cauchy's residue theorem to derive the existence of periodic profiles and obtain explicit formulae involving the Fourier coefficients of the boundary data (see Theorem 3.1, Theorem 3.5 and Theorem 3.6). Since the method is very similar to the one used in the previous section (Section 2), the results are presented in a more concise form. We first look at the linear heat equation in Section 3.1. Then we investigate the case of the linear KdV equation in Section 3.2, providing an explicit description of the periodic profile whose existence was stated in [START_REF] Bona | Initial-boundary-value-problems for linear and integrable nonlinear dispersive partial differential equations[END_REF].

The linear heat equation on the half-line

In this section, we investigate the long time behaviour of the solutions of the linear heat equation on the half-line with periodic boundary data. This problem is of the form (1.1) with ω(k) = k 2 . Using the linearity of the problem, we restrict ourselves to the case of homogeneous initial data q 0 ≡ 0. Moreover, we assume that f 0 is a smooth function of t and is 2π-periodic.

The transformation method

We are interested in the long time behaviour of the solution of the following problem:

     ∂ t q -∂ xx q = 0 t > 0, x > 0 q(0, x) = 0 x > 0 q(t, 0) = f 0 (t) t > 0, (3.1)
where f 0 is a smooth 2π-periodic function. Applying the Fokas' transformation method (see Section 2.2), we obtain the following representation formula for the exact solution of Problem (3.1):

∀t, x > 0, q(t, x) = 1 2iπ +∞ -∞ 2ke ikx-k 2 t f0 (t, ω(k))dk, (3.2)
where ω(k) = k 2 and f0 (t, κ) = t 0 e κs f 0 (s)ds.

Periodic boundary data with zero mean

Theorem 3.1 We have the following asymptotic periodic profile for the exact solution of Problem (3.1) with smooth 2π-periodic boundary data f 0 written in the form (2.12) and satisfying (2.13) and (2.14): for all x > 0,

q(t, x) ∼ t→+∞ ∞ n=1 e -x √ n √ 2 2 fn e -ix √ n √ 2 2 e int + f-n e ix √ n √ 2 2 e -int . (3.3)
Moreover, the convergence is uniform with respect to x.

Proof. By a straightforward calculation, we get, for n ∈ Z \ {0},

f0 (t, k 2 ) = e int e k 2 t -1 k 2 + in .
Hence, defining for all t, x > 0,

F 1 (t, x, k) = 2ke ikx -e -k 2 t k 2 + in and F 2 (t, x, k) = 2ke ikx e int k 2 + in , we get, 2ke ikx e -k 2 t f0 (t, ω(k)) = F 1 (t, x, k) + F 2 (t, x, k).
Using (3.2), we derive that for all t, x > 0,

q(t, x) = 1 2iπ +∞ -∞ F 1 (t, x, k)dk + 1 2iπ +∞ -∞ F 2 (t, x, k)dk. One easily checks that 1 2iπ +∞ -∞ F 1 (t, x, k)dk -→ t→+∞ 0,
uniformly in n and x. Moreover, using Jordan's lemma, we have

1 2iπ +∞ -∞ F 2 (t, x, k)dk = e -x √ n √ 2 2 e -ix √ n √ 2 2 e int , if n ≥ 1 and 1 2iπ +∞ -∞ F 2 (t, x, k)dk = e -x √ -n √ 2 2 e ix √ -n √ 2 2 e int ,
if n ≤ -1. This completes the proof.

Remark 3.2 One observes in Theorem 3.1 the following smoothing effect of the heat equation: for any periodic boundary datum f 0 with (sufficient) Sobolev-type regularity, for all x > 0, the asymptotic profile of the exact solution of Problem (3.1) at point x has a Gevrey-type regularity.

General periodic boundary data

As we did for the Schrödinger equation in Section 2.4, we investigate the asymptotic behaviour of the solution of Problem (3.1) with boundary datum g 0 ≡ 1.

Lemma 3.3 Let Γ be the complex path defined for r ∈ R by Γ(r) = r -i. One has

1 2iπ Γ e -z 2 z dz = 1 2 .
Proof. The proof is similar to that of Lemma 2.11.

Lemma 3.4 For all x > 0, 1 2iπ +∞ -∞ 2ke ikx e -k 2 t t 0 e k 2 s dsdk -→ t→+∞ 1.
Proof. The proof is similar to that of Lemma 2.12. The above integral reads by contour deformation

1 2iπ +∞ -∞ 2e ikx 1 -e -k 2 t k dk = 1 2iπ Γ 2e ikx 1 -e -k 2 t k dk,
where Γ is the complex path defined for r ∈ R by Γ(r) = r -i. Note that

1 2iπ Γ 2 e ikx k dk = 2,
by Jordan's lemma and

1 2iπ Γ 2e ikx e -k 2 t k dk -→ t→+∞ 1,
by Lemma 3.3. This completes the proof.

These lemmas allow us to generalise Theorem 3.1 to the case of general periodic boundary conditions: Theorem 3.5 We have the following asymptotic expansion for the exact solution of Problem (3.1) with smooth 2π-periodic boundary data f 0 written in the form (2.12) and satisfying (2.13): For all x > 0,

q(t, x) ∼ t→+∞ f0 + ∞ n=1 e -x √ n √ 2 2 fn e -ix √ n √ 2 2 e int + f-n e ix √ n √ 2 2 e -int . (3.4)
Proof. The proof consists in adding Lemma (3.4) to Theorem (3.1) by linearity.

The linearised KdV equation

In this section, we investigate the long time behaviour of the solutions of the linear Korteweg-de Vries (KdV) equation on the half-line with periodic boundary data. This problem is of the form (1.1) with ω(k) = -ik 3 . Using the linearity of the problem, we restrict ourselves to the case of homogeneous initial data q 0 ≡ 0. Moreover, we assume that f 0 is a smooth function of t and is 2π-periodic.

The transformation method

We are interested in the long time behaviour of the solution of the following problem:

     ∂ t q + ∂ xxx q = 0 t > 0, x > 0 q(0, x) = 0 x > 0 q(t, 0) = f 0 (t) t > 0, (3.5)
where f 0 is a smooth 2π-periodic function. Applying the Fokas' transformation method, we obtain the following representation formula for the exact solution of Problem (3.5):

∀t, x > 0, q(t, x) = 1 2π +∞ -∞ -3k 2 e ikx+ik 3 t f0 (t, ω(k))dk, (3.6)
where ω(k) = -ik 3 and f0 (t, κ) = t 0 e κs f 0 (s)ds.

Periodic boundary data with zero mean

We denote by α the complex number e i 2π 3 . We define the two following complex paths for θ ∈ (0, π 3 ):

∀r ∈ (-∞, 0), Γ 1 (r) = re -i( π 3 -θ) , ∀r ∈ (0, +∞), Γ 2 (r) = re i( π 3 -θ) .
Theorem 3.6 We have the following asymptotic expansion for the exact solution of Problem (3.5) with smooth 2π-periodic boundary data f 0 written in the form (2.12) and satisfying (2.13) and (2.14): for all x > 0,

q(t, x) ∼ t→+∞ ∞ n=1 fn e iα 3 √ nx e int + f-n e -iα 2 3
√ nx e -int .

(3.7)

Remark 3.7 Note that for all λ ∈ R, the complex conjugate of e iαλ is e -iα 2 λ . If the boundary datum f 0 is real valued then for all n ∈ Z, fn = f-n . Hence, the asymptotic profile given by formula (3.7) is real valued as well.

Remark 3.8 Using the definition of α, relation (3.7) also reads

q(t, x) ∼ t→+∞ ∞ n=1 e - √ 3 2 3 √ nx fn e -i 1 2 3 √ nx e int + f-n e i 1 2 3 √ nx e -int . (3.8)
Proof. By a straightforward calculation, we get, for n ∈ Z \ {0},

f0 (t, -ik 3 ) = e int e -ik 3 t -1 i(n -k 3 ) .
Hence, defining for all t, x > 0,

F 1 (t, x, k) = 3k 2 e ikx -e ik 3 t (k 3 -n)
and

F 2 (t, x, k) = 3k 2 e ikx e int (k 3 -n) ,
we get, -3ik 2 e ikx e ik 3 t f0 (t, ω(k)) = F 1 (t, x, k) + F 2 (t, x, k).

Note that for all t, x > 0, F 1 + F 2 is an analytic function of k. Using (3.6) and contour deformation, we derive that for all t, x > 0,

q(t, x) = 1 2iπ Γ 1 +Γ 2 F 1 (t, x, k)dk + 1 2iπ Γ 1 +Γ 2 F 2 (t, x, k)dk. One easily checks that 1 2iπ Γ 1 +Γ 2 F 1 (t, x, k)dk -→ t→+∞ 0,
uniformly in n and x. Moreover, using Jordan's lemma, we have

1 2iπ Γ 1 +Γ 2 F 2 (t, x, k)dk = e iα 3 √ nx e int , if n ≥ 1 and 1 2iπ Γ 1 +Γ 2 F 2 (t, x, k)dk = e iα 2 3
√ nx e int , if n ≤ -1. This completes the proof.

IBVP over the finite interval 4.1 Introduction

In this section, we consider linear evolution PDEs with constant coefficients on a finite interval (0, L) where L is a given positive real number. These problems are of the form (1.2). Such problems have been studied in [START_REF] Fokas | A transform method for linear evolution PDEs on a finite interval[END_REF] by implementing the Fokas' transformation method. We are interested in the long time behaviour of the solutions of these problems when the boundary data are periodic in time. In particular, just as we did for problems on the half line in the previous sections (see Section 2 and Section 3), we investigate the existence of a periodic profile for the solutions of such problems. More precisely, if q(t, x) denotes the solution of such a problem for t ≥ 0 and x ∈ [0, L], we define an asymptotic profile q p (t, x) as a smooth periodic function of t satisfying for all x ∈ [0, L]

|q(t, x) -q p (t, x)| -→ t→+∞ 0.
With the linear Schrödinger equation as an example, we present in this section different types of asymptotic behaviours that occur when the boundary conditions are periodic and share the same period. In contrast to the case of linear problems on the half-line (see Theorem 2.14), it turns out that an asymptotic profile does not exist in general (see Theorem 4.1, Theorem 4.4 and Theorem 4.9).

The linear Schrödinger equation

We investigate the long time behaviour of the solution of the following homogenous linear initial boundary value problem :

         i∂ t q + ∂ xx q = 0 t > 0, x ∈ (0, L) q(0, x) = 0 x ∈ (0, L) q(t, 0) = f 0 (t) t > 0 q(t, L) = g 0 (t) t > 0, (4.1)
where L and T are given positive real numbers and f 0 and g 0 are given smooth T -periodic functions. Note that this problem is of the form (1.2) with ω(k) = ik 2 and N 1 = N 2 = 1. Our analysis is based on the following classical representation formula of the solution of (4.1) (see [START_REF] Fokas | A transform method for linear evolution PDEs on a finite interval[END_REF], equation (3.2)): for all t > 0 and x ∈ (0, L),

q(t, x) = i 4L m∈Z sin(k m x)e -ik 2 m t N (t, k m ) -N (t, -k m ) , k m = mπ L , (4.2)
where the function N is given by

N (t, k) = k( f0 (t, k) -e -ikL g0 (t, k)),
and the functions f0 and g0 are the following t-transforms of f 0 and g 0 : for all k ∈ C and t ≥ 0, f0 (t, k) = t 0 e ik 2 s f 0 (s)ds and g0 (t, k) = t 0 e ik 2 s g 0 (s)ds. (4.3)

Although the exact solution of Problem (4.1) is a smooth function of (t, x), the representation formula (4.2) involves a series whose convergence can be weak (and is usually not uniform on (0, L) for example). However, as an exact representation of the solution in terms of transforms of the boundary data, it provides much information in particular about the long time behaviour. It turns out that the typical behaviour of the exact solution q of (4.1) depends on whether T and L 2 /π are linearly independent on Q or not. As we shall see, this condition corresponds to the possibility for a periodic signal of frequency 2π/T given at a boundary to get reflected on the other boundary in a specific way. Subsection 4.3 is devoted to the case T and L 2 /π are Q-linearly dependent. Subsection 4.4 is devoted to the case T and L 2 /π are Q-linearly independent.

In the following, we set b = 2π/T . Moreover, for all measurable complex functions h defined on (0, L), we define

h L 2 (0,L) = 1 L L 0 |h(x)| 2 dλ(x) 1/2
, and

h L ∞ (0,L) = inf y ∈ R λ x ∈ (0, L) |h(x)| > y = 0 ,
where λ denotes the Lebesgue measure on R. Of course, L 2 (0, L) (respectively L ∞ (0, L)) denotes the space of (classes of) complex functions h on (0, L) such that h L 2 (0,L) < +∞ (resp. h L ∞ (0,L) < +∞).

The dependent case 4.3.1 A simple example

Before dealing with the general dependent case, we investigate a simple example.

In this subsection, we assume that T and L 2 /π are linearly dependent on Q in the following way:

2 T = π L 2 . (4.4)
Hence, we have b = π 2 /L 2 . Moreover, we take as boundary data in Problem (4.1) the following functions:

f 0 (t) = sin(bt) and g 0 (t) = 0, (4.5) for all t ∈ R. Our result on the long time behaviour of q in that case is the following:

Theorem 4.1 Assume the period T of the boundary conditions f 0 and g 0 defined in (4.5) satisfies relation (4.4). Then the exact solution of the linear Schrödinger equation (4.1) on (0, L) satisfies

q(t) L ∞ (0,L) -→ t→+∞ +∞.
Remark 4.2 In particular, the solution of Problem (4.1) with periodic boundary conditions (4.5) does not become asymptotically periodic.

Proof. Note that formula (4.2) reads

q(t, x) = i 2L m∈Z k m sin(k m x)e -ik 2 m t f0 (t, k m ), k m = mπ L . (4.6)
As in Section 2.3, using (4.5), we derive that

f0 (t, k) = b + ik 2 e ik 2 t sin(bt) -be ik 2 t cos(bt) b 2 -k 4 and g0 (t, k) = 0. Hence, f0 (t, ± π L ) = L 2 1 -e 2i π 2 L 2 t 4π 2 + i t 2 , (4.7)
and for all m ∈ Z such that |m| ≥ 2,

f0 (t, m π L ) = L 2 π 2 1 + e im 2 π 2 L 2 t im 2 sin( π 2 L 2 t) -cos( π 2 L 2 t) 1 -m 4 . (4.8)
Using Lemma 4.3, we get that q(t)

2 L 2 (0,L) is the sum of 1 2 π 2 L 2 | f0 (t, k -1 )| 2 + | f0 (t, k 1 )| 2 + f0 (t, k -1 ) f0 (t, k 1 ) + f0 (t, k 1 ) f0 (t, k -1 ) , (4.9) and 1 2 |m|≥2 |k m f0 (t, k m )| 2 - |m|≥2 k m f0 (t, k m )k m f0 (t, -k m ) . (4.10)
Using Hölder's inequality, we have that the real number (4.10) is non-negative. Moreover, since for all k ∈ C, f0 (t, k) = f0 (t, -k), we have that the real number (4.9) is equal to

2 π 2 L 2 | f0 (t, k 1 )| 2 .
Finally, (4.7

) implies that | f0 (t, k 1 )| -→ t→+∞ +∞. Hence, q(t) L 2 (0,L) -→ t→+∞ +∞,
and so q(t) L ∞ (0,L) -→ t→+∞ +∞.

Lemma 4.3 For all sequences of complex numbers (a k ) k∈Z such that k∈Z |a k | 2 < +∞, the series

k∈Z a k sin kπ x L ,
converges in L 2 (0, L) and its sum f satisfies

1 L L 0 |f (x)| 2 dx = 1 2 k∈Z |a k | 2 - 1 2 k∈Z a k a -k .
Proof. Essentially, the proof is based on the identity valid for all k 1 , k 2 ∈ Z:

1 L L 0 sin k 1 π x L sin k 2 π x L dx = 1 2 δ k 1 ,k 2 -δ k 1 ,-k 2 ,
where δ is the Kronecker's symbol (i.e. δ a,b = 1 if a = b and δ a,b = 0 otherwise).

The general dependent case

In this subsection, we still assume that T and L 2 /π are linearly dependent on Q. Since Problem (4.1) is linear and homogeneous, one can easily deduce the asymptotic behaviour of the general case (where both f 0 and g 0 are non-zero) from the one where f 0 ≡ 0 or g 0 ≡ 0. Of course, while doing so, one has to keep in mind that the possible explosions as t → +∞ (see Theorem 4.1) arising from each boundary datum can interact in a destructive way and indeed lead to a bounded solution. Consider for example the case

f 0 (t) = g 0 (t) = sin(bt) with b = π 2 /L 2 = 2π/T . In that case, N (t, k ±1 ) = 0 and for |m| ≥ 2, N (t, k m ) is a 2L 2 /π-periodic function of t (see (4.8)). Therefore, q is itself a 2L 2 /π-periodic function of t.
Hence, we can focus without too much loss of generality on the asymptotic behaviour of the solution of (4.1) with g 0 ≡ 0. Any linear relation on Q between T and L 2 /π yields a relation of the form If there exists (n, m) ∈ R such that fn = 0, then the solution q of Problem (4.1)

α 2π T + β π 2 L 2 = 0, ( 4
satisfies q(t) L ∞ (0,L) -→ t→+∞ +∞,
and hence is not asymptotically periodic. Otherwise, q is a periodic function of period 2L 2 /π × max(1, -α/β).

Proof. Formula (4.6) can be rewritten

q(t, x) = i L m≥1 k m sin(k m x)e -ik 2 m t f0 (t, k m ), k m = mπ L . (4.12)
Note that for all (n, m)

∈ Z 2 , (n, m) ∈ R ⇐⇒ k 2 m + nb = 0 . For all m ≥ 1, all n ∈ Z and all t ≥ 0, t 0 e ik 2 m s e ibns ds =      e ik 2 m t e inbt -1 i(k 2 m + nb) if (n, m) ∈ R t if (n, m) ∈ R.
Note that, since f 0 is smooth, the coefficients fn tend to zero rapidly as |n| → +∞. Using (4.12), we write

q(t, x) = q p (t, x) + q r (t, x),
where

q p (t, x) = i L m≥1 k m sin(k m x)e -ik 2 m t n∈Z|(n,m) ∈R fn e ik 2 m t e inbt -1 i(k 2 m + nb)
,

and q r (t, x) = t i L m≥1 k m sin(k m x)e -ik 2 m t n∈Z|(n,m)∈R fn .
Note that q p is a periodic function and 2L 2 /π × max(1, -α/β) is one if its periods. In particular, q p (t) L 2 (0,L) is bounded. Moreover, for all m ≥ 1, n∈Z|(n,m)∈R fn contains at most 1 term. Using Lemma 4.3, we derive that 1

L L 0 |q r (t, x)| 2 dx = t 2 2L 2 k≥1 k 2 m n∈Z|(n,m)∈R fn 2 .
This proves the result.

The independent case

In this section, we consider the long time behaviour of the solution q of Problem (4.1) on (0, L) with g 0 ≡ 0 and f 0 a T -periodic smooth function such that T and L 2 /π are linearly independent on Q. Under these assumptions, we shall prove that q is not asymptotically periodic (see Section 4.1 for a definition) in general (see Theorem 4.9). In order to prove this result, we recall the following:

Theorem 4.5 (Kronecker's approximation theorem) For all n ∈ N ⋆ , for all (α 1 , . . . , α n ) ∈ R n , for all ε, T > 0 and all linearly independent real numbers θ 1 , . . . , θ n , there exists a real number t > T and n integers h 1 , . . . , h n such that ∀i ∈ {1, . . . , n},

|tθ i -h i -α i | < ε.
See [START_REF] Apostol | Modular functions and Dirichlet series in number theory[END_REF], Theorem 7.9 and Exercise 7 of Chapter 7 for a proof of this result.

Theorem 4.6 Assume that f , g and h are periodic complex functions defined on R with h = f + g.

Assume that f and g have positive smallest periods. We denote these periods by T 1 and T 2 respectively. Assume that T 1 and T 2 are linearly independent on Q. Assume that f is bounded. Then f does not have any open interval of continuity.

See [START_REF] Olmsted | On the sum of two periodic functions[END_REF], Theorem 5 for a proof of this result.

Preliminary results

Lemma 4.7 Assume that f , g and h are periodic complex-valued functions. Let us denote respectively T 1 , T 2 and T 3 any of their positive periods. Assume that f and g are continuous functions on R and that the function

d(t) = f (t) + g(t) -h(t),
tends to 0 when t tends to +∞. Then h is continuous on R.

Proof. Assume ε > 0 and t 0 ∈ R are given. Since d(t) -→ t→+∞ 0 and T 3 > 0, there exists n ∈ N such that

∀s ≥ t 0 -1, |d(s + nT 3 )| < ε/4. (4.13)
We infer that for all s ≥ t 0 -1,

|h(t 0 ) -h(s)| = |h(t 0 + nT 3 ) -h(s + nT 3 )| ≤ |f (t 0 + nT 3 ) -f (s + nT 3 )| + |g(t 0 + nT 3 ) -g(s + nT 3 )| + ε/2.
Since f and g are continuous at t 0 + nT 3 , we derive that there exists η > 0 such that for all s ∈ (t 0 -η, t 0 + η),

|f (t 0 + nT 3 ) -f (s + nT 3 )| < ε/4 and |g(t 0 + nT 3 ) -g(s + nT 3 )| < ε/4.
Hence, for all s ∈ R such that |s -t 0 | < min(1, η), we have

|h(t 0 ) -h(s)| ≤ ε.
Therefore, h is continuous in t 0 .

The next theorem states that, if a sum of two continuous periodic functions with uncommensurable smallest positive periods is asymptotically periodic, then it is periodic. Theorem 4.8 Assume that f , g and h are periodic complex-valued functions defined on R. Assume that f and g are continuous on R and that they admit smallest positive periods denoted by T 1 and T 2 respectively. Assume that T 1 and T 2 are linearly independent on Q and that the function

d(t) = f (t) + g(t) -h(t), (4.14)
tends to 0 when t tends to +∞. Then d ≡ 0. In other words,

∀t ∈ R, f (t) + g(t) = h(t). (4.15)
Proof. Let T 3 denote any of the positive periods of h. Using Lemma 4.7, the function h is continuous on R. Since f , g and h are periodic continuous functions on R, they are uniformly continuous on R. Since T 1 and T 2 are positive Q-linearly independent real numbers, so are 1/T 1 and 1/T 2 . This proof is divided in two parts, depending on whether

1/T 3 ∈ Span Q (1/T 1 , 1/T 2 ) or not. Assume that 1/T 3 ∈ Span Q (1/T 1 , 1/T 2 ).
There exists a, b, c ∈ Z such that

a 1 T 3 = b 1 T 2 + c 1 T 1 ,
and one can assume without loss of generality that a > 0. Hence, (4.18) Since 1/T 1 and 1/T 2 are linearly independent, so are T 1 and T 2 . Hence, the Kronecker's approximation theorem (Theorem 4.5) ensures that there exists t > T and x 1 , x 2 ∈ Z such that

aT 1 T 2 = bT 1 T 3 + cT 2 T 3 . (4.16)
|tT 1 -x 1 | < η max(aT 2 , |b|T 3 , |c|T 3 ) and |tT 2 -x 2 | < η max(aT 1 , |b|T 3 , |c|T 3 ) . (4.19) Relation (4.14) ensures that d(s) = d(s + atT 1 T 2 ) + δ 1 t (s) + δ 2 t (s) -δ 3 t (s)
, where In order to get a bound on δ 3 t (s), we note that, using (4.16), we have

δ 1 t (s) = f (s) -f (s + atT 1 T 2 ), δ 2 t (s) = g(s) -g(s + atT 1 T 2 ), ( 4
atT 1 T 2 = btT 1 T 3 + ctT 2 T 3 .
Therefore, using (4.21), Since s ∈ R and ε > 0 are arbitrary, this proves that (4.15) holds true when

δ 3 t (s) = h(s) -h(s + bT 3 tT 1 + cT 3 tT 2 ). Hence, defining δ 4 t (s) = h(s) -h(s + bT 3 x 1 + cT 3 tT 2 ),
1/T 3 ∈ Span Q (1/T 1 , 1/T 2 ). Assume now that 1/T 3 ∈ Span Q (1/T 1 , 1/T 2 ).
In that case, 1/T 1 , 1/T 2 and 1/T 3 are linearly independent on Q. After multiplication by T 1 T 2 T 3 , this implies that T 1 T 2 , T 2 T 3 and T 1 T 3 are linearly independent on Q. Assume ε > 0 is given. We choose η > 0 as we did in (4.18) from the uniform continuity of f , g and h on R. Fix s ∈ R. As before, since d(t) -→ t→+∞ 0, we can choose T > 0 big enough to ensure that |d(s + tT 1 T 2 T 3 )| < ε. Using Kronecker's approximation theorem (Theorem 4.5), we get the existence of t > T and (x 1 , x 2 , x 3 ) ∈ Z 3 such that This proves (4.15) when 1/T 3 ∈ Span Q (1/T 1 , 1/T 2 ).

|tT 1 T 2 -x 3 | < η/T 3 , |tT 1 T 3 -x 2 | < η/T 2 ,
T 2 T 3 -x 1 T 1 | < η, so that | δ1 t (s)| = |f (s) -f (s + tT 1 T 2 T 3 )| = |f (s + x 1 T 1 ) -f (s + tT 1 T 2 T 3 )| < ε, using the T 1 -periodicity of f . Similarly, we have | δ2 t (s)| = |g(s) -g(s + tT 1 T 2 T 3 )| = |g(s + x 2 T 2 ) -g(s + tT 1 T 2 T 3 )| < ε, and
| δ3 t (

Long time asymptotics in the independent case

Theorem 4.9 Let us denote by q the smooth solution of problem (4.1) corresponding to a smooth periodic boundary datum f 0 with smallest period T > 0 and g 0 ≡ 0. Assume that T and L 2 /π are linearly independent on Q. For all x ∈ (0, L), the function t → q(t, x) is not asymptotically periodic.

See Section 4.1 for a definition of asymptotic periodicity. Proof. For all t > 0, we have, using (4.2)

q(t, x) = i 2L k∈Z k m sin(k m x)e -ik 2 m t f0 (t, k m ), (4.28)
in L 2 (0, L), where f0 is defined in (4.3). Since f 0 is T -periodic, we write the Fourier expansion of f 0 , setting as before b = 2π/T :

∀t ∈ R, f 0 (t) = n∈Z fn e inbt .
Since T and L 2 /π are Q-linearly independent, we have for all m, n ∈ Z

k 2 m + nb = π 2 L 2 m 2 + n 2π T = π 2 T L 2 m 2 T + 2n L 2 π = 0.
We derive that for all m ∈ Z,

f0 (t, k m ) = n∈Z fn e ik 2 m t e inbt -1 i(k 2 m + bn) .
Hence, using (4.28), we can define the functions

f x (t) = 1 2L k∈Z k m sin(k m x) n∈Z fn e inbt k 2 m + nb , and
g x (t) = - 1 2L k∈Z k m sin(k m x) n∈Z fn e -ik 2 m t k 2 m + nb ,
to have the following decomposition of the function q:

∀x ∈ (0, L), ∀t ∈ R, q(t, x) = f x (t) + g x (t).
Note that f x is T -periodic and g x is 2L 2 /π-periodic so that q is an almost periodic function (see [START_REF] Besicovitch | Almost periodic functions[END_REF] and [START_REF] Bohr | Almost periodic functions[END_REF] for a definition). Note that, since f 0 is smooth, f x is a smooth function of time. Since q is also a smooth function of time, so is g x . Assume that for all x ∈ (0, L) there exists a continuous periodic profile q p (t, x) such that d x (t) = f x (t) + g x (t) -q p (t, x) -→ t→+∞ 0.

Fix x 0 ∈ (0, L). Since T and 2L 2 /π are linearly independent on Q, so are the smallest positive periods of the continuous functions f x 0 and g x 0 . Hence, Theorem 4.8 ensures that d ≡ 0. This yields ∀t ∈ R, f x 0 (t) + g x 0 (t) = q p (t, x 0 ).

Finally, Theorem 4.6 ensures that f x 0 does not have any open interval of continuity. This is a contradiction. We derive that t → q(t, x 0 ) is not asymptotically periodic.

5 Numerical experiments 5.1 Linear IBVP on the half-line

The Schrödinger equation

We first compute the numerical solution q of the homogeneous linear initial boundary value problem (2.1) on the half-line, with the boundary datum f 0 = sin. This fits the framework of Section 2.3 with b = 1 in formula (2.9). For any given positive values of t and x, the numerical value of q(t, x) is computed by integrating the integrand in relation (2.8) on the line γ 1 + γ 3 (see (2.15),(2.16)). Since this function and all its derivatives have exponentially decaying moduli, we truncate the integral to a bounded interval with an error of at most ε/2 where ε is a given positive (small) given real number. Then, we use the trapezoidal method to compute an approximation of the integral on the bounded interval with sufficiently many steps to ensure an error of at most ε/2.

We plot the results obtained for t ∈ [0, 20] and s ∈ [0, 7] on Figure 1 (real part of q) and Figure 2 (imaginary part of q). In order to illustrate the convergence of the solution q to the periodic profile provided by Theorem 2.4 and Theorem 2.14, we plot the real part of q and of the profile on Figure 3 at point x = 2.0 as functions of time. Of course, the numerical values of q are computed as before.

The heat equation

We now compute the numerical solution q of the homogeneous linear initial boundary value problem (3.1) on the half-line, with the boundary datum f 0 = sin. This fits the framework of Section 3.1.

To compute numerical values for q, we use contour deformation just as we did for the Schröginger equation above. We plot the solution we obtain on Figure 4 for t ∈ [0, 20] and x ∈ [0, 7].

In order to illustrate the convergence of the solution q to the periodic profile provided by Theorem 3.5, we plot q and the profile on Figure 5 at point x = 2.0 as functions of time.

The linearised KdV equation

As a last illustration for the asymptotics of linear IBVP problems on the half-line with periodic boundary data, we compute the solution q of the linear problem (3.5) with the boundary datum ∀t ≥ 0, f 0 (t) = sin(t) + 1 10 sin(10t).

(5.1)

On Figure 6, we plot the numerical results obtained for t ∈ [0, 30] and x ∈ [0, 4].

In order to illustrate the convergence of the solution q to the periodic profile provided by Theorem 3.6, we plot q and the profile on Figure 7 at point x = 2.0 as functions of time.

The linear Schrödinger equation on bounded intervals

The dependent case

We investigate the asymptotic behaviour of the solutions of the IBVP (4.1) on the finite interval (0, L) with L = π. We compute the numerical values of the solution q using a finite difference method scheme in both time and space, with a CFL number 0.05 and 50 discretisation points in space.

Firstly, we consider periodic boundary data f 0 and g 0 defined in (4.5) with b = 1 (i.e. T = 2π). Note that relation (4.4) holds true. We plot on Figure 8 the numerical values of q L 2 (0,L) as a function of time for t ∈ [0, 50]. We observe that the L 2 -norm of q is equivalent to a linear function of time in that case and hence tends to +∞. This is consistent with Theorem 4.1. Note that, with the notations of Theorem 4.4, we have α = 1 and β = -1 in that case. Moreover, f 0 is 2π-periodic and its only non-zero Fourier coefficients are f±1 = ±1/(2i). Hence, (1, -1) ∈ R and f-1 = 0. Therefore, the numerical results of Figure 8 also illustrate part of Theorem (4.4).

Secondly, we consider periodic boundary data f 0 and g 0 defined in (4.5) with b = 2 (i.e. T = π). Note that relation (4.11) holds, with α = 1 and β = -2. The only non-zero Fourier coefficients of f 0 are f±1 = ±1/(2i) (see their definition in Theorem 4.4). Hence, for all (n, m) ∈ R, fn = 0. Theorem 4.4 ensures that q is a periodic function of time in that case, with period 2π. This is illustrated by Figure 9 where the numerical values of q L 2 (0,L) are plotted as a function of time for t ∈ [0, 50] in that case.

The independent case

We investigate the asymptotic behaviour of the solutions of the IBVP (4.1) on the finite interval (0, L) with L = 1. We compute the numerical values of the solution q using a finite difference method scheme in both time and space, with a CFL number 0.05 and 50 discretisation points in space.

We consider the boundary data f 0 and g 0 given by (4.5) with b = 1 (i.e T = 2π). Theorem 4.9 ensures that the solution q is not asymptotically periodic in that case. This is illustrated on Figure 10 where the L 2 -norm of q is plotted as a function of time for t ∈ [0, 50]. The global behaviour of the solution is illustrated on Figure 11 where the real part of the solution is plotted for t ∈ [0, 50] and x ∈ [0, 1]. Moreover, we plot on Figure 12 the real part of q(t, x) evaluated at x = 0.5 as a function of time. Figure 10, Figure 11 and Figure 12 show that, even if the solution has some stability property (the L 2 -norm seems to remain bounded in that case), no periodic asymptotic behaviour seems to take place. This illustrates the result of Theorem 4.9.

  .11) with (α, β) ∈ (Z \ {0}) 2 relatively primes, α ≥ 1, β ≤ -1. Such a couple (α, β) is uniquely determined. Note that the hypotheses of Theorem 4.1 correspond to (α, β) = (1, -1). Theorem 4.4 Let T > 0 denote a period of the smooth function f 0 , and set g 0 ≡ 0. Define b = 2π/T . Denote ( fn ) n∈Z the Fourier coefficients of f 0 : ∀t ∈ R, f 0 (t) = n∈Z fn e inbt Assume that T and L 2 /π are linearly dependent on Q so that (4.11) holds. Let us denote R = (n, m) ∈ Z × Z αm 2 -βn = 0 .

Fix

  s ∈ R. Assume ε > 0 is given. Since d(t) -→ t→+∞ 0 and aT 1 T 2 > 0, we can choose T > 0 such that ∀t > T, |d(s + taT 1 T 2 )| < ε. (4.17) Moreover, since f ,g and h are uniformly continuous on R, there exists η > 0 such that for all t, u ∈ R with |t -u| < η, one has |f (t) -f (u)| < ε and |g(t) -g(u)| < ε and |h(t) -h(u)| < ε.

  s)| = |h(s) -h(s + tT 1 T 2 T 3 )| = |h(s + x 2 T 2 ) -h(s + tT 1 T 2 T 3 )| < ε.Using the triangle inequality in (4.27), we obtain |d(s)| < 4ε.

Figure 1 :

 1 Figure 1: Real part of the solution of Problem (2.1) with boundary datum (2.9) and b = 1.

Figure 2 :

 2 Figure 2: Imaginary part of the solution of Problem (2.1) with boundary datum (2.9) and b = 1.

Figure 3 :

 3 Figure 3: Solution of Problem (2.1) at point x = 2.0 with boundary datum (2.9) and b = 1 (full line) and the asymptotic profile given by Theorem 2.14 (dotted line).

Figure 4 :

 4 Figure 4: Solution of Problem (3.1) with boundary datum (2.9) and b = 1.

Figure 5 :

 5 Figure 5: Solution of Problem (3.1) at point x = 2.0 with boundary datum (2.9) and b = 1 (full line) and the asymptotic profile given by Theorem 3.5 (dotted line).

Figure 7 :

 7 Figure 7: Solution of Problem (3.5) with boundary datum f 0 (t) = sin(t) + sin(10t)/10 (full line) and the corresponding periodic profile provided by Theorem 3.6 (dotted line).

Figure 8 :

 8 Figure 8: L 2 -norm of the solution of Problem (4.1) on (0, π) with boundary data f 0 (t) = sin(t) and g 0 ≡ 0 as a function of time.

Figure 9 :

 9 Figure 9: L 2 -norm of the solution of Problem (4.1) on (0, π) with boundary data f 0 (t) = sin(2t) and g 0 ≡ 0 as a function of time.

Figure 10 :

 10 Figure 10: L 2 -norm of the solution of Problem (4.1) on (0, 1) with boundary data f 0 (t) = sin(t) and g 0 ≡ 0 as a function of time.

Figure 11 :

 11 Figure 11: Real part of the solution of Problem (4.1) on (0, 1) with boundary data f 0 (t) = sin(t) and g 0 ≡ 0.

Figure 12 :

 12 Figure 12: Real part of the solution of Problem (4.1) on (0, 1) with boundary data f 0 (t) = sin(t) and g 0 ≡ 0 evaluated at x = 0.5 as a function of time.

  (s + ax 2 T 1 ) -f (s + atT 1 T 2 ). + ax 2 T 1 ) -(s + atT 1 T 2 )| = aT 1 |x 2 -tT 2 | < η. + ax 1 T 2 ) -g(s + atT 1 T 2 ),and, using (4.19),|(s + ax 1 T 2 ) -(s + atT 1 T 2 )| = aT 2 |x 1 -tT 1 | < η,

	Since f is T 1 -periodic, we have
	δ 1 t (s) = f (4.22)
	Moreover, using (4.19), we have
	|(s Hence, using (4.18) and (4.20), we have
	|δ 1 t (s)| < ε.	(4.23)
	Similarly, we have
	δ 2 t (s) = g(s and hence
	|δ 2 t (s)| < ε.	(4.24)
		.20)
	and
	δ 3 t (s) = h(s) -h(s + atT 1 T 2 ).	(4.21)

  and δ 5 t (s) = h(s + bT 3 x 1 + cT 3 tT 2 ) -h(s + bT 3 tT 1 + cT 3 tT 2 ), Note that, since h is T 3 -periodic and b, c, x 1 , x 2 ∈ Z,δ 4 t (s) = h(s + (bx 1 + cx 2 )T 3 ) -h(s + bT 3 x 1 + cT 3 tT 2 ).

	we have,
	δ 3 t (s) = δ 4 t (s) + δ 5 t (s).
	Using (4.18) and (4.19), we derive that
	|δ 4 t (s)| < ε	and	|δ 5 t (s)| < ε.
	This yields
		|δ 3 t (s)| < 2ε.	(4.25)
	We derive from the inequalities (4.17), (4.23), (4.24) and (4.25) that
		|d(s)| < 5ε.

  and|tT 2 T 3 -x 1 | < η/T 1 . -h(s + tT 1 T 2 T 3 ).From (4.26), we derive that |tT 1

			(4.26)
	Using (4.14), we derive that
		d(s) = d(s + tT 1 T 2 T 3 ) + δ1 t (s) + δ2 t (s) -δ3 t (s),	(4.27)
	with
		δ1 t (s) = f (s) -f (s + tT 1 T 2 T 3 ),	δ2 t (s) = g(s) -g(s + tT 1 T 2 T 3 ),
	and	δ3 t = h(s)

(r) = 2 √ b + re i π

, ∀r > 0, Γ - 3 (r) = 2 √ b + re -i π 4 ,

Acknowledgments

The author thanks A.S. Fokas and P.A. Markowich for their ideas and comments on this work. This publication is based on work supported by Award No. KUK-


B Nedjar

D T Nguyen

On a simple cyclic plasticity modeling with implicit kinematic hardening restoration

Keywords: Kinematic hardening restoration, Viscoplasticity, Continuum thermodynamics, Large number cylic loading, Numerical implementation

This paper presents a fully three-dimensional plastic constitutive modeling framework suitable for the prediction of cyclic loading at large number of cycles. It can require only one yield surface and it is motivated by a simple rheological model where a restoration of the kinematic hardening is introduced. The classical kinematic hardening rules are then simply adapted leading to time-dependent evolution laws that are consistent with continuum thermodynamics requirements. The resulting behavior is physically motivated by many man-made materials of engineering interest such as bituminous material. This framework allows all types of yield functions to be easily implemented numerically. This is first illustrated with algorithmic details through a simple associative pressure-insensitive model example of the von Mises type. Then a more elaborated model is given where the present framework is applied to the description of bituminous materials submitted to triaxial static creep and to large number of cyclic loadings. Of particular interest is the ratcheting and the mean stress relaxation. The responses agree well with some experimental test results found in the literature.

Introduction

Plastic accumulation is one of the major distresses of materials and structures under cyclic loadings. This phenomenon designates in general the apparition and accumulation of permanent deformations with the increasing number of loading cycles. For many materials, this accumulation may not cease until fracture. In the particular domain of pavement engineering, the ratcheting of asphalt concrete is the origin of the rutting phenomenon, an important distress mode under repeated traffic loading.

During the last years, many efforts have been devoted to model the ratcheting responses of various materials. These works concentrated especially on metal-like materials, see for example [START_REF] Lemaitre | Mechanics of Solid Materials[END_REF][START_REF] Ohno | Current state of the art in constitutive modeling for ratcheting behavior[END_REF][START_REF] Ohno | Kinematic hardening rules with critical state of dynamic recovery, part I: formulations and basic features for ratcheting behavior[END_REF][START_REF] Yoshida | A constitutive model of cyclic plasticity[END_REF] among many others. In general, these models postulate single surface elastic domain with complex nonlinear kinematic and isotropic hardening rules and have proven their relevance on the modeling of such metal-like materials. On the other hand, the multimechanisms and multi-criteria type models, generalized by Mandel in [START_REF] Mandel | Généralisation de la théorie de w.t. koiter[END_REF], are also able to reproduce complex phenomena under cyclic loadings. Here again, the relevance of these approaches has been clearly demonstrated, for instance in [START_REF] Cailletaud | Study of plastic/viscoplastic models with various inelastic mechanisms[END_REF][START_REF] Khoei | On the implementation of multi-surface kinematic hardening plasticity and its applications[END_REF][START_REF] Nguyen | Cyclic elasto-viscoplastic model for asphalt concrete materials[END_REF] to mention a few. However, the application of multi-surface models to complex problems is not obvious, mostly because of the numerical difficulties due to singular points when multiple pressure-sensitive yield criteria are used. This was experienced for instance by the authors in [START_REF] Nguyen | Prédiction des déformations irréversibles des couches de surface des chaussées bitumineuses[END_REF][START_REF] Nguyen | Cyclic elasto-viscoplastic model for asphalt concrete materials[END_REF]. Note that, in general, the use of pressure-sensitive yield criteria is necessary when modeling geomaterials, concretes or bituminous materials. In particular for these latters, advances have been reported by a number of authors, see for example [START_REF] Scarpas | Finite element simulation of damage development in asphalt concrete pavements[END_REF][START_REF] Seibi | Constitutive relations for asphalt concrete under high rates of loading[END_REF][START_REF] Sides | A comprehensive visco-elastoplastic characterization of sand-asphalt under compression and tension cyclic loading[END_REF][START_REF] Tashman | A microstructure-based viscoplastic model for asphalt concrete[END_REF] among others. The relevant long literature has contributed significantly to the understanding of the mechanisms influencing the asphalt concrete performance as well, and most of the proposed models are more adapted to the predictive material responses under monotonic loadings. A model that is able to predict at the same time the material responses under monotonic and large number cyclic loadings is still a challenging task.

The constitutive modeling framework proposed in this paper is precisely motivated by physical observations on bituminous materials. In such mixtures, the permanent deformations occur as a result of slight movements between the aggregates at the meso-scale. These movements are accompanied by the viscous flow of the asphalt binder. As a consequence, the material response is viscoplastic at the macroscopic level. Moreover, it is also observed that the viscous flow of the binder relaxes the auto-stresses due to the incompatibility of the plastic deformations. This relaxation phenomenon constitutes the point of departure of the present formulation resulting in the so-called concept of kinematic hardening restoration. That is, the kinematic hardening stresses (the back-stresses in the continuum thermodynamics terminology) still relax even when the plastic straining ceases. The relaxation process developed in this work is based on a differential evolution equation not resorting to history variables, as that was the case in earlier models found in the literature. In fact, the concept of kinematic hardening restoration is not new. It has already been used for metals, precisely, by using history variable models, see for instance [START_REF] Lemaitre | Mechanics of Solid Materials[END_REF]. The framework we develop in this paper may be applied not only to bituminous materials, but also to other materials such as polymers and even for metals as well.

An outline of the remainder of the paper is as follows. In section 2, the formulation is firstly motivated by a simple rheological model where the governing equations are derived by elementary considerations. Of particular interest, it is shown that this model is able to exhibit plastic accumulation under both static and cyclic loadings. A straightforward three-dimensional generalization is given in Section 3 following the basic concepts of the continuum thermodynamics. Then, in section 4, a simple J 2 flow-like model example of the von Mises type is detailed from the theoretical and numerical points of view to show how plastic accumulations at large number of cycles can be simply implemented by standard numerical schemes. In section 5, a more elaborated model is briefly presented. This latter is adapted to bituminous materials and a set of numerical results are compared with some experimental responses of [START_REF] Brown | A fundamental study of the stress-strain characteristics of a bituminous materials[END_REF]. Finally, conclusions and perspectives are drawn in section 6. Noteworthy remarks and comments are given throughout this paper.

Motivation: a one-dimensional rheological model

We first illustrate the basic idea behind the present modeling framework through a one-dimensional study. For this, let us examine the different rheological models depicted in Figure 1.

• The classical Prager model of Figure 1(a) is the simplest one that is capable to produce time-independent permanent strains with kinematic hardening. It consists of a spring of constant E > 0 in series with a plastic device composed by a friction element with constant σ y > 0 in parallel with a spring of constant H > 0 (the kinematic hardening modulus in our case). • In Figure 1(b), a parallel dashpot of constant η > 0 is added to the last plastic device. The resulting model is then no more than an extension of the Prager model that renders the plastic evolution time-dependent. • Finally, in the model depicted in Figure 1(c), we simply replace the last spring of constant H by a Maxwell element. That is, the stress in the spring of constant H is relaxed by a dashpot of constant η X > 0. Notice that this replacement can also be done on the Prager model of Figure 1(a), i.e. without the dashpot of constant η. In all cases, the global constitutive behavior will always be time-dependent by construction. Now the governing equations for the rheological model of Figure 1(c) are derived by completely elementary considerations. Let σ denotes the total stress acting on it, ε the total strain, and ε p the plastic strain. Moreover, let α be the strain in the spring of constant H, and β the strain in the dashpot of constant η X . One has then: where the (back-)stress X represents the kinematic hardening. The elastic domain f , represented by the friction element of constant σ y , is here written as f = |σ -X| -σ y , and the plastic flow is given by:

σ = E(ε -ε p ), X = Hα, ε p = α + β (1) E E E H H H σ y σ y σ y η η η X ε ε p α β (a) (b) (c)
εp = γ ∂f ∂σ (2)
where the scalar γ can be given by a viscoplastic regularization equation of the Perzyna type [START_REF] Perzyna | Fundamental problems in viscoplasticity[END_REF]:

γ = f η (3)
where . defines the positive part function as f = 1 2 (f + |f |). Now, assuming a linear relation connecting the stress X with its corresponding strain rates β, i.e. X = η X β, the evolution of the strain-like variable α is obtained from the rate form of Equation. (1) 3 :

α = εp - X η X (4)
Analyzing the above evolution equations, it follows that for situations where loadings are such that f ≤ 0, Equation (3) gives γ = 0 leading to εp = 0 by Equation (2), and Equation (4) becomes then:

α = - X η X (5)
That is, the kinematic hardening continues to evolve even when the stress state belongs to the interior of the elastic domain, i.e. with f ≤ 0. In this case, the back-stress X is restored through the following evolution equation obtained by combining Equation [START_REF] Germain | Continuum thermodynamics[END_REF] with rate form of the constitutive relation [START_REF] Brown | A fundamental study of the stress-strain characteristics of a bituminous materials[END_REF]

2 : Ẋ + H η X X = 0 (6)
Remark 1. If the dashpot of constant η is removed from the rheological model of Figure 1(c), the scalar γ in Equation ( 2) becomes a plastic multiplier determined by the classical Kuhn-Tucker loading/unloading conditions: γ ≥ 0, f ≤ 0 and γf = 0, together with the consistency condition γ ḟ = 0. In this case the evolution equations ( 4)-( 6) remain the same.

⊓ ⊔

As an illustration, Figure 2 shows the results of a static and a cyclic creep semi-analytical computations obtained with the constitutive equations ( 1) -(4). The used material parameters are: E = 3000 M P a, σ y = 0.1 M P a, H = 80 M P a, η = 2500 M P a s and η X = 10 6 M P a s.

For the static creep of Figure 2(a), a stress of value σ = 0.25 M P a has been quickly applied at time t = 0 and then maintained. One can observe that the plastic strain ε p continue growing with time. And in Figure 2(b), a sinusoidal loading has been applied with an amplitude of σ max = 0.25 M P a at a frequency of 1 Hz. The evolution of the plastic strain ε p is plotted versus the number of cycles. One can observe that the plastic accumulation grows contiuously even at very large number of cycles.

Hence, both the static and cyclic creep are simultaneously described by the same rheological model where the kinematic hardening restoration is in turn controlled by the viscous parameter η X through the evolution equation ( 4).

Three-dimensional generalization and continuum thermodynamics

The extension of the precedent one-dimensional developments to the threedimensional constitutive modeling is straightforward. As usual in the infinitesimal theory, the total strain tensor ε is additively split into an elastic part ε e and a plastic part ε p : ε = ε e + ε p . From the continuum thermodynamics point of view, the state of the material is characterized by its free energy function denoted here by ψ. In our context, this latter is chosen of the form ψ ≡ ψ(ε e = εε p , α) where the tensor quantity α is the strain-like internal variable related to kinematic hardening. Without loss of generality, isotropic hardening will not be considered for the sake of clarity. This latter can in fact be included in the formulation in a usual manner.

In the simplest case, the elastic domain E is defined by a single-surface yield criterion of the form

E = {(σ, X) / f (σ, X) ≤ 0} (7)
T ime (s) where f is a convex function, σ is the stress tensor, and X denotes the stresslike internal variable conjugated to α in the sense that

X = ∂ψ ∂α (8)
this latter being in fact the set of back-stress tensors. By taking the time derivative of ψ, inserting the result into the Clausius-Duhem inequality, i.e. D = σ : εψ ≥ 0, and making use of standard arguments in continuum thermodynamics [START_REF] Germain | Continuum thermodynamics[END_REF], one obtains the classical state law

σ = ∂ψ ∂ε e (9)
together with the reduced dissipation inequality

D = σ : εp -X : α ≥ 0 ( 10
)
where the symbol : denotes the double contraction, i.e. A :

B = tr [AB T ] ≡ A ij B ij (
summation on repeated indices is assumed), and tr[.] designates the trace operator. Now motivated by the rheological model of Figure 1(c), the plastic strain rate is in turn additively split in the form εp = α + β where β is the viscous kinematic hardening strain rate tensor. The inequality [START_REF] Nguyen | Cyclic elasto-viscoplastic model for asphalt concrete materials[END_REF] becomes then:

D = (σ -X) : εp + X : β ≥ 0 ( 11
)
In view of the expression of the first term in [START_REF] Ohno | Current state of the art in constitutive modeling for ratcheting behavior[END_REF], the yield surface defining the elastic domain in ( 7) is restricted to the general form f ≡ f (σ -X). For the kinematic hardening restoration, we assume the existence of a potential denoted by Ω r (X), i.e. a convex and positive function such that Ω r = 0 ⇔ X = 0. Then, the evolution equations can be given by

       εp = γ ∂f ∂(σ -X) ≡ γ ∂f ∂σ β = ∂Ω r ∂X =⇒ α = εp - ∂Ω r ∂X (12)
With this flow rule, the thermodynamic requirement ( 11) is satisfied for any admissible process thanks to the convexity of the functions f and Ω r together with the positiveness of the multiplier γ. The evolution equation ( 12) 1 is the three-dimensional extension of its one-dimensional counterpart given in [START_REF] Cailletaud | Study of plastic/viscoplastic models with various inelastic mechanisms[END_REF].

Likewise, with a quadratic potential Ω r = X : X/2η X , (12) 2 extends (4) to three-dimensions.

A J 2 -flow model example

To make matters as concrete as possible, we consider in detail a simple plastic model of the von Mises type with a pressure insensitive yield criterion. We denote by s and e the deviatoric stress and strain tensors, respectively

s = σ -p1, e = ε - 1 3 tr[ε]1 (13)
where 1 denotes the second order unit tensor and p = 1 3 tr[σ] the hydrostatic stress. Considering again only kinematic plastic hardening, the yield criterion is chosen as

f (σ, X) = ||s -X|| - 2 3 σ y ( 14
)
with σ y being the flow stress and where the norm of a second order tensor is defined as ||.|| = [(.) : (.)] 1/2 . By construction, the back-stress tensor is purely deviatoric, i.e. tr[X] = 0. The constitutive relations ( 8) -(9) are given by

σ = C : (ε -ε p ), X = Hα ( 15
)
where C is the rank four Hooke's elasticity tensor (λ and µ being the Lamé coefficients), and the scalar H is the plastic kinematic hardening modulus considered constant in this example. In other words, the free energy function ψ is given by the simple quadratic form

ψ = 1 2 ε e : C : ε e + 1 2 Hα : α ( 16
)
The model is completed by appending the evolution equations of the form given by ( 12) for the internal variables:

εp = γ s -X ||s -X|| = n α = γ n - X η X ( 17
)
where n is the unit normal to the yield surface, and where again the potential Ω r (X) = X : X/2η X has been used. In view of the numerical resolution, instead to deal with the evolution equation ( 17) 2 , we prefer its equivalent conjugate form by use of the rate form of the constitutive relation (15

) 2 : Ẋ = γH n - H η X X (18)

Numerical approximation

The algorithmic approximation of the local evolution equations (17) 1 and ( 18) is based on the nowadays well-known elastic predictor/plastic corrector scheme. From the finite element point of view, this procedure is accomplished locally at each integration point. Let [t n , t n+1 ] be a typical time interval. For an arbitrary material point, the known (internal) variables {ε p n , X n } are considered as prescribed initial data. Then, within the elastic predictor step, the trial state is obtained by freezing the plastic flow ignoring the constraint placed by the plastic yield criterion. The trial state is obtained by mere function evaluation as

s trial n+1 = s n + 2µ∆e n+1 , p trial n+1 = p n + κ∆ε n+1 : 1 (19)
where κ = λ + 2 3 µ is the bulk modulus, ∆ε n+1 is the total strain increment and ∆e n+1 its deviatoric part. The yield criterion is then evaluated as

f trial n+1 = ||s trial n+1 -X n || - 2 3 σ y (20)
• On the one hand, if f trial n+1 > 0, then γ > 0 and an incrementally plastic straining process is taking place. A plastic correction has to be performed to get the final solution as

   s n+1 = s trial n+1 -2µ∆γ nn+1 X n+1 = X n + H∆γ nn+1 - H η X ∆tX n+1 (21)
where ∆γ = γ∆t with ∆t = t n+1 -t n , and where an implicit backward Euler finite difference scheme has been employed for the evolution equations (17) 1 and (18). Equation ( 21) 2 can be conveniently rewritten as

X n+1 χ = X n + H ∆γ nn+1 with χ = η X η X + H∆t (22)
Here the tensor nn+1 at the final state has the remarkable property that it is also given by the known trial state as

nn+1 = s n+1 -X n+1 ||s n+1 -X n+1 || ≡ s trial n+1 -χX n ||s trial n+1 -χX n || (23)
Finally, the unknown multiplier ∆γ is obtained by enforcing the consistency condition written as f n+1 = (η/∆t)∆γ to avoid numerical difficulties when η → 0, see for exemple [START_REF] Simo | Computational Inelasticity[END_REF]. It is explicitely given by ∆γ =

||s trial n+1 -χX n || - 2 3 σ y 2µ + χH + η ∆t (24)
where the relation ||s n+1 -X n+1 || = ||s trial n+1 -χX n || -∆γ(2µ + χH) has been used.

• On the other hand, if f trial n+1 ≤ 0, then γ = 0 and the trial state is the final solution with also ε p n+1 = ε p n , except for the back-stress tensor where the kinematic hardening restoration takes place through the following evolution equation, see Equation ( 18):

Ẋ = - H η X X, with X| t=tn = X n (25)
This latter can be approximated numerically by an exponential scheme as

X n+1 = X n exp - H η X ∆t . ( 26
)
For the sake of clarity, Table 1 summarizes the conceptual steps involved during this local resolution procedure.

Representative numerical examples

As an illustration, we consider a (5 × 5 × 10) mm 3 parallelepipedic sample discretized with (4 × 4 × 4) classical 8-nodes hexahedral elements. The specimen is firstly submitted to a cyclic uniform pressure on the top face with an amplitude of 200 KP a at a period T = 40 s. The rest of the boundary conditions are such that the bottom face is fixed vertically (along the e 3 direction) while leaving the displacements free on its plane, and the lateral faces of the sample are completely free. Figure 3 shows the resulting stress-strain curves Table 1 Local algorithm for the plastic evolution.

1.

Trial state:

e n+1 = ε n+1 -1 3 (ε n+1 : 1)1 e p n = ε p n -1 3 (ε p n : 1)1 ≡ ε p n s trial n+1 = 2µ(e n+1 -e p n ) p trial n+1 = κ(ε n+1 -ε p n ) : 1 ≡ p n+1 f trial n+1 = ||s trial n+1 -Xn|| -2 3 σy 2.
Check the yield condition:

IF f trial n+1 > 0 THEN let χ = η X η X + H ∆t compute:                  nn+1 = s trial n+1 -χXn ||s trial n+1 -χXn|| ∆γ = ||s trial n+1 -χXn|| -2 3 σy 2µ + χH + η ∆t update:        ε p n+1 = ε p n + ∆γ nn+1 s n+1 = s trial n+1 -2µ∆γ nn+1 X n+1 = χXn + χH ∆γ nn+1 go to step 3 ELSE IF f trial n+1 ≤ 0 THEN update:          ε p n+1 = ε p n s n+1 = s trial n+1 X n+1 = Xn exp - H η X ∆t go to step 3 END IF

3.

Compute the total stress tensor:

σ n+1 = s n+1 + p n+1 1
(in the e 3 direction) of two different computations: the first one with η = 0 and the second one with η = 750 M P a s. The common material parameters used for both the computations are: E = 7500 M P a, ν = 0.3, σ y = 0.06 M P a, H = 250 M P a and η X = 45 10 3 M P a s. Notice that even for η = 0, the material response is still time-dependent through the kinematic hardening restoration process (see also the comments of Remark 1). Irrespective to the value of the parameter η, plastic accumulation grows with the number of cycles.

To show again the ability of the present model for more complex loading conditions, a non-homogeneous cyclic loading is this time simulated on the same sample. Here a pressure of amplitude 400 KP a at a period T = 40 s is applied on one fourth of the top face as shown in conditions on the other faces and the same material parameters as for the precedent computations are used with η = 750 M P a s. Figure 4 depicts the evolution of resultant force versus the maximum displacement of the upper right nodes.

Remark 2. Notice that from all the preceding results, at constant loading amplitudes and periods, plastic accumulation evolves linearly with the number of cycles and can therefore grow indefinitely. However, the materials of engineering interest show that this is not the case at very large numbers of cycles. For instance, one can add isotropic hardening to reach accomodation when a certain plastic equivalent strain is reached. Another way is to sophisticate the model by using non linear hardening rules. This latter is partly what is done in the next application to model the behavior of bituminous materials.

⊓ ⊔

A model adapted to bituminous materials

In this section, a more elaborated model is proposed within the preceding framework. Having in mind the behavior of bituminous materials, experimental observations have shown that these latters are compressible and the permanent deformations are usually caused by volumetric densification, specifically for insufficiently compacted materials, followed by shearing at constant volume, see for example [START_REF] Brown | A fundamental study of the stress-strain characteristics of a bituminous materials[END_REF][START_REF] Eisenmann | Influence of wheel load and inflation pressure on the rutting effect at asphalt-pavements -experiments and theoretical investigations[END_REF]. Hence, a two-mechanisms model is proposed via two distinct kinematic stress-like internal variables: a purely deviatoric tensor X 1 for shearing and a scalar X 2 for the volumetric plasticity. We postulate the following yield criterion which takes the form of a Camclay-like model adapted here for our purposes:

f = ||s -X 1 || 2 + 2 3 κ 2 (p -X 2 + δR 0 ) 2 - 2 3 R 0 ( 27
)
This criterion has the form of an elliptic closed surface of which the initial properties are entirely determined by the parameters R 0 , κ and δ. R 0 is the flow stress, κ characterizes the contribution of the confining stresses, and δ is introduced in order to ensure the dissimmetry of the behavior between tension and compression. For an illustration, Figure 5 shows this elastic domain at the initial state where, in the (p, q) plane drawing, the quantity q is the so-called von Mises equivalent stress. Associated plasticity is not suitable for geomaterials as it could not correctly capture the influence of the confining stresses on the volumetric strain evolution, see for example [START_REF] Florea | Nonassociated elastic/viscoplastic model for bituminous concrete[END_REF][START_REF] Zienkiewicz | Associated and non-associated viscoplasticity in soils mechanics[END_REF]. We postulate then a non-associated plasticity by introducing a flow potential F given by

F = ||s -X 1 || 2 + 2 3 ̟ 2 (p -X 2 + δR 0 ) 2 + a ||s -X 1 || + b |p -X 2 | ( 28
)
where ̟ is a newly introduced parameter for the confining stresse part, and a and b are two parameters that control the nonlinearity of the kinematic hardenings. The associated plasticity could be reached by setting ̟ = κ and a = b = 0. Finally, the kinematic hardening restoration is chosen to depend solely on the (deviatoric) back-stress X 1 via the following potential:

Ω r = (1 -H(f )) 1 2η X X 1 : X 1 ( 29
)
where H(.) is the Heaviside function whose argument is the value of the yield function ( 27). That is, Ω r = 0 when f ≥ 0, and when f < 0, Ω r is positive and the restoration mechanism happens only for pure shear.

The local evolution of the internal variables are then given by the following flow rules of the Perzyna type:

               εp = f η ∂F ∂σ Ẋ1 = H 1 ėp -H 1 ∂Ω r ∂X 1 Ẋ2 = H 2 3 tr[ εp ] (30)
where use has been made of the decomposition ε p = e p + tr[ε p ]1. H 1 and H 2 are the kinematic hardening moduli for the deviatoric and volumetric parts, respectively.

• For f > 0, we have Ω r = 0 from Equation (29), and the precedent evolution equations has a classical form. • For f ≤ 0, the local evolution equations (30) reduce to

εp = 0, Ẋ1 + H 1 η X X 1 = 0, Ẋ2 = 0 (31)
For the sake of clarity, details about the numerical discretization are here omitted this model, see for this [START_REF] Nguyen | Prédiction des déformations irréversibles des couches de surface des chaussées bitumineuses[END_REF]. Its finite element implementation follows the main steps adopted for the model example of Section 4.

A numerical parametric study

In order to demonstrate the performance of the present model and to justify the choices for the material parameters, a set of simple numerical examples are presented below. For this, we consider a cylindrical sample under triaxial loading with lateral confining pressure:

σ =   σ a σ c σ c   (32)
where σ a and σ c are respectively the axial and the confining stresses. The volumetric plastic strain ε p v and the shear plastic strain ε p s are then computed as follows:

ε p =   ε p a ε p c ε p c   , ε p v = tr[ε p ], ε p s = 2 3 (ε p a -ε p c ) (33)
The material parameters are chosen as follows: the Young's modulus E = 3000 M P a, the Poisson ratio ν = 0.35, and for the yield criterion (30), R 0 = 0.1 M P a, κ = 0.9 and δ = 0.75.

Both static and cyclic creep tests are considered. The volumetric and shear plastic strains as functions of the loading durations (time or number of loading cycles) are studied. For the static creep, the applied stresses are σ a = 0.267 M P a and σ c = 0.167 M P a. And for the cyclic creep, the precedent values are the amplitudes at a frequency of 1 Hz (that is, the mean stresses are σ a = 0.133 M P a and σ c = 0.083 M P a).

Figure 6 shows the responses for two particular cases of kinematic hardening: case 1 with H 1 = 0 and H 2 = 80 M P a, and case 2 with H 1 = 65 M P a and H 2 = 0. In both cases, the parameters appearing in the plastic potential F in Equation (28) are ̟ = 1.15, a = b = 0, the viscous parameters are η = 0 and η X = 2.65 10 2 M P a s. We observe that when activating only one kinematic hardening, the model is unable to maintain simultaneously the evolution of volumetric and shear strains under static as well as under cyclic creep loadings. The introduction of both of the two kinematic hardening mechanisms is then necessary for the modeling.

On the other hand, fixing the values of the kinematic hardening moduli H 1 and H 2 (here with H 1 = 65 M P a and H 2 = 80 M P a), the nonlinearities of both of the kinematic hardenings are required. This is proved by the responses of Figure 7 for particular cases of the parameters a and b, i.e. the two parameters controlling the nonlinear evolutions of X 1 and X 2 , respectively (again see Equation ( 28)).

One can observe as well from Figure 7 that these nonlinearities induce a linear evolution of the permanent strain versus the loading duration. But this evolution trend is not conform to the experimental observations in general. It is then necessary to introduce one or several parameters which evolve with the loading history. As the permanent strain rates are inversely proportional to the viscous parameter η, η is thus chosen to depend on the cumulative permanent strain ēp which is widely used to characterize the plastic strain history, see for example [START_REF] Simo | Computational Inelasticity[END_REF][START_REF] Simo | Numerical analysis and simulation of plasticity[END_REF]:

ε p v (%) ε p v (%) ε p s (%) ε p s (%)
ε p v (%) ε p v (%) ε p s (%) ε p s (%) a = b = 0 a = 1, b = 0 a = 0, b = 20 a = 1, b = 20
η = η 0 (1 + η 1 ēp exp[η 2 ēp ]) , ēp (t) = t τ =0 3 2 εp : εp dτ (34)
where η 0 is the initial viscous parameter, η 1 and η 2 are positive parameters.

Along with loading duration, ēp increases, η increases too leading then to the decrease of the permanent strain rates. Likewise, to control the evolution of cyclic permanent strains, we also postulate a progressive change for the restoration parameter η X as follows:

η X = η X0 (1 + η X1 (exp[η X2 ēp ] -1)) ( 35
)
where η X0 , η X1 and η X2 are positive parameters. We use these new enhancements in the next modeling example.

Parameter identification

Of concern in this section, some indications of how the various material parameters can be identified. By first noticing that the three parameters R 0 , κ and δ characterize the elastic domain, see Figure 5, they can then be determined in a usual manner from elastic limits for different loading paths.

Now by considering the model response from the precedent parametric study and its further extensions given by the constitutive relations (34) and (35), one finds that the ε p v (t) and ε p s (t) curves for the static creep and the ε p v (N ) and ε p s (N ) curves for the cyclic creep can be divided into two parts. The first part, that we denote by Phase 1, corresponds to the early response where the permanent strains increase quickly. The next part, that we denote by Phase 2, consists to a somehow stabilized phase where the permanent strains increase slowly.

Parameters H 1 , H 2 , a, b, η 1 and η 2 mostly influence the whole response in static creep. In particular, H 1 , a, H 2 and b strongly influence the Phase 1 responses, while Phase 2 is mostly influenced by η 1 and η 2 . On the other hand, for cyclic creep, the volumetric response is mostly influenced by the parameters ̟ and η 0 in Phase 1, and by the parameter η 2 in Phase 2.

At last, as the constitutive relation (35) solely influences the cyclic creep responses, mostly for the shear response, the paramaters η X0 , η X1 and η X2 can be determinined from the ε p s (N ) curves. With these observations at hand, we propose the following identification procedure by curve fitting of different permament evolution responses:

-Parameters ̟ and η 0 are determined from Phase 1 of static and cyclic creep tests.

- -Parameters η X0 , η X1 and η X2 are determined by curve fitting of the whole cyclic creep response ε p s (N ). These identifications need then triaxial experimental tests under static and cyclic creep, preferably at loading levels of the same order.

A comparison with selected experimental data

In this section, some simulations of creep tests are compared to some experimental results by [START_REF] Brown | A fundamental study of the stress-strain characteristics of a bituminous materials[END_REF]. These authors have in fact realized a complete test campaign to study permanent behaviors of asphalt mixtures. All the tests were processed on cylindrical specimens of length 240 mm. Both static and cyclic creep tests were performed under triaxial loadings with the presence of lateral confining stresses in order to limit the instability of the samples. In cyclic creep tests, the square waves stress form was involved with the stress on for 0.5 s and off for 0.5 s. The stress levels in the corresponding static creep tests were equal to the peak values in the cyclic creep tests. Axial and radial permanent strains were measured and used to calculate the volumetric and shear strains after Equations (33).

The static and cyclic creep tests at stress levels of σ a = 0.267 M P a and σ c = 0.167 M P a were chosen to validate the proposed model. The Young's modulus and the Poisson ratio of the materials were [START_REF] Brown | A fundamental study of the stress-strain characteristics of a bituminous materials[END_REF]: E = 830 M P a and ν = 0.25. Within the present framework, the following set of parameters has been identified, see [START_REF] Nguyen | Prédiction des déformations irréversibles des couches de surface des chaussées bitumineuses[END_REF] for details: κ = 0.9; δ = 0.75; R 0 = 0.10 M P a; ̟ = 1.15 H 1 = 65 M P a; H 2 = 80 M P a; a = 1.80; b = 19.95 η 0 = 265 M P a s; η 1 = 15; η 2 = 1700 η X 0 = 25 M P a s; η X 1 = 0.05; η X 2 = 1750 With these material parameters, Figure 8 shows the numerical responses of the model under static creep loading. Thanks to the presence of two different kinematic hardening mechanisms, both the volumetric and shear evolutions are in satisfactory good agreement with the experimental results. The non-linearity of the kinematic hardening stresses allow the strains rates to be maintained until a high loading duration. The evolution of the viscous parameter η brings a diminution of permanent strains rates, and thus a nonlinear evolution of ε p v and ε p s versus time. Also, a comparison of the numerical responses and the experimental results in the case of cyclic creep tests is presented in Figure 9. Here again, a good correlation confirms the capacities of the proposed model. The restoration of shear kinematic hardening stress X 1 ensures a higher permanent strain level in cyclic creep tests than in static creep tests.

Conclusion and perspectives

In this paper we have presented a viscoplastic constitutive framework suitable for the prediction of the cyclic loadings at high number of cycles. The model is based on the notion of kinematic hardening restoration that is physically observed in many man-made materials of engineering interest.

Motivated by a simple rheological model with elementary considerations, we have shown that the three-dimensional constitutive modeling is straightforward and in agreement with the continuum thermodynamics requirements. We have shown that using a single-surface yield criterion, one can get plastic accumulation at a very large number of cycles.

An application to describe the behavior of bituminous materials has been proposed. We have shown that this latter is richly featured and is capable of simulating not only cyclic creep loadings, but also static creep loadings at the same time. However, as bituminous materials are very complex materials in general, the proposed model should be improved to account for many other characteristics such as the omnipresent temperature dependency in the asphalt concrete pavements, and damage in terms of crack and air void growth.

Fig. 1

 1 Fig. 1 One-dimensional rheological models: (a) the Prager model, (b) the Bingham model, (c) the model for the present formulation.

Fig. 2

 2 Fig. 2 Responses of the rheological model of Figure 1(c): (a) static creep, (b) cyclic creep.

Figure 4 .Fig. 3

 43 Fig. 3 Homogeneous cyclic loading of amplitude 200 KP a and period T = 40 s with: (a) η = 0, and (b) η = 750 M P a s.

Fig. 4

 4 Fig.4Non-homogeneous cyclic loading of amplitude 400 KP a and period T = 40 s.

Fig. 5

 5 Fig. 5 Initial elastic domain.

Fig. 6

 6 Fig.6Responses of the model for two particular cases of H 1 and H 2 (case 1: with H 1 = 0 and H 2 = 80 M P a, and case 2: with H 1 = 65 M P a and H 2 = 0).

Fig. 7

 7 Fig. 7 Responses of the model for particular cases of the parameters a and b.

Fig. 8

 8 Fig. 8 Static creep test. Comparison with experimental results of Brown and Cooper [1].

Fig. 9

 9 Fig. 9 Cyclic creep test. Comparison with the experimantal results of Brown and Cooper [1].

  Paramaters H 1 , a, H 2 and b are determined by curve fitting ε p v (t) and ε p s (t) under static creep. -Parameter η 1 is determined from Phase 1 of the volumetric response ε p v (N ) under cyclic creep. -Parameter η 2 is determined from Phase 2 responses of the ε p

	ε p v (N ) curves.	v (t), ε p s (t) and