D Gasperini

H P Beise

U Schroeder

X Antoine

C Geuzaine

A multi-harmonic finite element method for scattering problems with small-amplitude boundary deformations

Keywords: high frequency scattering, moving boundary, Doppler effect, multi-harmonic resolution, finite element method

A finite element method in the frequency domain is proposed for solving scattering problems with moving or, more generally, deforming boundaries. First, the original problem is rewritten as an equivalent weak formulation set in a fixed domain. Next, this formulation is approximated as a simpler weak form based on asymptotic expansions when the amplitude of the movements or the deformations is small. Fourier series expansions of some geometrical quantities and of the solution are next introduced to obtain a coupled multi-harmonic frequency domain formulation. Standard finite element methods can then be applied to solve the resulting problem and a block diagonal preconditioner is proposed to accelerate the Krylov subspace solution of the linear system for high frequency problems. The efficiency of the resulting method is demonstrated on a radar sensing application for the automotive industry.

Introduction

In the framework of wave scattering problems, the Doppler effect [START_REF] Chen | The Micro-Doppler Effect in Radar, 2nd Ed[END_REF][START_REF] Gill | The Doppler Effect: An Introduction to the Theory of the Effect[END_REF] characterizes the property that the motion of a target modulates the frequency of the reflected signal initially emitted by a source. For motions with uniform velocity, the Doppler frequency shift can be easily obtained [START_REF] Censor | Scattering of electromagnetic waves in uniformly moving media[END_REF] while for more general movements it is estimated by combining simple models and signal processing techniques [START_REF] Cammenga | Micro-Doppler target scattering[END_REF][START_REF] Cammenga | High range resolution micro-Doppler analysis[END_REF][START_REF] Chen | The Micro-Doppler Effect in Radar, 2nd Ed[END_REF][START_REF] Chen | Micro-Doppler effect in radar: phenomenon, model, and simulation study[END_REF][START_REF] Gu | A hybrid radar-camera sensing system with phase compensation for random body movement cancellation in Doppler vital sign detection[END_REF][START_REF] Li | Radar remote monitoring of vital signs[END_REF][START_REF] Liu | Analysis of phase noise influence on micro-Doppler feature extraction on vibrating target[END_REF]. During the last years, the radar detection of non uniformly moving scatterers was applied to a wide variety of problems because of the availability of newly designed high-frequency sensors and devices. For example, the two radar frequency ranges 24-24.5 GHz and 61-61.5 GHz (ISM bands) are now standard in many applications, and the new 77-81 GHz band is being used for applications in the automotive industry. Furthermore, the THz frequency range will be used in the next coming years (e.g. at the frequency 140 GHz). A crucial advantage of high-frequency radar sensing is its sensitivity to micro movements (see e.g. [START_REF] Chen | The Micro-Doppler Effect in Radar, 2nd Ed[END_REF][START_REF] Chen | Micro-Doppler effect in radar: phenomenon, model, and simulation study[END_REF]) of obstacles involving several moving parts, which results in the micro-Doppler effect [START_REF] Chen | The Micro-Doppler Effect in Radar, 2nd Ed[END_REF][START_REF] Chen | Micro-Doppler effect in radar: phenomenon, model, and simulation study[END_REF] used nowadays in many applications [START_REF] Agnihotri | Effect of frequency on micro-Doppler signatures of a helicopter[END_REF][START_REF] Brooks | Temporal deep learning for drone micro-Doppler classification[END_REF][START_REF] Cammenga | High range resolution micro-Doppler analysis[END_REF][START_REF] Chen | Micro-Doppler effect in radar: phenomenon, model, and simulation study[END_REF][START_REF] Garcia-Rubia | Analysis of moving human micro-Doppler signature in forest environments[END_REF]. In the automotive industry, micro-Doppler sensing has recently been applied with success [START_REF] Dias Da | A theoretical investigation of the detection of vital signs in presence of car vibrations and radar-based passenger classification[END_REF][START_REF] Gu | A hybrid radar-camera sensing system with phase compensation for random body movement cancellation in Doppler vital sign detection[END_REF][START_REF] Karahasanovic | Radar-based detection of thoracoabdominal asynchrony during breathing using autocorrelation function analysis[END_REF][START_REF] Li | Radar remote monitoring of vital signs[END_REF] to the contactless detection of vital signs, in particular for the breathing of children left alone on the back seat of overheating cars. Difficult challenges are then related to this kind of applications as for example the analysis of random body movements or vehicle vibrations [START_REF] Dias Da | A theoretical investigation of the detection of vital signs in presence of car vibrations and radar-based passenger classification[END_REF][START_REF] Gu | A hybrid radar-camera sensing system with phase compensation for random body movement cancellation in Doppler vital sign detection[END_REF][START_REF] Li | Radar remote monitoring of vital signs[END_REF][START_REF] Liu | Analysis of phase noise influence on micro-Doppler feature extraction of vibrating target[END_REF][START_REF] Liu | Analysis of phase noise influence on micro-Doppler feature extraction on vibrating target[END_REF][START_REF] Peng | A sinusoidal frequency modulation Fourier transform for radar-based vehicle vibration estimation[END_REF] that can be classified thanks to their radar signature by deep learning techniques [START_REF] Brooks | Complex-valued neural networks for fully-temporal micro-Doppler classification[END_REF][START_REF] Brooks | Temporal deep learning for drone micro-Doppler classification[END_REF][START_REF] Chen | High-resolution sparse representation of micro-Doppler signal in sparse fractional domain[END_REF]. In the development life cycle of these new sensors, a full realistic simulation of the high frequency scenarios is therefore needed. The goal of the present paper is to design a suitable numerical modelling to this problem and to propose efficient computational solutions.

In practice, the physical model leads to solving a time-dependent wave propagation problem in a complex environment (the interior of a car here) with a complex geometry made up of several dielectric materials strongly interacting with the high frequency emitted wave. Moreover, the moving targets, such as a child located on the back seat of the car, can have nontrivial shapes, consist of different materials, often show small amplitude displacements, and usually happen at frequencies several orders of magnitude below the frequencies of the emitter. Therefore, the natural mathematical modelling involves a system of partial differential equations (PDEs) which needs to be efficiently solved numerically for complex geometries and materials. Finally, more complicated situations like random vibrations can be included by adapting the PDE system.

The solution of moving target problems formulated in a PDE framework has already received some attention from both the mathematical and engineering sides. Among them, analytical methods for simple motions were designed in [START_REF] Censor | Scattering of electromagnetic waves by a cylinder moving along its axis[END_REF][START_REF] Censor | Non-relativistic scattering: pulsating interfaces[END_REF][START_REF] Chen | Time-varying Doppler analysis of electromagnetic backscattering from rotating object[END_REF][START_REF] Chuang | Backscatter of a large rotating conducting cylinder of arbitrary crosssection[END_REF][START_REF] Harfoush | A numerical technique for analyzing electromagnetic wave scattering from moving surfaces in one and two dimensions[END_REF][START_REF] Liu | Analysis of phase noise influence on micro-Doppler feature extraction on vibrating target[END_REF][START_REF] Van Bladel | Electromagnetic fields in the presence of rotating bodies[END_REF]. Moreover, numerical approximation schemes based e.g. on FDTD [START_REF] Harfoush | A numerical technique for analyzing electromagnetic wave scattering from moving surfaces in one and two dimensions[END_REF][START_REF] Zheng | Analysis of micromotion characteristics from moving conical-shaped targets using the Lorentz-FDTD method[END_REF][START_REF] Zheng | Analysis of Doppler effect of moving conducting surfaces with Lorentz-FDTD method[END_REF] or fast integral equation solvers [START_REF] Zhang | Efficient analysis of scattering by multiple moving objects using a tailored MLFMA[END_REF] were developed. More theoretical contributions essentially linked to onedimensional moving boundary problems have also been considered [START_REF] Christov | On mechanical waves and Doppler shifts from moving boundaries[END_REF][START_REF] Fokas | Method for solving moving boundary value problems for linear evolution equations[END_REF]. To the best of the authors' knowledge, however, numerically solving the micro-Doppler PDE modeling problem has not been considered in the literature until now.

In the present paper, we propose a numerical approach to solve two-and three-dimensional micro-Doppler problems by the finite element method. The method is inspired by the ap-proach introduced in [START_REF] Gasperini | A frequency domain method for scattering problems with moving boundaries[END_REF] for the one-dimensional case, but has also major differences since the numerical method must be able to tackle efficiently much more complicated engineering problems (higher spatial dimension, complex geometries and materials). After stating the problem and introducing some notations in Section 2, we derive in Section 3 an approximate multi-harmonic Helmholtz-type weak coupled formulation well-adapted to the finite element method. To this end, in subsection 3.1, we rewrite the moving target problem written in weak form for a fixed domain by means of a mapping that can be computed by the finite element method. Since the resulting weak formulation in the fixed domain remains complicated to solve, we propose in subsection 3.2 an approximation of this formulation by a simplified weak form for small deformations. To justify this approximation, we derive some asymptotic error estimates of the bilinear forms according to the small deformation parameter setting. Then, the coupled multi-harmonic Helmholtz-type system is derived in subsection 3.3 from the approximate weak formulation. Considering the time Fourier series expansions of the geometrical quantities (e.g. the Jacobian tensor) for periodic movements and the unknown wavefield, we obtain an approximate truncated finite system of coupled weak form Helmholtz equations at different equally spaced frequencies. The finite element implementation is next developed in Section 4 for the coupled system. An efficient preconditioned solver that exploits the structure and properties of the resulting linear system is proposed in Section 5. Section 6 presents the application of the resulting method to a two-dimensional problem of the detection of a breathing baby placed on the rear seat of a car. Numerical results illustrate the behavior of the approach and its main features. Finally, we conclude our work in Section 7.

Problem statement and notations

Let t > 0 be the time variable and x := (x 1 , ..., x d ) T a point of R d (d = 2, 3). We assume that a wave is emitted at a source with localized fixed finite support Ω s ⊂ R d of boundary Γ s = ∂Ω s , and is scattered by a bounded obstacle Ω obst (t) ⊂ R d , with C s (s ≥ 2) boundary Γ(t) = Γ t = ∂Ω obst (t) smoothly deformed around an equilibrium position Γ 0 := Γ(t = 0), oscillating at a period T = 1/ν for a frequency ν . Let us remark here that some of the results stated in the paper do not need the periodicity condition of the deformation. We define Ω ext (t) as the d-dimensional unbounded homogeneous domain of propagation with boundary Γ s ∪ Γ(t), where we assume that Γ s ∩ Γ(t) = ∅ for all t > 0. Denoting by ∆ x := d i=1 ∂ 2 x i the spatial Laplace operator, the unknown total wave field u := u(x, t) is searched as the solution to the scalar wave equation

1 c 2 ∞ ∂ 2 t u -∆ x u = 0, (1) for
(x, t) ∈ Ω ext (t) × R + * (setting R + * := {t > 0})
, where c ∞ is the constant wave speed. (The extension to space-dependent wave speeds will be considered later, at the discrete level.) Here, we impose a T f -periodic time-dependent source term supported on Γ s , i.e. u(x, t)| Γ s = Ae iω f t , with real-valued amplitude A > 0, periodicity T f = 1/ν f for a frequency ν f , and angular frequency ω f := 2πν f > 0. In practice, ν f is supposed to be very large compared to ν . In what follows, we set a homogeneous Dirichlet boundary condition u(x, t) = 0 on Γ t , but the extension to other boundary conditions (e.g. Neumann or Fourier-Robin boundary condition) is relatively direct. In addition, the two initial conditions in Ω ext (t = 0) are set to u(x, 0) = 0 and ∂ t u(x, 0) = 0.

(

) 2
To apply later a numerical method, e.g. the finite element method, we need to bound the exterior domain Ω ext . To this end, we introduce a fixed fictitious smooth boundary Σ enclosing Ω s and Ω obst . The resulting finite computational domain Ω(t) = Ω t is then defined with boundary ∂Ω t = Σ ∪ Γ s ∪ Γ t . In the present paper, we simply apply the Sommerfeld absorbing boundary condition on Σ

1 c ∞ ∂ t u + ∂ n u = 0, (3)
with the normal derivative ∂ n u := ∇ x u • n, where n designates the exterior unit normal vector to Ω t and the spatial gradient operator is given by ∇ x := (∂ x 1 , ..., ∂ x d ) T . However, higherorder absorbing boundary conditions or PMLs could also be applied without any specific difficulty. In the sequel, we keep on denoting by u the solution to the wave equation ( 1) but set in Ω t × R + * , with the source term u(x, t) = Ae iω f t on Γ s , the homogeneous Dirichlet boundary condition u(x, t) = 0 on Γ t , the absorbing boundary condition (3) on Σ and with homogeneous initial data (2) in the compact set Ω(t = 0) = Ω 0 . We schematically illustrate the configuration in Figure 1. Let us remark that, in practical applications, e.g. Ω t being the interior of a car with shape Σ, impedance boundary conditions are usually prescribed on some parts of Σ to model transparent, absorbing or reflecting materials (see also Figure 2 in Section 6). Hence, the system does not have any sharp resonance and the useful frequency information is only related to the perturbation of the wave emitted by the radar Γ s , with high-frequency ν f ( ν ), by the local slowly oscillating smooth boundary Γ t , with frequency ν . Let us finally note that the limit situation of the fully reflecting domain is analyzed in details in the one-dimensional case in [START_REF] Gasperini | A frequency domain method for scattering problems with moving boundaries[END_REF].

n Γ s Σ Ω t Γ t

Approximate multi-harmonic Helmholtz-type weak coupled formulation

We now present an approximate multi-harmonic weak formulation based on coupled Helmholtz equations to solve the previous initial boundary-value problem. To this end, we first derive in subsection 3.1 an equivalent weak formulation of the initial problem in a fixed domain. An approximate variational formulation is then proposed in subsection 3.2 for small deformations. Next, in subsection 3.3, a truncated Fourier series expansion of the unknown scalar field and some other geometrical quantities are considered around ν f for uniformly spaced sample frequencies with step ν . This results in the solution of a system of coupled multi-harmonic Helmholtz-type equations written under a weak form. The finite element approximation of the formulation is then discussed in section 4.

Equivalent weak formulation in a fixed domain

The aim of this section is to write a weak form of the system in a fixed initial domain Ω 0 through a metric change parametrized by the time t > 0. For fixed t > 0, we designate by H m (Ω t ) the standard Sobolev space of solutions of finite energy on Ω t and introduce the Sobolev space

H m 0,A (Ω t ) := u ∈ H m (Ω t ) such that u| Γt = 0 and u| Γ s = Ae iω f t , (4)
for m = 1, 2. With our notations, we have in particular that, a t = 0,

H m 0,A (Ω 0 ) := {u ∈ H m (Ω 0 ) such that u| Γ 0 = 0 and u| Γ s = A} . ( 5
)
We also need the following Sobolev space

H m ω f t,A (Ω 0 ) := u ∈ H m (Ω 0 ) such that u| Γ 0 = 0 and u| Γ s = Ae iω f t . ( 6
)
All along the paper, we assume that the solution u := u(x, t) to the initial boundaryvalue problem stated in section 2 is an element of the space

C 0 ([0, T ], H 2 0,A (Ω t )), ∂ t u is in C 0 ([0, T ], H 1 (Ω t )) while ∂ 2 t u is C 0 ([0, T ], L 2 (Ω t ))
(see e.g. [START_REF] Lions | Non-Homogeneous Boundary Value Problems and Applications[END_REF] for regularity results of the solution of the wave equation for fixed domains).

According to these assumptions, the weak formulation in Ω t writes: find u(

•, t) ∈ H 1 0,A (Ω t ) such that ∀φ ∈ H 1 0,0 (Ω t ), Ωt ∇ x u • ∇ x φdΩ t + Ωt 1 c 2 ∞ ∂ 2 t uφdΩ t + Σ 1 c ∞ ∂ t uφ dΣ = 0, (7)
with the two initial conditions u(x, 0) = 0 and ∂ t u(x, 0) = 0. Here, we define φ as the complex conjugate of the function φ, a • b := d j=1 a j b j for two vectors a and b in C d , and the norm |a| of a is: |a| = √ a • a. We now consider a vector flow Φ t , for t ∈ R + , defined by

Φ t : Ω 0 → Ω t ⊂ R d x 0 → x = (Φ t1 (x 0 ), ..., Φ td (x 0 )) T = (x 1 , ..., x d ) T , ( 8
)
mapping the fixed initial domain Ω 0 to Ω t in such a way that we have

Φ t (Γ 0 ) = Γ t , Φ t (Σ) = Σ, Φ t (Γ s ) = Γ s and Φ 0 = Id, ( 9
)
where Id is the identity map from Ω 0 to Ω 0 . For fixed t ≥ 0, Φ t is assumed to be a C 2 -diffeomorphism with inverse map Φ -1 t : Ω t → Ω 0 . In addition, we suppose that the vector flow {Φ t } t>0 is of class C s in time, with s ≥ 2. Throughout the paper, x 0 = (x 01 , ..., x 0d ) T denotes a point of the fixed domain Ω 0 , image of x ∈ Ω t by the change of variables Φ -1 t , i.e. x 0 = Φ -1 t (x). Let u 0 be the unknown in Ω 0 after applying Φ -1 t and let us set u(x, t) = u 0 (Φ -1 t (x), t) = u 0 (x 0 , t) for t ∈ R + . We also consider the space-time vector X 0 := (Φ -1 t1 (x), ..., Φ -1 td (x), t) T = (x 01 , ..., x 0d , t) T ∈ Ω 0 × R + and define the two gradient operators ∇ X 0 := (∂ x 01 , ..., ∂ x 0d , ∂ t ) T and ∇ x 0 := (∂ x 01 , ..., ∂ x 0d ) T .

The time-dependent spatial Jacobian matrix J t of Φ -1 t at point x ∈ Ω t then writes

J t (x) = (∇ x Φ -1 t (x)) T .
By the change of variable Φ -1 t , we have the following relations:

∇ x u = J T t ∇ x 0 u 0 , ∂ t u = ∂ t X 0 • ∇ X 0 u 0 , ∂ 2 t u = ∂ 2 t X 0 • ∇ X 0 + (∂ t X 0 • ∇ X 0 ) 2 u 0 .
Let us explain a simple construction of Φ t , even if other possibilities could be alternatively considered. We introduce a smooth fixed surface Γ loc enclosing Γ t for all t ≥ 0, such that Γ loc ∩ Γ t = ∅. This defines a local neighborhood Ω loc t := Ω loc (t) of Γ t in Ω t . Let us set Ω loc 0 := Ω loc (t = 0). Therefore, we have: ∂Ω loc 0 = Γ 0 ∪ Γ loc , and the fixed complementary domain Ω fix := Ω 0 -Ω loc 0 is the part where the deformation is not active. Our aim is to model the impact of the local surface deformation to compute the volume coordinates changes, parametrized by t ≥ 0. By assumption, we know that the deformation δ t on Γ 0 is given by Φ t -Φ 0 ∈ C 2 (Γ 0 ) and by the knowledge of Γ t , which has C 2 regularity. In addition, since we assume that there is no effect of the deformation on Ω fix , we set: δ t = Φ t -Φ 0 = 0 on Ω fix , and in particular on Γ loc , for t ≥ 0. Now, we propose to build the volume deformation δ t := Φ t -Φ 0 through the local unique solution Φ t (a proof could be obtained by the Lax-Milgram theory for example at any time t > 0, accordingly to the smoothness properties of the data) to the following Dirichlet problem, parametrized by t > 0,

     -∆ x 0 Φ t = 0 in Ω loc 0 , Φ t = Id on Γ loc , Φ t = Φ t | Γ 0 on Γ 0 , ( 10
)
where Id is the identity operator over Γ loc and Φ t | Γ 0 is the given deformation on Γ 0 . The Laplacian is defined by ∆ x 0 := (∇ x 0 ) 2 . The Dirichlet problem [START_REF] Censor | Non-relativistic scattering: pulsating interfaces[END_REF] splits into d separated subproblems over each component Φ tj of Φ t , 1 ≤ j ≤ d. For a fixed time t, the solution to [START_REF] Censor | Non-relativistic scattering: pulsating interfaces[END_REF] provides Φ t which has components in C 2 (Ω loc 0 )∩C 0 (Ω loc 0 ), and by extension in C 2 (Ω 0 )∩C 0 (Ω 0 ). Now, we can compute the Jacobian J t of Φ -1 t as

J t = Jac(Φ t ) -1 . ( 11
)
For a fixed t > 0, the coefficients of J t are C 1 (Ω 0 ) functions defined as spatial derivatives of the C 2 mapping Φ t based on the Jacobian expression. Practically, a finite collection of the local deformations at some uniformly distributed discrete times will be computed through the solution of the boundary-value problem [START_REF] Censor | Non-relativistic scattering: pulsating interfaces[END_REF] by the finite element method. Through [START_REF] Chen | The Micro-Doppler Effect in Radar, 2nd Ed[END_REF], we obtain a way to compute J t by derivation (again numerically realized later by the finite element method) in Ω loc 0 . The tensor J t is trivially equal to the identity matrix in Ω fix . Let us recall that the Sobolev space H 1 ω f t,A (Ω 0 ) is given from definition (6) for m = 1. Then, the variational formulation [START_REF] Cammenga | High range resolution micro-Doppler analysis[END_REF] writes in the fixed domain Ω 0 , for fixed t > 0: find

u 0 (•, t) ∈ H 1 ω f t,A (Ω 0 ) such that, ∀φ ∈ H 1 0,0 (Ω 0 ), Ω 0 J t J T t ∇ x 0 u 0 • ∇ x 0 φ|J -1 t |dΩ 0 + Ω 0 1 c 2 ∞ ∂ 2 t X 0 • ∇ X 0 + (∂ t X 0 • ∇ X 0 ) 2 u 0 φ|J -1 t |dΩ 0 + Σ 1 c ∞ ∂ t u 0 φ dΣ = 0, (12)
since (see also expression ( 21))

(∂ t X 0 • ∇ X 0 )u 0 = ∂ t u 0 , on Σ, (13)
and where |J -1 t | denotes the determinant of J -1 t , which in particular is equal to 1 on Σ. Finally, we also have the two initial conditions u 0 (x 0 , 0) = 0 and ∂ t u 0 (x 0 , 0) = 0 in Ω 0 . From the assumptions on the smoothness of ∂ t u and

∂ 2 t u, we deduce that ∂ t u 0 is an element of H 1 (Ω 0 ) while ∂ 2 t u 0 is in L 2 (Ω 0 ). Remark 1.
Conversely to the methodology introduced in [START_REF] Gasperini | A frequency domain method for scattering problems with moving boundaries[END_REF], we first write here the weak formulation and then apply a change of coordinates so that the problem is set in a time independent domain with a metric that depends on the parameter t ≥ 0. In [START_REF] Gasperini | A frequency domain method for scattering problems with moving boundaries[END_REF], the change of coordinates was first applied and then the weak form was stated, which may provide a different formulation according to the algebraic calculations. However, the weak formulation of this previous approach is equivalent to [START_REF] Chen | Micro-Doppler effect in radar: phenomenon, model, and simulation study[END_REF] when applied to specific test functions of the form φ = |J -1 t |φ. Remark 2. Let us remark that we choose to model the deformation by [START_REF] Censor | Non-relativistic scattering: pulsating interfaces[END_REF]. More advanced elastic deformations could also be used. For the problem of interest here, this is nevertheless a sufficiently accurate model.

Approximate weak formulation for small deformations

For conciseness, let us set ∇ := ∇ x 0 since we are now working on the fixed domain Ω 0 . In practice, the exact formulation (12) could be used for a numerical computation. However, in the case of a small deformation of the moving interface, an accurate approximate formulation can be used to simplify the calculations. Indeed, a careful analysis developed in [START_REF] Gasperini | A frequency domain method for scattering problems with moving boundaries[END_REF] for the one-dimensional case, and formally extended to higher-dimensional problems, shows that the following simplified weak formulation [START_REF] Chen | High-resolution sparse representation of micro-Doppler signal in sparse fractional domain[END_REF] provides an accurate approximate solution v of u 0 that satisfies [START_REF] Chen | Micro-Doppler effect in radar: phenomenon, model, and simulation study[END_REF]. More precisely, the formulation writes: for a fixed time t > 0, find

v := v(•, t) ∈ H 1 ω f t,A (Ω 0 ) such that, for any test-function φ ∈ H 1 0,0 (Ω 0 ), Ω 0 J t J T t ∇v • ∇φ|J -1 t |dΩ 0 + Ω 0 1 c 2 ∞ ∂ 2 t vφ|J -1 t |dΩ 0 + Σ 1 c ∞ ∂ t vφdΣ = 0, (14)
with the two initial conditions v(x 0 , 0) = 0 and ∂ t v(x 0 , 0) = 0. The derivation of ( 14) is based on the following formal approximation

(∂ 2 t X 0 • ∇ X 0 + (∂ t X 0 • ∇ X 0 ) 2 )u 0 ≈ ∂ 2 t v in Ω 0 . ( 15
)
To give a more rigorous understanding on how these approximation lead to [START_REF] Chen | High-resolution sparse representation of micro-Doppler signal in sparse fractional domain[END_REF] as a correct approximation of ( 12), let us write the application Φ t with a small smooth perturbation 1 under the form

x = Φ t (x 0 ) = x 0 + δ t (x 0 ), (16)
where

(δ t ) t>0 ∈ C 2 (]0; T max ], C 2 (Ω 0 , R d ))
, where T max is a finite maximal time of computation. Then, for t ∈]0; T max ], the time derivatives of δ t up to the second-order are uniformly bounded, i.e. there exist a constant d such that, for any t ∈]0; T max ], for = 0, 1, 2, we have

d ,t := ∂ t δ t L ∞ d (Ω 0 ) ≤ d := sup t∈]0;Tmax] d ,t , (17)
with

L ∞ d (Ω 0 ) := (L ∞ (Ω 0 )) d .
Here, δ t may be oscillating but this is not yet mandatory. In the case of a periodic deformation, we can fix: T max = T . Now, since x 0 = Φ -1 t (x), we have from ( 16)

x 0 = Φ -1 t (x) = x -δ t (Φ -1 t (x)). ( 18
)
Based on this expression, we get

∂ t Φ -1 t (x) = -∂ t δ t (Φ -1 t (x)) =: -A(x 0 , t), (19)
and then

∂ t X 0 = ∂ t Φ -1 t (x) t = -A 1 . ( 20
)
Consequently, we have the following expansion in Ω 0

(∂ t X 0 • ∇ X 0 )u 0 = -A 1 • ∇ ∂ t u 0 = (∂ t -A • ∇)u 0 . ( 21
)
Now, we obtain

(∂ 2 t X 0 • ∇ X 0 + (∂ t X 0 • ∇ X 0 ) 2 )u 0 = -∂ t A 0 • ∇u 0 ∂ t u 0 + -A 1 • ∇ ∂ t -A 1 • ∇u 0 ∂ t u 0 . ( 22
)
By using [START_REF] Garcia-Rubia | Analysis of moving human micro-Doppler signature in forest environments[END_REF], we deduce that

(∂ 2 t X 0 • ∇ X 0 + (∂ t X 0 • ∇ X 0 ) 2 )u 0 = -∂ t A • ∇u 0 + (∂ t -A • ∇) 2 u 0 = ∂ 2 t u 0 -2 ∂ t (A • ∇u 0 ) + 2 A • ∇(A • ∇u 0 ). (23)
Now, to prove Proposition 1 which justifies the approximation ( 14), let us consider the second term in [START_REF] Chen | Micro-Doppler effect in radar: phenomenon, model, and simulation study[END_REF] and introduce the following bilinear forms from the decomposition ( 23)

a ,t (u 0 , φ) = a 0,t (u 0 , φ) + a 1 (u 0 , φ) + 2 a 2 (u 0 , φ), ( 24
)
with a ,t (u 0 , φ) :=

Ω 0 1 c 2 ∞ ∂ 2 t X 0 • ∇ X 0 + (∂ t X 0 • ∇ X 0 ) 2 u 0 φ|J -1 t |dΩ 0 , ( 25
)
a 0,t (u 0 , φ) := Ω 0 1 c 2 ∞ ∂ 2 t u 0 φ|J -1 t |dΩ 0 , ( 26
)
a 1,t (u 0 , φ) := -2 Ω 0 1 c 2 ∞ ∂ t (A • ∇u 0 )φ|J -1 t |dΩ 0 , ( 27
)
and

a 2,t (u 0 , φ) := Ω 0 1 c 2 ∞ A • ∇ (A • ∇u 0 ) φ|J -1 t |dΩ 0 . ( 28
)
Let us first consider a 1,t given by [START_REF] Karahasanovic | Radar-based detection of thoracoabdominal asynchrony during breathing using autocorrelation function analysis[END_REF]. Then, the Cauchy-Schwarz inequality implies that

|a 1,t (u 0 , φ)| ≤ 2γ t ||∂ t (A • ∇u 0 )|| L 2 (Ω 0 ) ||φ|| L 2 (Ω 0 ) , ( 29
)
where

γ t := |J -1 t | c 2 ∞ L ∞ (Ω 0 ) . ( 30
)
If we develop the right hand side derivative term we obtain

||∂ t (A • ∇u 0 )|| L 2 (Ω 0 ) = ||∂ t A • ∇u 0 + A • ∇∂ t u 0 || L 2 (Ω 0 ) ≤ ||∂ t A • ∇u 0 || L 2 (Ω 0 ) + ||A • ∇∂ t u 0 || L 2 (Ω 0 ) . ( 31
)
Since A has continuous first-order time derivative, then we can deduce that there exist a constant C 0,t which depends on

d 2,t = ||∂ t A|| L ∞ d (Ω 0 ) such that ||∂ t A • ∇u 0 || L 2 (Ω 0 ) ≤ C 0,t ||∇u 0 || L 2 d (Ω 0 ) ≤ C 0,t ||u 0 || H 1 (Ω 0 ) (32)
and a constant C 1,t depending on

d 1,t = ||A|| L ∞ d (Ω 0 ) such that ||A • ∇∂ t u 0 || L 2 (Ω 0 ) ≤ C 1,t ||∇∂ t u 0 || L 2 d (Ω 0 ) ≤ C 1,t ||∂ t u 0 || H 1 (Ω 0 ) . ( 33
)
If one collects ( 29)- [START_REF] Peng | A sinusoidal frequency modulation Fourier transform for radar-based vehicle vibration estimation[END_REF], we obtain

|a 1,t (u 0 , φ)| ≤ 2γ t (C 0,t ||u 0 || H 1 (Ω 0 ) + C 1,t ||∂ t u 0 || H 1 (Ω 0 ) )||φ|| L 2 (Ω 0 ) . ( 34
)
For the term a 2,t , let us first set:

A := c -2 ∞ |J -1 t |A, which leads to a 2,t (u 0 , φ) = Ω 0 ( Aφ) • ∇(A • ∇u 0 )dΩ 0 = ∂Ω 0 ( A • n)(A • ∇u 0 )φd(∂Ω 0 ) - Ω 0 ((∇ • A)φ + A • ∇φ)(A • ∇u 0 )dΩ 0 . ( 35
)
By assumption, we have:

∂Ω 0 = Σ ∪ Γ s ∪ Γ 0 . On Σ ∪ Γ s ,
since there is no deformation, A is identically equal to 0. In addition, the test-functions φ are in the space H 1 0,0 (Ω 0 ) and vanish on Γ 0 (as well as on Γ s ). Therefore, the boundary term in [START_REF] Saad | GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF] is equal to zero. For the second term, we have

|a 2,t (u 0 , φ)| = Ω 0 ((∇ • A)φ + A • ∇φ)(A • ∇u 0 )dΩ 0 ≤ Ω 0 ((∇ • A)φ)(A • ∇u 0 )dΩ 0 + Ω 0 ( A • ∇φ)(A • ∇u 0 )dΩ 0 . ( 36
)
For the first term on the right hand side of the above inequality, one gets

Ω 0 ((∇ • A)φ)(A • ∇u 0 )dΩ 0 ≤ C 3,t u 0 H 1 (Ω 0 ) φ L 2 (Ω 0 ) , ( 37
)
where C 3,t is a constant that depends on ∇• A and A. Concerning the second term, we obtain

Ω 0 ( A • ∇φ)(A • ∇u 0 )dΩ 0 ≤ C 4,t u 0 H 1 (Ω 0 ) φ H 1 (Ω 0 ) , ( 38
)
where the positive constant C 4,t is related to A and A. Combining ( 36)-( 38), we deduce that

|a 2,t (u 0 , φ)| ≤ C 2,t u 0 H 1 (Ω 0 ) φ H 1 (Ω 0 ) , ( 39
)
with C 2,t = C 3,t + C 4,t . Finally, this directly proves the following proposition.

Proposition 1. Let us define a ,t a 0,t , a 1,t and a 2,t , by [START_REF] Harfoush | A numerical technique for analyzing electromagnetic wave scattering from moving surfaces in one and two dimensions[END_REF], ( 26), ( 27) and [START_REF] Lahaye | Modern Solvers for Helmholtz Problems[END_REF], respectively, such that (24) holds. Then, for the solution

u 0 (•, t) ∈ H 2 0,A (Ω 0 ), ∂ t u 0 (•, t) ∈ H 1 (Ω 0 ), and
∂ 2 t u 0 (•, t) ∈ L 2 (Ω 0 )
, and for any test-function φ ∈ H 1 0,0 (Ω 0 ), we have the following error bounds, for t > 0,

|a ,t (u 0 , φ) -a 0,t (u 0 , φ)| ≤ |a 1,t (u 0 , φ)| + 2 |a 2,t (u 0 , φ)| ≤ (C 1 0,t u 0 H 1 (Ω 0 ) + C 1 1,t ∂ t u 0 L 2 (Ω 0 ) ) φ H 1 (Ω 0 ) + 2 C 2,t u 0 H 1 (Ω 0 ) φ H 1 (Ω 0 ) (40) and |a ,t (u 0 , φ) -a 0,t (u 0 , φ) -a 1,t (u 0 , φ)| ≤ 2 |a 2,t (u 0 , φ)| ≤ 2 C 2,t u 0 H 1 (Ω 0 ) φ H 1 (Ω 0 ) , ( 41
)
where we define:

C 1 0,t := 2γ t C 0,t , C 1 1,t := 2γ t C 1,t and C 2,t := γ t d 1,t (d 1,t + d 3,t ).
From this we conclude that solving the weak formulation ( 14) provides an approximation v of u 0 solution to [START_REF] Chen | Micro-Doppler effect in radar: phenomenon, model, and simulation study[END_REF] up to a bounded error of order O( ). Let us remark that the estimates (40) and (41) also hold uniformly in time since one can take the maximum of C 1 0,t , C 1 1,t and C 2,t over [0, T max ]. Furthemore, an improved first-order formulation can be obtained thanks to the bilinear form a 0,t + a 1,t up to an error O( 2 ). More precisely, one could replace the bilinear form in ( 14) by

(a 0,t + a 1,t )(v, φ) = 0, ( 42
) with a 1,t (v, φ) = 2 Ω 0 1 c 2 ∞ ∂ t ∂ t Φ -1 t • ∇v φ|J -1 t |dΩ 0 .
However, in the present paper, we restrict our analysis to [START_REF] Chen | High-resolution sparse representation of micro-Doppler signal in sparse fractional domain[END_REF]. In the following, we provide a suitable formulation to [START_REF] Chen | High-resolution sparse representation of micro-Doppler signal in sparse fractional domain[END_REF] when the motion of the boundary is T -periodic in time.

Coupled multi-harmonic Helmholtz-type system

We now assume that the parametrization of Γ t is periodic with periodicity T . Then, the smooth Jacobian J t is also a T -periodic tensor of class C s in time, with ν ν f . Consequently, both J T t , J t J T t , |J -1 t | and |J -1 t |J t J T t are also time-dependent function/tensors with time periodicity T with C s regularity according to t > 0 and C 1 smoothness in Ω 0 (since the quantities can be expressed as sums and products of T -periodic smooth functions). In particular, we have the Fourier series expansions

|J -1 t |J t J T t =: C(x 0 , t) = n∈Z C n (x 0 )e inω t and |J -1 t | =: c(x 0 , t) = n∈Z c n (x 0 )e inω t , ( 43
)
for x 0 ∈ Ω 0 , where

C n (x 0 ) := 1 T T /2 -T /2 |J -1 t |J t J T t e -inω t dt and c n (x 0 ) := 1 T T /2 -T /2 |J -1 t |e -inω t dt. ( 44
)
The tensor fields C n satisfy C n = C -n since |J -1 t |J t J T t is a real-valued matrix. In addition, we have c n = c -n because |J -1 t | is real-valued. Finally, these two quantities have a C 1 (Ω 0 ) regularity.

In view of a numerical implementation of the coefficients C n and c n , we discretize the T -periodic motion of Γ t at 2N + 1 equally spaced sampling times t k = k∆t , 0 ≤ k ≤ 2N , with time step ∆t := T /(2N + 1). This is a natural choice since it avoids any aliasing effect and provides a uniform discretization over one period T . Moreover, to satisfy the Nyquist-Shannon sampling theorem, (2∆t ) -1 has to be larger than the maximal frequency arising in |J -1 t |J t J T t and |J -1 t |, represented as (43). Thus, we obtain a sequence {Γ k } 0≤k≤2N of boundaries sampling Γ t at the discrete times t k , 0 ≤ k ≤ 2N , described through the maps Φ k = Φ t k . Therefore, based on the numerical solution of the 2N + 1 local boundary value problems [START_REF] Censor | Non-relativistic scattering: pulsating interfaces[END_REF], we can compute both C k := C(t k ) and c k := c(t k ) from J k := J(t k ) evaluated at the discrete times t k , 0 ≤ k ≤ 2N . Consequently, the approximations of C and c (defined by ( 43)) are given by the following truncated tensor C N and function c N

C N (x 0 , t) := n∈N C n (x 0 )e inω t and c N (x 0 , t) := n∈N c n (x 0 )e inω t , (45)
which are all well-defined at any point x 0 in Ω 0 . Here, we set

N := {n ∈ Z; -N ≤ n ≤ N }.
Since we assumed that Φ t is a T -periodic diffeomorphism of class C s in time, with s ≥ 2, then one gets for example for the function c (but a similar property arises for the tensor C), for j = 0, ..., s, that it is such that

∂ j t c(x 0 , t) := n∈Z 2iπn T j c n (x 0 )e inω t ,
for x 0 ∈ Ω 0 . This implies that c n L ∞ (Ω 0 ) = o(|n| -s ), justifying hence that (45) are suitable approximations of (43) if N is chosen large enough. In particular, if the mapping is C ∞ , then one gets rapidly decaying coefficients since the o(|n| -s ) asymptotic holds for any s > 0 (see section 6 for an illustrative numerical example). Concretely, in a finite element context, the tensor fields C n and functions c n are obtained by using some FFTs, locally on each element of the finite element triangulation, then leading to piecewise constant tensors and functions (see e.g. sections 4 and 6.1). Now and similarly to [START_REF] Gasperini | A frequency domain method for scattering problems with moving boundaries[END_REF], we develop the solution v of ( 14) as a Fourier series expansion with complex Fourier coefficients a j v(x 0 , t) = j∈Z a j (x 0 )e iω j t = e iω f t j∈Z a j e ijω t , (46) for x 0 ∈ Ω 0 , where ω j = ω f + jω = 2πν j , with ν j := ν f + jν . Here, we assume that a j ∈ H 2 0,0 (Ω 0 ), meaning that a j (x 0 ) = 0 on Γ s for j ∈ Z -{0}, and a 0 ∈ H 2 0,A (Ω 0 ), i.e. a 0 (x 0 ) = A on Γ s . In view of a numerical procedure, we have to truncate the Fourier series expansion of the unknown v. To this end, let us introduce the finite sets of integers J := {j ∈ Z ; -J ≤ j ≤ J} and J * := J -{0}, for a fixed truncation index J ≥ 0. We denote by w J the Fourier series expansion (46) truncated with J := 2J + 1 terms, i.e.

w J (x 0 , t) = j∈J a j (x 0 )e iω j t , ( 47
)
with a j ∈ H 2 0,0 (Ω 0 ), for j ∈ J * , and a 0 ∈ H 2 0,A (Ω 0 ). The parameter J must be carefully chosen to include the significant contributions to the solution [START_REF] Gasperini | A frequency domain method for scattering problems with moving boundaries[END_REF]. Let us remark that the Fourier transform of (47) in the sense of distributions is given as a finite sum of Dirac distributions δ ν j located at the discrete frequencies {ν j } j∈J , modulated by the functions a j , i.e.

w J = j∈J a j δ ν j , ( 48
)
where the Fourier transform of a given function f is defined by

f (ξ) := ∞ -∞ f (t)e -2iπξ dt.
Setting w J (x 0 , 0) = 0 and ∂ t w J (x 0 , 0) = 0 as initial conditions, we obtain, by derivation of (47),

∂ t w J (x 0 , t) = j∈J iω j a j e iω j t and ∂ 2 t w J (x 0 , t) = - j∈J ω 2 j a j e iω j t . ( 49
)
Considering the approximations (45) and (49), we have

c N ∂ 2 t w J = -e iω f t n∈N c n e inω t j∈J ω 2 j a j e ijω t = -e iω f t J+N m=-(J+N ) J j=-J -N ≤m-j≤N c m-j ω 2 j a j e imω t = - J+N m=-(J+N ) J j=-J -N ≤m-j≤N c m-j ω 2 j a j e iωmt , ( 50
)
since the following identity holds

N n=-N α n J j=-J β j = J+N m=-(J+N ) J j=-J -N ≤m-j≤N α m-j β j . ( 51
)
Similarly, one gets in Ω 0

C N ∇w J = e iω f t n∈N C n e inω t j∈J ∇a j e ijω t = J+N m=-(J+N ) J j=-J -N ≤m-j≤N C m-j ∇a j e iωmt . ( 52
)
Based on these truncated Fourier series expansions and from the previous calculations, we obtain the following approximation of the weak form ( 14): find a j ∈ H 1 0,0 (Ω 0 ), for j ∈ J * , and a 0 ∈ H 1 0,A (Ω 0 ) such that, for any φ ∈ H 1 0,0 (Ω 0 ),

J+N m=-(J+N ) J j=-J -N ≤m-j≤N Ω 0 C m-j ∇a j • ∇φdΩ 0 -ω 2 j Ω 0 c m-j c 2 ∞ a j φdΩ 0 +δ m-j Σ iω j c ∞ a j φdΣ e iωmt = 0, (53)
where δ m-j is the delta Kronecker symbol, i.e. δ m-j = 1 if m = j, and zero otherwise. Finally, by projecting (53) onto the Fourier modes -J ≤ n ≤ J in the outer sum over m, we obtain the square set of (2J + 1) equations with (2J + 1) unknowns {a} j∈J , for n ∈ J ,

J j=-J -N ≤n-j≤N Ω 0 C n-j ∇a j • ∇φdΩ 0 -ω 2 j Ω 0 c n-j c 2 ∞ a j φdΩ 0 + δ n-j Σ iω j c ∞ a j φdΣ = 0. ( 54
)
A first possibility consists in assuming that we fix N = 2J to get access to the (4J + 1) Fourier coefficients (C n ) -2J≤n≤2J and (c n ) -2J≤n≤2J . Then, (54) writes, for n ∈ J ,

J j=-J -2J≤n-j≤2J Ω 0 C n-j ∇a j • ∇φdΩ 0 -ω 2 j Ω 0 c n-j c 2 ∞ a j φdΩ 0 + δ n-j Σ iω j c ∞ a j φdΣ = 0. (55)
This first choice leads to a global weak form with (2J + 1) × (2J + 1) blocks involving sparse matrices. Let us now assume that we compute the (2J + 1) first Fourier coefficients in (45), i.e. N = J. This has the advantage of considering less modal coupling between the geometry and the unknown expansions when computing their products. As a consequence, the resulting (2J + 1) × (2J + 1) system is given by a 2J + 1 band limited coupled system, for n ∈ J , J j=-J -J≤n-j≤J

Ω 0 C n-j ∇a j • ∇φdΩ 0 -ω 2 j Ω 0 c n-j c 2 ∞ a j φdΩ 0 + δ n-j Σ iω j c ∞ a j φdΣ = 0, (56)
involving 3J 2 + 3J + 1 blocks, instead of 4J 2 + 4J + 1 for the first approach.

More generally, if one only retains 2N + 1 terms in the Fourier expansions (45), leading to indices -N ≤ n-j ≤ N in (54), then the global system is still block banded, with band length 2N -1, and requires n (N,J) := (4N + 2)J + (1 + N -N 2 ) blocks. From these remarks, we will say that the Fourier approximation is of order (N, J) when we consider an approximation of the unknown w J given by (47) with 2J + 1 Fourier terms, and 2N + 1 (0 ≤ N ≤ 2J) Fourier coefficients for the geometrical quantities, based on (45). The (2J + 1) × (2J + 1) system to be solved is then given by, for n ∈ J , J j=-J -N ≤n-j≤N

Ω 0 C n-j ∇a j • ∇φdΩ 0 -ω 2 j Ω 0 c n-j c 2 ∞ a j φdΩ 0 + δ n-j Σ iω j c ∞ a j φdΣ = 0. ( 57
)
When solving (57), the approximation of order (N, J) of w J solution to (47) is defined by

w (N,J) (x 0 , t) = j∈J a (N,J) j (x 0 )e iω j t , ( 58
)
where the a (N,J) j

, j ∈ J , are the 2J + 1 unknown functions associated to (57). Finally, let us remark that the case of inhomogeneous media, i.e. when the speed is space dependent (c ∞ = c ∞ (x 0 )), is considered in Sections 4 and 6. To this end, we directly adapt the previous formulation by considering a piecewise constant finite element approximation of the velocity, according to the local properties of the media.

Finite element approximation

Let us consider a fixed covering Ω 0,h = Ω h (0) of Ω 0 with n 0,h triangular/tetrahedral elements in 2D/3D. We also set n loc 0,h as the number of elements involved in the discrete domain Ω loc 0,h

interpolating Ω loc 0 . The basis functions of the P 1 finite element approximation are denoted by {φ j } 1≤j≤n dof , where n dof is the number of degrees of freedom for the Lagrange finite element approximation of H 1 0,0 (Ω 0 ). For j ∈ J , we designate by a (N,J) j,h the P 1 finite element approximation of a (N,J) j and define the interpolated field w

(N,J) h in Ω 0,h approximating (58) w (N,J) h (x 0,h , t) = j∈J a (N,J) j,h (x 0,h )e iω j t . ( 59
)
Higher-order finite element methods (or alternative spatial discretization methods) could be used to increase the accuracy as well as to minimize the dispersion/pollution error related to high frequency time-harmonic wave problems [START_REF] Ihlenburg | Finite Element Analysis of Acoustic Scattering[END_REF]. The representation in the finite element basis of the unknown wave field w (N,J) h built from the (N, J)-approximation is given by the vector w (N,J) = (a (N,J) -J , ..., a (N,J) 0 , ..., a (N,J) J ) T , with vector components a

(N,J) j ∈ C n dof , j ∈ J .
Therefore, w (N,J) is an element of C n J dof , with n J dof = (2J + 1)n dof . Let us briefly explain how to build the linear system related to (57). To compute J k,h , which is the P 0 finite element approximation of J k on Ω 0,h , we solve [START_REF] Censor | Non-relativistic scattering: pulsating interfaces[END_REF] under its weak form. This leads to the solution of (2N + 1) variational problems with unknown discrete ddimensional deformations Φ k,h , 0 ≤ k ≤ 2N , in the small local fixed domain Ω loc 0,h , where only the boundary condition on Γ 0,h changes at index k. This means that the matrix associated with the linear system resulting from the discretization of (10) remains unchanged while the contribution of the boundary conditions modifies the right-hand side of the linear system to solve. To this end, a Choleski factorization of the sparse symmetrical positive definite matrix (here, the stiffness matrix for the Laplacian) can be used before the time loop over k while only 2d sparse triangular linear systems can next be efficiently solved for each new right-hand side at step k, with k = 0, ..., 2N . Based on P 1 finite elements, one gets the finite sequence of discretized deformations Φ k,h , 0 ≤ k ≤ 2N . Then J k,h is obtained by a simple numerical derivation through the discrete gradients and based on [START_REF] Chen | The Micro-Doppler Effect in Radar, 2nd Ed[END_REF], corresponding to P 0 approximations of the sequence of Jacobians, and piecewise constant tensors/functions C k h and c k h . As a consequence, on each element in Ω loc 0,h , a Fourier series expansion can be computed locally to get C n,h and c n,h , -N ≤ n ≤ N , from the sampling C k h and c k h , 0 ≤ k ≤ 2N . The main computational cost in this procedure is related to the FFT computations locally on each element of the domain Ω loc 0,h to obtain C n,h and c n,h , -N ≤ n ≤ N , which means O(n loc 0,h N log(N )) operations.

Let us now introduce the generalized finite element stiffness matrix K n and mass matrix M n , for -N ≤ n ≤ N , respectively given by the n dof × n dof matrices

(K n ) ı, := Ω 0,h C n,h ∇φ ı • ∇φ  dΩ 0,h , (M n ) ı, := Ω 0,h c n,h c 2 ∞ φ ı φ  dΩ 0,h , 1 ≤ ı,  ≤ n dof .
These two symmetrical complex-valued sparse matrices are such that:

K n = K -n , M n = M -n (if c ∞ is real-valued).
In addition, we need the real-valued highly sparse generalized mass matrix M Σ on the fixed interpolated fictitious boundary Σ h (but considered as a n dof × n dof matrix)

(M Σ ) ı, := Σ h 1 c ∞ φ ı φ  dΣ h , 1 ≤ ı,  ≤ n dof ,
where c ∞ can be assumed to be defined constant by element, i.e. as the P 0 approximation of the space-dependent velocity c ∞ . Let us remark that considering a non-homogeneous medium inside the domain (i.e. different physical materials) leads to non constant velocity c ∞ according to the space variable, and therefore to additional coupled terms due to the metric change. For the P 0 approximation, both approaches are identical. We now introduce some block banded matrix notations. We define S (N,J) as the block Toeplitz matrix, with (2J + 1) × (2J + 1) blocks and bandwidth equal to 2N + 1, based on a collection of n dof × n dof sparse complex-valued elementary matrices S n , -N ≤ n ≤ N , such that

S (N,J) =              S 0 • • • S -N 0 • • • 0 . . . .
0 • • • 0 S N • • • S 0              . ( 60
)
Let us consider the global stiffness matrix K (N,J) associated with K n , n ∈ N , which is a block Toeplitz matrix, with bandwidth 2N + 1. In addition, it is generated by the N + 1 elementary matrices K n , n = 0, ..., N , since K -n = K n . If the storage of one elementary stiffness matrix is γn dof , with γ > 0 a small constant that depends on the mesh interconnections, then storing K (N,J) needs γ(N + 1)n dof coefficients. Similarly, we can introduce the global mass matrix M (N,J) , based on the elementary blocks M n , n = 0, ..., N , with similar properties as K (N,J) . To simplify the explanations, we assume that storing M n also needs γn dof complex values. Now, let us consider the two diagonal matrices: ω (J) = diag(ω j I n dof ) j∈J and ω 2,(J) = diag(ω 2 j I n dof ) j∈J , where I n dof designates the unitary matrix of size n dof × n dof . We also need the matrix M (J) Σ , where M Σ is repeated (2J + 1) times on the diagonal. Both its construction cost and memory storage are negligible compared with K (N,J) and M (N,J) since M Σ is only related to the boundary nodes on Σ h . Finally, the finite element solution of (57) implies that the boundary conditions generate some contributions to the right hand side vector b (N,J) ∈ C n J dof . To conclude, the frequency domain finite element method yields the block linear system

A (N,J) ω (J) w (N,J) = b (N,J) , (61) setting A
(N,J)

ω (J) := K (N,J) -M (N,J) ω 2,(J) + iω (J) M (J) Σ , (62) where A
(N,J)

ω (J) := (A (N,J) ω (J) ,j,n ) (j,n)∈J
×J is a block banded matrix but is not block Toeplitz. In addition, we have the block hermitian property: A (N,J)

ω (J) ,j,n = A (N,J) ω (J) ,n,j , for (j, n) ∈ J * × J * .

Efficient solution of the linear system

Let us now focus on the efficient numerical solution of the linear system (61). The storage of the matrix A (N,J) ω (J) is given by the matrices K (N,J) , M (N,J) and M (J) Σ , and the bruteforce application of a direct solver is not computationally tractable in the high-frequency regime even for moderate values of N and J, as each single high-frequency Helmholtz-type problem is already known to be extremely challenging, requiring the use of advanced solvers [START_REF] Ernst | Why it is difficult to solve Helmholtz problems with classical iterative methods[END_REF][START_REF] Lahaye | Modern Solvers for Helmholtz Problems[END_REF]. Here, our aim is to show that solving (61) scales as (2J + 1) calls to robust single high-frequency Helmholtz-type solvers through an iterative process based on a preconditioned GMRES [START_REF] Saad | Iterative methods for sparse linear systems[END_REF][START_REF] Saad | GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF]. One of the key ingredients of GMRES is to compute the Matrix-Vector Products (MVPs) y ← A (N,J) ω (J) x, where x and y are some complex-valued vectors in C n J dof . In our situation, this can be easily achieved block-wise, based on the Toeplitz storage. Since A (N,J)

ω (J) has n (N,J) := (4N + 2)J + (1 + N -N 2 )
blocks, then each global MVP has a cost 2γn dof n (N,J) . In practice, N is often smaller than J (see section 6). Therefore, asymptotically, the cost of a MVP is about 8γn dof N J, which means that it scales linearly with J. Another crucial point of the acceleration of the GMRES is to use a preconditioner. Here, we build a left preconditioner P such that PA (N,J) ω (J) ≈ I := I n dof (2J+1) to improve the convergence rate of the Krylov method. A robust preconditioner P for GMRES must yield a clustering of the eigenvalues of PA (N,J) ω (J) around (1, 0) in the complex plane [START_REF] Saad | Iterative methods for sparse linear systems[END_REF]. Since we have a Toeplitz storage of A (N,J) ω (J) , we cannot easily build an approximation of (A (N,J) ω (J) ) -1 . Let us instead consider the splitting: A (N,J)

ω (J) = A (0,J) ω (J) + ∆A (N,J) ω (J) , where A (0,J) ω (J) := diag(A (0,J) ω j ) j∈J is the static part of A (N,J) ω (J) , i.e. A (0,J) ω j = K 0 -ω 2 j M 0 + iω j M Σ , for j ∈ J .
Let us assume that we have access to an efficient solver for each linear system associated to the j-th block matrix A (0,J) ω j . Then, we could consider a preconditioner P based on the "inversion" of the 2J + 1 diagonal blocks, which means solving an associated linear system for each j ∈ J . Nevertheless, since j is changing, the linear systems are j-dependent. Let us recall that: ω j = ω f +jω , with ω f ω . This means that we can consider the high-frequency approximation ω j ≈ ω f = ω 0 with associated block diagonal preconditioner P := diag((A (0,J) ω 0 ) -1 ) j∈J . Since we are solving the preconditioned linear system PA (N,J)

ω (J) w (N,J) = Pb (N,J) , ( 63
)
we then have PA (N,J)

ω (J) = I + R (N,J) ω (J)
, where R (N,J) ω (J) is expected to be a small perturbation of the identity matrix I. The application of the preconditioner P at each GMRES iteration requires to solve the linear system

A (0,J) ω 0 x j = b j , (64)
for 2J + 1 right hand sides b j , j ∈ J . Therefore, the cost for the application of the preconditioner in each iteration is equal to (2J + 1) times the cost of solving (64). In practice, a sparse LU factorization of P is computed before entering into the GMRES iteration loop after reordering the unknowns by a Cuthill-MacKee algorithm. This then leads to an exact factorization A (0,J) ω 0 = L ω 0 U ω 0 and block LU preconditioner P. Let us remark that incomplete LU factorizations A (0,J) ω 0 L ω 0 U ω 0 (or possibly other single-frequency more advanced preconditioners) can alternatively be used [START_REF] Saad | Iterative methods for sparse linear systems[END_REF]. Applying such LU preconditioners P requires the solution to 2(2J + 1) decoupled upper/lower sparse complex-valued triangular systems. Remark 3. Let us remark that, under the high-frequency assumption ω j ≈ ω 0 = ω f , j ∈ J , for ω ω f , the linear system ( 61)-( 62) can be approximated by

A (N,J) ω (0) w (N,J) ω (0) = b (N,J) , (65) setting A
(N,J)

ω (0) := K (N,J) -ω 2 f M (N,J) + iω f M (J) Σ . The matrix A (N,J) ω (0) has (2J + 1) × (2J + 1
) blocks, each one being ω j independent (unlike A (N,J) ω (J) ), leading to a block Toeplitz structure. More precisely, the matrix writes

A (N,J) ω (0) =              T 0 • • • T -N 0 • • • 0 . . . . . . . . . . . . . . . T N . . . . . . 0 0 . . . . . . T -N . . . . . . . . . . . . 0 • • • 0 T N • • • T 0              , with T 0 = A (0,J) ω (0) , T n = K n -ω 2 f M n and T -n = T n , 0 ≤ n ≤ N . Therefore, one gets A (N,J) ω (J) = A (N,J) ω (0) + ω ω f c 2 ∞ R (N,J) ,
where R (N,J) can be expressed in terms of mass matrices that can then be bounded in norms.

A direct computation shows that we have w (N,J)

ω (0) -w (N,J) w (N,J) ω (0) ≤ ω ω f c 2 ∞ cond(A (N,J) ω (J) ) R (N,J) A (N,J) ω (J) (66)
for a given matrix norm A of a matrix A and where the condition number is defined by cond(A) = A × A -1 . Storing A (N,J) ω (0) still requires γ(N + 1)n dof complex coefficients. Solving (65) can be expected to provide a suitable solution to the problem, under the assumption

ω ω f c 2 ∞
1, by using an algorithm based on fast block Toeplitz solvers [START_REF] Alonso | An efficient parallel algorithm to solve block-Toeplitz systems[END_REF]. In addition, one may also use A (N,J) ω (0) to build a banded limited preconditioner by dropping some off-diagonal blocks T n for solving (62), then generalizing P = (A (0,J)

ω (0) ) -1 .
6 Application to a model problem

Description of the problem

We now apply the proposed method to a simplified two-dimensional test-case motivated by the high-frequency radar detection of a breathing baby installed in a child seat located on the rear seat of a car. Let us consider the initial geometry at t = 0 depicted in Figure 2, which represents the two-dimensional longitudinal section of the interior of a car. The computational domain Ω 0 is delimited by the shape of the car, including the surface of the body of the baby in his seat, facing the road. We assume that the red parts of the car, called Σ 0 , are Perfectly Electrical Conductors (PEC), i.e. we impose a homogeneous Dirichlet boundary condition. In addition, the belly of the baby Γ 0 (black part on Figure 2) is also considered as a PEC assuming e.g. that the baby is wearing a perfectly conducting jacket which reflects the wave field created by an antenna modeled by an extruded circle Γ s placed near the roof of the car, above the rear seat. This antenna emits a transverse magnetic electromagnetic field of amplitude A > 0 and high-frequency ν f . The belly of the baby is the part Γ t of the boundary that moves in the modeling. At t = 0, Γ 0 is a curved segment with two endpoints p 1 and p 2 . Finally, the rest of the surface of the car, called Σ 1 (blue part on Figure 2), is considered to be non-reflecting as a rough approximation, assuming that it is constituted e.g. of glass or synthetic materials. This is modelled by the zeroth-order absorbing boundary condition [START_REF] Brooks | Complex-valued neural networks for fully-temporal micro-Doppler classification[END_REF]. Therefore, within our notations, Ω 0 is the domain with boundary Γ 0 ∪ Γ s ∪ Σ 0 ∪ Σ 1 , while

∂Ω t := Γ t ∪ Γ s ∪ Σ 0 ∪ Σ 1 .
In Ω 0 , we consider that, except in the seats (grey parts), we have some air. Therefore, we model the propagation of the wave in Ω air 0 by the velocity c air of the light in the air. In the seats Ω seat 0 , we have specific dielectric materials where the velocity takes different values, depending on the refraction index. More precisely, Ω seat 0 models the seat foam (gray volume in Figure 2) with complex-valued constant velocity c seat and refractive coefficient (with the air) η seat

c seat = c air η seat , with η seat := µ seat r ε seat r + i θ seat ν f ε air , ( 67
)
where ε seat r , µ seat r , θ seat and ε air are the relative permittivity and permeability of the material, its conductivity, and the permittivity of the vacuum, respectively. Here, √ • denotes the principal determination of the square-root with branch-cut along the negative real-axis. Let us remark that these physical parameters are chosen only according to the emission frequency ν f . Therefore, we can define the velocity c ∞ as piecewise constant in our model, i.e. c ∞ := c air,seat in Ω air,seat 0 , and Ω 0 = Ω air 0 ∪ Ω seat 0 . Now, we need to introduce the contour Γ loc to define the domain of computation of the metric change. Here, we choose a fixed polygonal open curve included inside Ω 0 with endpoints p 1 and p 2 (dashed black line on Figure 2). Therefore, this delimits the domain Ω loc 0 with boundary Γ 0 ∪ Γ loc . When the baby is breathing, Γ 0 is deformed to Γ t , but the two points p 1 and p 2 stay fixed. In the paper, we consider a small deformation of size 1 described by a simple T -periodic sine mapping deforming the initial configuration Γ 0 into Γ t in the normal direction. More precisely, for x 0 ∈ Γ 0 and t ≥ 0, we define

x = Φ t (x 0 ) = x 0 + δ t (x), (68)
with the deformation

δ t (x 0 ) = sin π l(p 1 x 0 ) l(Γ 0 ) sin(ω t)n(x 0 ), (69)
assuming ω ω f . Trivially, one gets δ t (x 0 ) ≤ 1 for all x 0 ∈ Γ 0 and t ≥ 0. Here, l(p 1 x 0 ) is the arclength over Γ 0 from p 1 to x 0 , and l(Γ 0 ) denotes the length of Γ 0 . The unit normal vector n(x 0 ) at x 0 ∈ Γ 0 is outwardly directed to Ω 0 . More complex breathing models are described in [START_REF] Dias Da | A theoretical investigation of the detection of vital signs in presence of car vibrations and radar-based passenger classification[END_REF][START_REF] Karahasanovic | Radar-based detection of thoracoabdominal asynchrony during breathing using autocorrelation function analysis[END_REF]. In our framework, considering general movements is easy since only the moving of the surface finite element interpolation nodes on Γ 0 is needed. In addition, we have: Ω fix 0 := Ω 0 -Ω loc 0 , which is the part of Ω 0 where the deformation is not active. Finally, the problem under consideration is given by: find the total wave field u in Ω t , t > 0, such that

                         1 c 2 ∞ ∂ 2 t u -∆ x u = 0 in Ω t , 1 c ∞ ∂ t u + ∂ n u = 0 on Σ 1 , u = 0 on Σ 0 ∪ Γ t , u = Ae iω f t at Γ s , u(x, 0) = 0 in Ω 0 , ∂ t u(x, 0) = 0 in Ω 0 . ( 70
)
We now consider the discrete times {t k } 0≤k≤2N , where N ≥ 1. Since Φ t is given following (68)-(69), we can alternatively compute simultaneously the sequence of discrete deformations δ k in the layer Ω loc 0 as the solutions to

         -∆δ k = 0 in Ω loc 0 , δ k = 0 on Γ loc , δ k = sin π l(p 1 x 0 ) l(Γ 0 ) sin( 2πk 2N + 1 )n(x 0 ) for x 0 ∈ Γ 0 , (71)
for 0 ≤ k ≤ 2N . The solution δ k,h to (71) is realized by a P 1 finite element method for the corresponding weak form in Ω loc 0,h (with n loc 0,h triangles), providing hence Φ k,h (x 0,h ) = x 0,h + δ k,h (x h ). A second step allows to calculate the P 0 approximation of J k,h by a weak derivation. The solution to the linear systems is obtained by computing once the Choleski factorization of the stiffness matrix related to the problem (71) and then by a forwardbackward solution for each of the 2(2N + 1) right hand sides.

Let us remark that the size L of the layer Ω loc 0,h along the normal direction must be sufficiently large compared to so that the deformation does not flatten too much the triangulation related to Ω loc 0,h . At the same time, the mesh in Ω 0,h must also be sufficiently refined to correctly capture the high oscillations of the solution due to ω f . The coefficients C n and c n , -N ≤ n ≤ N , are finally evaluated as piecewise constant on each triangle from the Jacobian J k,h , 0 ≤ k ≤ 2N , using [START_REF] Chen | The Micro-Doppler Effect in Radar, 2nd Ed[END_REF] and FFTs.

Numerical study of the method

We consider the following physical parameters corresponding to realistic applications of radar detection of car occupants [START_REF] Anthony | Dielectric properties of polymers at microwave frequencies: a review[END_REF][START_REF] Elhawil | A quasi-optical free-space method for dielectric constant characterization of polymer materials in mm-wave band[END_REF]: c ∞ = 3 × 10 8 in Ω air 0 , µ r = 1, ε r = 1.2, ε air = 8.85 × 10 -12 and θ = 2.7 × 10 -3 . These values yield c seat = 2.7 × 10 8 -1.2 × 10 6 i according to (67). The emission frequency and amplitude, and the breathing motion frequency are respectively ν f = 3 × 10 9 , A = 1 and ν = 1. The layer Ω loc 0,h is chosen with thickness L = 0.1. The tuning parameters , N and J are varying parameters of the study.

Let us first illustrate through an example the construction of the function c n,h and tensor C n,h from the metric computation, which is a pre-processing step of the method. We fix = 0.01 and N = 4 (see below for this choice), which means that the time interval [0, T ] = [0, 1] is discretized by 2N +1 = 9 equally spaced points. We report on Figure 3a the triangulation Ω 0,h (gray mesh) zoomed around the rear seat. For the sake of clarity, we consider a relatively large element size (triangle edge length) of 6 × 10 -3 (leading to n 0,h = 115702). We superimpose on Ω 0,h (colored circles) the nodes of the initial triangulation mapped by Φ k,h , for k = 6 (i.e. t k = 0.67), resulting from the finite element approximation to (71). From this calculation, we can deduce the piecewise constant function c n,h and tensor C n,h . We represent in Figures 3b-3e the amplitude of (C n,h ) 1,1 , n = 0, • • • , 3. In particular, we observe that (C 0,h ) 1,1 = 1 and (C n,h ) 1,1 = 0 (n = 1, 2, 3) outside Ω loc k,h . One of the main points of the method is to correctly adjust the value of N in (45) to have an accurate representation of the moving geometry but also to avoid enlarging the bandwidth of the linear system (61) more than necessary. Now, for a given calculation, the mesh size must be fixed to capture the small oscillations of the wave field related to the high frequency signal (with usual wavelength in the range 0.01-0.1). For ν f = 3 × 10 9 , the element size in the triangulation Ω 0,h is fixed at 7 × 10 -4 , which leads to n 0,h = 8001827 triangular elements and n dof = 4000778 nodes. For n ∈ N, we set

m n := c n,h L ∞ (Ω loc 0,h ) c 0,h L ∞ (Ω loc 0,h ) , M n := max 1≤ı,≤2 (C n ) ı, L ∞ (Ω loc 0,h ) max 1≤ı,≤2 (C 0 ) ı, L ∞ (Ω loc 0,h )
.

Then, we fix N as the smallest positive integer such that: m N ≤ τ and M N ≤ τ , where τ is a small tolerance parameter. Usually, because of the expressions of c n and C n , and from numerical simulations (see also Figure 4 (Left)), the most important point is to check the criterion on M n . We report M n and m n , for 0 ≤ n ≤ 8, in Figure 4 and see that N = 4 (respectively N = 5) for τ = 10 -2 and = 0.01 (respectively = 0.02). In practice, we define an initial value for N and adjust it to fulfill the above criterion. The associated local cost per element in Ω loc 0,h is however low, for a global cost of about N log(N )n loc 0,h . This pre-processing step is important since it impacts the block structure of the matrix A (N,J) ω (J) . In addition, we observe that the ratios M n (as well as m n ) are rapidly decaying and behave as e -κ n when We now consider the solution of the linear system (61) and its preconditioned version (63) by using the GMRES without restart. The computations were performed using Matlab on a laptop with Intel Core i7-10510U CPU at 2.30 GHz. The LU decomposition A (0,8) ω 0 = L ω 0 U ω 0 defining P is used to solve (64) for j ∈ J , as explained in Section 5. The resulting upper and lower n dof × n dof triangular matrices are highly sparse for nnz/n 2 dof = 2.2 × 10 -3 %. We report in Figure 7 the residual history for solving (61) and (63), for 5 ≤ J ≤ 8, with = 0.01 and = 0.02. The parameter settings are the same as before. Without the preconditioner, the convergence is very slow (the maximum number of iterations is fixed to 50) and may even breakdown for some configurations. When using the preconditioner P, we obtain a very good convergence rate which seems to be independent of J. In addition, the slope of the convergence rate increases with which is natural since the boundary movement perturbates more the static configuration, which penalizes the preconditioner efficiency. In Table 2, we can observe that the CPU time for the iterative solution scales linearly with the parameter J and is faster than the direct Matlab solver for = 0.01. We expect that this difference will be much more important for higher frequencies or/and for three-dimensional problems.

Ω loc k,h Γ k,h Ω 0,h Γ loc h • • • • • •

Conclusion and perspectives

This paper proposes a general approach for computing the solution to two-and threedimensional scattering problems by moving boundaries. To this end, the initial boundary value problem is rewritten in a fixed domain thanks to a well-adapted mapping that can be computed by a finite element method. Next, approximations are introduced according to the small deformation of the boundary, and some geometrical quantities as well as the solution are expanded in terms of truncated Fourier series. This leads to the solution of a finite coupled system of Helmholtz-type equations, that can be solved by the finite element method. The numerical implementation of the method is fully detailed on a simplified model example arising from the automotive industry. This shows the potentiality of the approach for solving moving boundary problems for wave scattering, with applications to Doppler effects. Since the method is now validated, further improvements remain to be investigated. Among them, let us mention the use of high-order finite element methods and the development of robust solvers for very high frequency three-dimensional problems. Finally, all the material presented here can be extended to the three-dimensional Maxwell's equations.

Fig. 1 .

 1 Fig. 1. Example of a two-dimensional (d = 2) computational domain Ω(t) with source term on the boundary Γ s , with moving scattering surface Γ t and fictitious boundary Σ.

1 Fig. 2 .

 12 Fig. 2. Schematic description of the 2D model problem.

( a )Fig. 3 .

 a3 Fig. 3. Computation of the local deformation and piecewise constant tensor C n,h for 0 ≤ n ≤ 3.

Fig. 4 .

 4 Fig. 4. M n , m n and α (N,8) j , for 0 ≤ n, j ≤ 8 (threshold τ = 10 -2 ).

Fig. 5 .

 5 Fig. 5. Amplitudes |a (4,8) j,h (x 0 )|, 0 ≤ j ≤ 3, in Ω 0,h . The parameters are: c ∞ = 3 × 10 8 , ν f = 3 × 10 9 , ε r = 1.2, θ = 2.7 × 10 -3 , A = 1, ν = 1, = 0.01 and L = 0.1.

Fig. 6 .

 6 Fig. 6. Amplitudes |a (4,8) j,h (x 0 )|, 0 ≤ j ≤ 3, in Ω 0,h . The parameters are: c ∞ = 3×10 8 , ν f = 1.5×10 10 , ε r = 1.2, θ = 2.7 × 10 -3 , A = 1, ν = 1, = 0.01 and L = 0.1.

Fig. 7 .

 7 Fig. 7. Convergence history of GMRES for solving (61) (black) and (63) with increasing J, for = 0.01 (red) and = 0.02 (blue).

Table 2 :

 2 CPU time for solving (61) with a direct solver and (63) with the preconditioned GMRES.

		.01	1124	1276	1450	1719
		0.02	1709	1893	2304	2627
	Direct Matlab solver	0.01 0.02	1667 1680	1924 2028	2371 2489	2743 2858

Acknowledgements. The authors thank the support of the Luxembourg National Research Fund (FNR) (2017-1 PPP 11608832). This research was funded in part through the ARC grant for Concerted Research Actions (ARC WAVES 15/19-03), financed by the Wallonia-Brussels Federation of Belgium.

n > n 0 , for some constants n 0 ∈ N and κ > 0, i.e. max 1≤ı,≤2

which is due to the fact that we use a smooth mapping.

For N fixed, since the computation of c n,h and C n,h is spatially localized near the moving boundary as corrections of the static case (n = 0) and their maximal amplitude decays with n, we can expect that the blocks (A (N,J) ω (J) ) j,n of A (N,J) ω (J) , (j, n) ∈ J * × J * , have a sparser profile with smaller coefficients when n increases. We report on Table 1 the number of nonzero elements (nnz) and the Frobenius norm • F of these blocks, for j = 0 and n = 0, ..., 4, fixing J = 8 and for = 0.01 and = 0.02 (N is obtained thanks to the above truncation criterion). We can see that the number of elements and the Frobenius norm of the blocks get smaller as k grows. In addition, a larger value of implies that N must be taken larger and that both the number of elements and the norm of the blocks is slightly larger.

Let us now analyze the selection of the truncation parameter J ≥ N for the unknown expansion (59). In Figure 4, we report α

) , for 0 ≤ j ≤ J = 8, and = 0.01 (N = 4) and = 0.02 (N = 5). We observe that the behaviour of α (N,J) j is similar to the behavior of M n in both situations, i.e. we have rapidly decaying coefficients according to j. We recommend in practice to take J a little bit larger than N , e.g. J = N + 3. While it would be very interesting to have a priori and a posteriori error estimators for the expansion of the solution, this is not trivial as it depends on the global regularity of the procedure, and in particular of the mapping. For completeness we report the amplitudes of the four first functions |a (4,8) j,h (x 0 )|, 0 ≤ j ≤ 3, in Ω 0,h when = 0.01, for ν f = 3 × 10 9 on Figure 5 and ν f = 1.5 × 10 10 on Figure 6. Let us remark that when the radar frequency increases, the first modes contribute more to the scattering problem [START_REF] Gasperini | A frequency domain method for scattering problems with moving boundaries[END_REF], leading to the Doppler effect.