Thomas Brion

Pascal Fossat

Mohamed Ichchou

Olivier Bareille

Abdel-Malek Zine

Christophe Droz

Identification of partial differential equations in structural mechanics theory through k-space analysis and design

Keywords: Inverse method, k-space, Wave propagation, Partial Differential Equations, Stiffened plates

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

Introduction

Composite structures and meta-materials are used in a wide range of applications. Those structures can show exotic performances unobtainable with classic materials [START_REF] Del | Dynamic problems for metamaterials: review of existing models and ideas for further research[END_REF][START_REF] Mukhopadhyay | Stochastic mechanics of metamaterials[END_REF]. They can even be built such that they provide original properties in other fields like acoustic, dielectric or optic properties. Even smart structures are now developed [START_REF] Montgomery | Recent advances in additive manufacturing of active mechanical metamaterials[END_REF]. Composite materials are made of the assembly of at least two non miscible components and are therefore heterogeneous structures. Meta-materials are structures mostly made of spatially periodic patterns of assemblies of different materials. The contrast between the components is often responsible for the original properties. The mechanical behavior of this type of structure can be dominated by one of its components, depending on the test or the observation conducted. For these reasons, this type of structure needs non classical formulations. Due to their complexity, their modeling can be using computationally expensive models. And because of the contrasts, a not high enough resolution model can lead to the loss of some contributions to the behavior. However, the behavior can also be modelled using equivalent formulations [START_REF] Timoshenko | Theory of plates and shells[END_REF] or non classical analytic formulations built by homogenisation process [START_REF] Panasenko | Multi-scale modelling for structures and composites[END_REF][START_REF] Hans | Dynamics of discrete framed structures: a unified homogenized description[END_REF][START_REF] Boutin | Microstructural effects in elastic composites[END_REF]. These analytic approaches require the knowledge of equivalent parameters that can't always be directly calculated. One can mention stiffness, Young modulus in a direction and apparent mass among others.

In structural mechanics, one needs the PDE describing the motion in order to push forward the design and the analysis. The issue is very important and is still considered in the literature with regards to new materials and meta-structures specifically. Precisely, a great number of equation of motion are often used in structural engineering. These PDE rely on one-dimensional, two-dimensional or even tree-dimensional behavior. Models like Euler-Bernoulli, Timoshenko, Kirchhoff are probably the most established and employed ones. All correspond to the simplification of the 3D elastodynamic theory for specific kinematic behaviors. Such simplification could lead to discrepancies when the frequency increases or for specific materials. The issue is tricky for newly developed concepts of metamaterials and metastructures where the heterogeneiticy scale becomes a key factor. Two solutions could then be investigated. The first one makes use of multiscale homogenisation theories to get the PDE representing the media. The second solution is investigated in this paper: it makes use of an inverse technique to identify the main parameter of the PDE.

A lot of work has already been done in the field of identification and characterisation of materials and structures from different measurements. The Oberst beam method [START_REF] Renault | Characterization of elastic parameters of acoustical porous materials from beam bending vibrations[END_REF][START_REF] Wojtowicki | New approach for the measurement of damping properties of materials using the oberst beam[END_REF] can be cited for beams to regain parameter values but only at resonance frequencies. Some techniques to calculate or estimate elastic parameters in the case of porous materials are presented in [START_REF] Jaouen | Elastic and damping characterizations of acoustical porous materials: Available experimental methods and applications to a melamine foam[END_REF], often the Young modulus or the damping properties. From a kinematic full field measurement, five methods are presented in [START_REF] Avril | Overview of identification methods of mechanical parameters based on full-field measurements[END_REF]: the finite element model updating method, the constitutive equation gap method, the virtual fields method, the equilibrium gap method and the reciprocity gap method. In [START_REF] Shirota | Adjoint method for the problem of coefficient identification in linear elastic wave equation[END_REF], a spatial dependent Young modulus is identified based on the knowledge of the surface traction and using an adjoint method. Moreover some approaches make use of waves propagation to extract information from different type of structures. [START_REF] Every | Determination of the elastic constants of anisotropic solids from acoustic-wave group-velocity measurements[END_REF] makes use of group velocity measurements to estimate elastic constants in anisotropic structures. The elastic properties are identified from a transitory displacement field known on a set of points in [START_REF] Liu | Material characterization of fgm plates using elastic waves and an inverse procedure[END_REF][START_REF] Liu | An inverse procedure for determination of material constants of composite laminates using elastic waves[END_REF]. [START_REF] Liu | Material characterization of fgm plates using elastic waves and an inverse procedure[END_REF] deals with functionally graded materials and [START_REF] Liu | An inverse procedure for determination of material constants of composite laminates using elastic waves[END_REF] with laminate composites but improves the search of optimal values using a genetic algorithm. [START_REF] Deng | Characterization of elastic anisotropy of a solid plate using nonlinear lamb wave approach[END_REF] proposes a nonlinear Lamb wave approach to characterise solid plates. In [START_REF] Sale | Semi-analytical formulation for the guided waves-based reconstruction of elastic moduli[END_REF] an identification of elastic constants is proposed, based on phase velocity of guided waves using a semi analytic finite element formulation. And in the case of layered composites, in [START_REF] Chronopoulos | Accurate structural identification for layered composite structures, through a wave and finite element scheme[END_REF], the parameters are identified comparing time of flight from measurements and from wave finite element resolutions.

The wave propagation characteristics provide a lot of information about the behavior of a structure [START_REF] Leissa | Vibration of plates[END_REF]. They are already used in a lot of other applications: the time of flight of ultrasonic waves for non destructive testing [START_REF] Chien | An acousto-ultrasonic nde technique for monitoring material anisotropy[END_REF] or the propagation of Lamb waves for damage detection in laminate composite materials [START_REF] Kessler | Damage detection in composite materials using lamb wave methods[END_REF][START_REF] Su | Guided lamb waves for identification of damage in composite structures: A review[END_REF].

In particular the wavenumbers can be studied. Wavenumbers are complex quantities which characterise the propagation of the waves in a structure. They can be interpreted as an angular spatial frequency and they can be linked with the frequency and the geometric and mechanical parameters of the structure with the dispersion relation. A k-space can provide a lot of information on the vibration behavior of the structure. In the case of 2-dimensional structures like plate or panel, the wavenumber is a propagation direction-dependent function. The set of the wavenumbers in all the direction is the k-space. Some methods exist to get the k-space from a measured displacement field. As it deals with the spatial frequency, one of the oldest tool to study them is the Fourier transform [START_REF] Heckbert | Fourier transforms and the fast fourier transform (fft) algorithm[END_REF]. It can be used to get the wavenumber but is limited to real wavenumbers, so the imaginary part is neglected. It is ideal for an infinite domain. Moreover, in the discrete case, the spatial sample step has to be below enough the studied wavelength (Nyquist-Shannon Theorem).

That is why other methods should be used like the Inhomogeneous Wave Correlation method (IWC) or INverse COnvolution MEthod (INCOME). IWC, [START_REF] Berthaut | K-space identification of apparent structural behaviour[END_REF][START_REF] Ichchou | Identification of effective sandwich structural properties via an inverse wave approach[END_REF] is based on the maximisation of a correlation function using an inhomogeneous wave. In two dimension, this wave is parameterized by the angle and the wavenumber at this angle. INCOME, [START_REF] Boukadia | An inverse convolution method for wavenumber extraction (income): Formulations and applications[END_REF][START_REF] Boukadia | Vibroacoustics of periodic structures: model order reduction, characterization, optimization[END_REF] is an other method to get the k-space from a displacement field. The aim of the method is to find a convolution kernel such that the convolution between it ant the discrete field is minimized. The form of the kernel is determined by the dimension of the field and the hypothesis on the equation. After finding the coefficients of the convolution kernel, the wavenumber in a direction come from the resolution of a polynomial equation made of the coefficients the convolution kernel.

The purpose of this work is to use the propagation characteristics to identify PDE coefficients from a solution. The process presented is adapted to find information on wave propagation type equation from solution fields on 2 dimension domains. It is based on the comparison between analytic formulations and kspaces extracted from solution fields at different frequencies. It could be very useful to determine parameters of equivalent models like homogenisation models and therefore simplify the resolution of those problems in the case of meta-structure. In particular, this process can provide frequency dependent values for the equivalent parameters. The second main objective is the resolution of all types of linear PDE, not only wave based equations, using the finite element solver of COMSOL and the mathematics module. The dimensions of the problem to solve are not greater than 3. A lot of boundary conditions can be defined, like Dirichlet conditions or Neumann conditions for 2 nd -order equations but a lot of others combinations are possible as well. The limitation in order to the 2 nd -order is overridden by replacing the high order equation by system of equations with a maximal of two derivations in each equation. The idea is near writing a 1 dimension differential equation under matrix form.

This paper presents a method to use the finite element software solver to get the solution of linear PDE of any order, with a lot of different possibilities for boundary conditions and loading. It also present an identification process for PDE coefficient of 2 dimensional waves propagation equations, based on the comparison between analytic formulation of and reference k-spaces extracted from solutions field at different frequencies. Both the implementation of the PDE resolution and the identification method have been tested on plate and membrane in isotropic cases and non isotropic cases including a multi scale case.

The paper is structured as follows. Section 2 presents a method to implement on COMSOL the resolution of linear partial differential equations of any order with possibly unusual terms, on potentially complicated domain. Then, section 3 presents a wavenumbers based identification method for PDE coefficient which is then tested on 2 nd -and 4 th -order cases and in particular on a homogenised model of a bi-directional ribbed plate [START_REF] Boutin | Dynamics of ribbed plates with inner resonance: Analytical homogenized models and experimental validation[END_REF].

Partial differential equations implementation in on finite element software

This section shows how to implement the resolution of PDE in COMSOL software [START_REF] Multiphysics | Introduction to comsol multiphysics®[END_REF] using the mathematics module. In this module, the general form of PDE Interfaces is used 1 . The PDE can be in dimension 1 to 3 and the resolution can be done on complicated domains. The PDE will have the following form in 1D (1), 2D (2) or 3D [START_REF] Montgomery | Recent advances in additive manufacturing of active mechanical metamaterials[END_REF]. The coefficients (a m ), (a l,m-l ) or (a l,n-l,m-n ) and a 0 are the equation parameters. ω is the pulsation such that these equations can be studied at different frequencies. N is the total order of the equations. m, n and l are the sum index such that all the derivative and cross derivatives of order not greater than N are considered.

N m=1 a m ∂ m u ∂x m + ω 2 a 0 u = 0 (1) N m=1 m l=0 a l,m-l ∂ m u ∂x l ∂y m-l + ω 2 a 0 u = 0 (2) N m=1 m n=0 n l=0 a l,n-l,m-n ∂ m u ∂x l ∂y n-l ∂z m-n + ω 2 a 0 u = 0 (3)
Moreover the coefficients of the equations can be dependent on the spatial variables or the frequency. In addition to the equation the boundary conditions have to be defined. They can be related to the unknown, its derivatives or a combination of them.

Equations formulation

The equations (1), ( 2) and (3) cannot be directly implemented in the software and have to be put in the following form.

∇.Σ = f ( 4
)
Where Σ is a conservative flux term and f a source term that have to be defined. In their definition, the unknown and only its first derivative can appear. So with the presence of the divergence operator ∇., only 2 nd -order equation seems to be implementable. To increase the order of the equations the equation has to be replaced by a system of equations. (1) can be written:

∂ ∂x         u ∂u ∂x . . . ∂ N -1 u ∂x N -1         =        0 1 0 . . . 0 . . . 0 1 ω 2 a 0 a N - a 1 a N • • • - a N -1 a N                u ∂u ∂x . . . ∂ N -1 u ∂x N -1         (5)
1 Some usual equations are already implemented but the general form allows more possibilities

3

(5) is a 1 st -order problem. Because second derivatives can be implemented, not all the intermediate deriva-

tives, but only the even ones need to be defined.

∂ n u ∂x n is the notation of the defined function which is the n th derivative of u. Without loss of generality, if the order N of the equation is even, equation ( 1) is equivalent to the following equation:

∂ ∂x ∇.                    N/2 p=1 a 2p ∂ ∂x ∂ 2p-2 u ∂x 2p-2 + a 2p-1 ∂ 2p-2 u ∂x 2p-2 ∂ ∂x [u] ∂ ∂x ∂ 2 u ∂x 2       N/2 p=1 2p l=2 a l,2p-l ∂ ∂x ∂ 2p-2 u ∂x l-2 ∂y 2p-l + a l-1,2p-l ∂ 2p-2 u ∂x l-2 ∂y 2p-l ∂ ∂x [u] ∂ ∂y [u] 0       N/2 p=1 a 0,2p ∂ ∂y ∂ 2p-2 u ∂y 2p-2 + a 1,2p-1 ∂ ∂x ∂ 2p-2 u ∂y 2p-2 + a 0,2p-1 ∂ 2p-2 u ∂y 2p-2 0 0 ∂ ∂y [u] . . . ∂ ∂y ∂ N -2 u ∂y N -2                      (8) f = -ω 2 a 0 u ∂ 2 u ∂x 2 ∂ 2 u ∂x∂y ∂ 2 u ∂y 2 . . . ∂ N u ∂y N T (9)
This way the equation is defined without any derivation greater than 2 directly written.

Software implementation

The physics PDE general form (g) from the mathematics module allows to work with an equation of the type:

e a ∂ 2 u ∂t 2 + d a ∂u ∂t + ∇.Σ = f (10)
In the studied cases, only the part ∇.Σ = f is used2 . There is a COMSOL syntax for the first derivatives:

ux for the first x derivative of u, uy for its first y derivative and uz its z derivative. So Σ and f have to be defined in PDE general form like explained in the previous subsection. To complete the problem, the boundary conditions and the source terms are needed.

Boundary conditions and source terms

In COMSOL, the Dirichlet boundary conditions allow to add a constraint on a function at a boundary of the domain. The constraint can be a constant or a function of the space and the frequency. However the functions on which there is a constraint have to be defined in the system of equation ∇.Σ = f like the even derivatives. The boundary can be a point, a line or a surface.

Source terms

The source terms, depending on their forms and dimensions can be defined directly in the equation definition, using the boundary conditions or by adding a point source, an edge source or a source to the model. When the dimension of the load is the same as the one of the domain, it can either be put in the formulation of the equation, in ∇.Σ = f or in a source added to the model.

Frequency and spatial dependent parameter

As a reminder it is possible to implement complicated equations, with frequency and spatial dependent coefficients. These terms are implemented in the variables. A lot of functions are pre implemented and can be used. These terms can be used in the formulation of the equation but also in the source term or in the boundary condition. "freq" is the reserved variable used for the frequency and "x", "y" or "z" are the reserved ones for the spatial coordinates. At the end, a frequency domain study is conducted. A summary of the implementation procedure of PDE resolution in the COMSOL software is presented in Algorithm 1.

Examples of implementations

Three examples are presented for equations of different orders in different dimensions. The first one is a beam equation which is a 4 th -order equation in 1D to illustrate the implementation of a more than 2 nd -order equation. The second is a membrane equation to illustrate the implementation for a 2D equation. And the last one is a isotropic plate equation for both 2D and more than 2 nd -order equation. For each case the implementation is presented and then the results are compared with analytic solutions.

Euler Bernoulli beam equation

The first example is the Euler Bernoulli beam equation with a simply supported condition at the boundary x = 0 and a unitary displacement in a simply supported way at x = L. It is a 4 th -order equation in one dimension.

EI ∂ 4 u ∂x 4 = ω 2 Λu (11) u(x = 0) = 0; u(x = L) = 1; ∂ 2 u ∂x 2 (x = 0) = 0 and ∂ 2 u ∂x 2 (x = L) = 0 (12)
With E the Young modulus, I = bh 3 12 the bending inertia and Λ = ρbh the linear mass. The implementation of this equation requires the definition of the second derivative of u or of the bending moment : M = EI ∂ 2 u ∂x 2 . A way to implement this problem is the following one:

Σ =      ∂M ∂x EI ∂u ∂x      f =   ω 2 Λu M  
The boundary conditions are set using "Dirichlet Boundary Condition", at the boundary x = 0: u = 0, M = 0 and at x = L, u = 1, M = 0. The solution can be compared with the analytic solution of this

u analytic (x) = sin (kx) 2 sin (kL) + sinh (kx) 2 sinh (kL) ( 13
)
The implementation is done with EI = 1 kg m 3 /s 2 , Λ = 1 kg/m and L = 1 m. The numerical resolution of the problem [START_REF] Avril | Overview of identification methods of mechanical parameters based on full-field measurements[END_REF] with the boundary conditions of ( 12) is implemented. The results of the resolution are 130 presented and compared to the analytic solution [START_REF] Every | Determination of the elastic constants of anisotropic solids from acoustic-wave group-velocity measurements[END_REF] in term of displacement function of the position figure 1b and in decibel function of the frequency figure 1c. It can be seen on this figure that the numerical solution agrees perfectly with the analytic solution. The equation and its resolution have been well implemented.

Isotropic membrane equation

The next problem is the Membrane deflection equation on a square domain of side L with a fixed boundary condition and a unitary punctual load at the center. The implementation of the membrane equation is easier even if it is a 2 dimensional problem because it is only a 2 nd -order equation. So no additional function has to be defined. The equation is the following one.

N ∆u + Λω 2 u = 0 ( 14
)
N is the surface tension and Λ the surface mass. The boundaries are considered fixed, so u = 0 on the 135 boundaries. The implementation on COMSOL using general PDE is direct.

Σ = N ∂u ∂x N ∂u ∂y f = ω 2 Λu
The boundary conditions need the use of "Dirichlet conditions" to set u to 0 at the boundaries. The load is put using a "punctual load" at the center of the domain (x = L/2, y = L/2). The solution can be compared with the analytic solution from the modal decomposition.

u analytic (x, y) = +∞ n=0 +∞ m=0 4(-1) n+m sin (2n + 1)πx L sin (2m + 1)πy L L 2 Λω 2 -π 2 N ((2n + 1) 2 + (2m + 1) 2 ) ( 15
)
The implementation is done with unitary values N = 1N/m and Λ = 1kg/m 2 . The resolution has been implemented in the software. Figures 3b and3a show the solution fields from the numerical resolution and from the analytic formula [START_REF] Liu | An inverse procedure for determination of material constants of composite laminates using elastic waves[END_REF] and figure 3c their comparison on a line. The two fields seem to be almost the same and the comparison of the sectional view shows a perfect agreement.

140

Moreover the comparison of the solutions at a point function of the frequency is presented figure 3d and shows a very good agreement. The finite element implementation solves the equation well.

Isotropic plate equation

Then the Kirchhoff Love equation for the flexion in an isotropic plate is implemented. The plate is supposed square with a length L. The boundary conditions are simply supported and a unitary punctual load is applied at the center. It is a 4 th -order equation. This case combines the dimension 2 with a higher than 2 nd -order.

D∆ 2 u = ω 2 Λu ( 16
)
Σ =            D ∂ ∂x ∂ 2 u ∂x 2 + 2 ∂ ∂x ∂ 2 u ∂y 2 D ∂ ∂y ∂ 2 u ∂y 2 D ∂u ∂x 0 0 D ∂u ∂y            f =           ω 2 Λu ∂ 2 u ∂x 2 ∂ 2 u ∂y 2          
A punctual load is set to 1 on the point at the center of the plate. And two "Dirichlet boundary conditions" are defined: u = 0, ∂ 2 u ∂x 2 = 0 at {x = 0}, {x = L} and u = 0, ∂ 2 u ∂y 2 = 0 at {y = 0}, {y = L}. The solution obtained can be compared to the analytic solution from modal decomposition.

u analytic (x, y) = +∞ n=0 +∞ m=0 4(-1) n+m sin (2n + 1)πx L sin (2m + 1)πy L L 2 π 4 D L 4 ((2n + 1) 2 + (2m + 1) 2 ) 2 -ω 2 Λ ( 18
)
The implementation is done with D = 1 kg m 2 /s 2 , Λ = 1 kg/m 2 and L = 1m. A good agreement between the numerical solutions from the software and the analytic ones can be seen figures 5a, 5b, 5c and 5d. In all the previous cases it seems that the implementations on COMSOL software give results very near the analytic solutions, for the form of the solution field and for the frequency response functions. The implementation of the resolution of PDE using its finite element solver seems to work well. Now that the implementation process has been validated, numerical resolutions can be conducted to provide solution fields of 2D equations. This got solution fields are used in the next section to test the identification process.

Partial differential equations identification by k-space analysis

A method for the identification of parameter of partial differential equation on 2 dimensions domains based on wavenumber analysis is presented. The idea is to use the k-spaces extracted from displacement fields as references and to compare them with the analytic formulation of the dispersion relation taken from the model. With this approach, a frequency dependent value of a parameter of the model can be obtained.

Types of equation

This method works for propagation equations in the frequency domain on two dimensional domains. The equations are contained in the form [START_REF] Mukhopadhyay | Stochastic mechanics of metamaterials[END_REF]. For this type of equation, the analytical relation dispersion can be found as a root of the following polynomial equation, function of the direction in the domain.

N m=1 (ik θ ) m m l=0 a l,m-l cos l θ sin m-l θ = -ω 2 a 0 (19)

Methodology

From given solution fields, the objective is to retrieve to a coefficient of the equation, which can have physical meaning. The method presented here uses wavenumber analysis to find an approximation of a coefficient when the other are known. The very first step is to determine the form of the equation, the order, the known and unknown coefficients and which of them are null or linked. From the form of the equation, the analytic dispersion relation or polynomial equation can be defined [START_REF] Leissa | Vibration of plates[END_REF]. It is function of the coefficients of the equation, in particular the missing one, but also the frequency and the direction.

Then wavenumbers have to be extracted from the given solution fields using k-space method. The methods used are the Inhomogeneous wave correlation method (IWC) [START_REF] Ichchou | Identification of effective sandwich structural properties via an inverse wave approach[END_REF][START_REF] Berthaut | K-space identification of apparent structural behaviour[END_REF] and the INverse COnvolution MEthod (INCOME) [START_REF] Boukadia | An inverse convolution method for wavenumber extraction (income): Formulations and applications[END_REF][START_REF] Boukadia | Vibroacoustics of periodic structures: model order reduction, characterization, optimization[END_REF] they are described in Appendix A. For each solution field, a set of n θ wavenumbers are calculated in different directions θ i : k ref (f, θ i ). Now the identification can be done. The value of the unknown parameter, named p, is searched such that the wavenumbers calculated from the analytic dispersion relation are as near as possible to the extracted one from the solution field, at each direction. That is why at a given frequency f, the identification of the value of the wanted parameter of the model use the resolution of the following least square minimization problem on all the calculated directions [START_REF] Chien | An acousto-ultrasonic nde technique for monitoring material anisotropy[END_REF].

J(p) = 1 n θ n θ i=1 |k ref (f, θ i ) -k model (f, θ i , p)| 2 (20)
This way the coefficient is identified. In the following results, the search of the minimum has been done by evaluation but other methods of extremum research could be use.

The method is tested on equations of 2 nd -and 4 th -order. First the equations are solved for different frequencies with a finite element solver as explained in section 2. Then the method presented previously is used to find one of the parameters of the equation from the solution.

2 nd -order equations

The identification of parameter begin with solution fields from the following type of 2 nd -order equation.

a x ∂ 2 u ∂x 2 + a y ∂ 2 u ∂y 2 + a xy ∂ 2 u ∂x∂y + ω 2 a 0 u = 0 ( 21
)
In this case, the analytic dispersion relation is immediate:

k 2 θ = ω 2 a 0 a x cos 2 θ + a y sin 2 θ + a xy cos θ sin θ ( 22
)
This equation is solved on a square domain with a side of 1m and with Dirichlet boundary conditions, or fixed boundary conditions for the interpretation on the behavior of an elastic membrane. A unitary punctual load is put at the center of the domain. The frequency are between 0 and 50 Hz, ω is 2π the frequency. The parameter a x , a y , a xy and a 0 are tested with different configurations. The grid uses 400 by 400 points.

In order to use INCOME on fields solution of this type of equation, the convolution kernel has to be at least of the following form. The choice of the kernel comes from [START_REF] Boukadia | Vibroacoustics of periodic structures: model order reduction, characterization, optimization[END_REF]. Three configurations for the coefficients a x , a y , a xy and a 0 are tested.

Isotropic membrane

The first case is the one of an isotropic membrane. a x = 1, a y = 1, a xy = 0 and a 0 = 1 The equation is symmetric in the x and y direction, a x = a y = a. The equation ( 21) can be written:

a∆u + a 0 ω 2 u = 0 ( 24
)
The solution at the frequency 5 Hz has been presented figure 3a. Using INCOME with the kernel (23) the wavenumbers are calculated in different directions but because of the values chosen, the wavenumbers are expected to be independent from the direction.

In this simple case, figure 6a the calculated wavenumbers are very near the analytic ones, except at low frequencies. It is because, at low frequencies the length of the domain is almost the same as the wavelength. For each frequency, the wavenumbers are calculated in 100 angles between 0 and π. The identification parameter method is then used to regain a parameter when the others are known. The values are compared with the expected ones entered in the COMSOL implementation. It is done for the parameters a 0 and a x . It can be seen in figures 7a and 7b that when the frequency is higher than 10 Hz, the results for the identification of a x and those for the identification of a 0 are near the expected values.

The bad values at low frequencies can come from the size of the domain which is not high enough compared to the wavelength. In order to verify that the equation has been solved at low frequency on a larger domain: L=5m in place of 1, all other things remaining the same, even the number of points in the grid. Wavenumbers are extracted the same way as before Figure 6b and then a x and a 0 are identified from this new solution field figures 8b and 8a.

It comes that the values found at low frequencies are far better than with a domain of side L=1. The limit at higher frequency comes from the spatial sample rate.

Figure 6c shows the dispersion relation when INCOME is used with a grid of 60 by 60 points, so a far lower sample rate. For frequencies higher than 30 Hz the wavenumbers calculated don't match with the analytic values so the identification process won't work.

Orthotropic membrane

For the second example, a x = 1, a y = 0.5, a xy = 0 and a 0 = 1. So a y becomes half of a x . This way the problem is not symmetric or isotropic. The problem can be interpreted as the behavior of an orthotropic membrane. The wavenumbers are calculated by INCOME with the same kernel [START_REF] Heckbert | Fourier transforms and the fast fourier transform (fft) algorithm[END_REF]. In contrast to the previous case, the values are not the sames depending on the angle. In figure 9a, in each direction, the calculated wavenumbers are very close to the analytic ones, except at low frequencies. Then again, the identification parameter method is then used to regain one parameter when the others are known. It is done for the parameter a 0 figure 10a and a y figure 10b. On this figure, it can be seen that as in the previous case, when the frequency is higher than 3, the results for the identification of a y and those for the identification of a 0 are near the expected values.

Case with cross derivative term

The last 2 nd -order equation tested adds the term a xy associated with a cross derivative term. Here a x = 1, a y = 1, a xy = 0.5 and a 0 = 1. The wavenumbers are calculated with the same kernel [START_REF] Heckbert | Fourier transforms and the fast fourier transform (fft) algorithm[END_REF] in different directions and presented figure 9b. For this case as well there is an angle dependency of the wavenumbers but the values at angles 0 and π 2 are the same. Except at low frequencies, the wavenumbers calculated are near the expected values even if the values calculated ay 0 and π 2 are not exactly the same. The parameters a 0 and a xy are identified using the method. The values identified for a 0 figure 11a fit very well but there is a bigger error for the identification of a xy even at higher frequencies, figure 11b.

4 th -order equations

Now, equations of the 4 th -order are studied. They have the following form.

D x ∂ 4 u ∂x 4 + D y ∂ 4 u ∂y 4 + D xy ∂ 4 u ∂x 2 ∂y 2 + a x ∂ 2 u ∂x 2 + a y ∂ 2 u ∂y 2 -ω 2 a 0 u = 0 (25)
The dispersion relation is solution of the the polynomial equation [START_REF] Boukadia | An inverse convolution method for wavenumber extraction (income): Formulations and applications[END_REF].

k 4 θ D x cos 4 θ + D y sin 4 θ + D xy cos 2 θ sin 2 θ -k 2 θ a x cos 2 θ + a y sin 2 θ = ω 2 a 0 (26)
Like for the 2 nd -order cases, the equation has been numerically solved on a unitary square domain, with a unitary load at the center. The boundary conditions are simply supported, a condition on the solution and its second derivatives. In order to get wavenumbers high enough compared to the dimensions of the domain, the frequencies go to 20 kHz. Two configurations are studied.

D∆ 2 u = ω 2 a 0 u (27)
And the dispersion relation is a direct relation.

k 4 = ω 2 a 0 D (28)
With, in this particular case a 0 and D having unitary values. The solution field of this equation at 500Hz has been presented figure 5a. An other convolution kernel has to be used to apply INCOME for this type of equation in order to get the wavenumbers.

Equation ( 27) is the equation of an isotropic plate, that is why there is no variation of wavenumbers with the direction in figure 12a. The calculated wavenumbers are very near the analytic values. For equation of 4 thorder, the wavenumber isn't a linear function of the frequency but a square root function. The identification results for a 0 and D x is presented figures 13a and 13b. They show a good match with the expected values.

Orthotropic plate

The next set of values is D x = 1, D y = 0.5, D xy = 2, a x = 0, a y = 1 and a 0 = 1. It adds an orthotropy and a 2 nd -order term. This way the problem is not symmetric and the dispersion relation needs a polynomial resolution. The wavenumbers are calculated with the same kernel [START_REF] Multiphysics | Introduction to comsol multiphysics®[END_REF]. Even if the wavenumbers are a little more distant from the expected values on figure 12b than in the previous case figure 12a, the identification results of the terms a 0 and D y presented figures 14a and 14b are still very near the exact values. This shows that the identification process seems to work well, limited in frequency by the size of the domain and the spatial sample rate. The last illustration is a 4 th -order equation from [START_REF] Boutin | Dynamics of ribbed plates with inner resonance: Analytical homogenized models and experimental validation[END_REF] that models the flexural behavior of a bi-ridded plate illustrated figure 15. In this equation there is an unusual frequency dependant term. A bi-directional ribbed plate is a periodic structure made of a plate on which two perpendicular rows of 245 stiffeners have been added. The beam grid of the stiffeners which is supposed to be stiffer than the little plate put between them, will dominate the dynamics of the structure. This presented model of the flexural behavior of this type of structure is a model of an equivalent orthotropic plate obtains from an asymptotic homogenisation process. Because the global dynamics is mostly linked to the stiffener grid, this orthotropic plate model come from the one of an orthotropic plate equivalent to a beam grid.

All the stiffeners in the same direction have the same geometrical and mechanical properties. In the x or y direction, h j is the height of the stiffeners, b j the width and l j the contact length with an internal plate. The mechanical parameters of the stiffeners E j , ν j and ρ j are respectively the Young Modulus, the Poisson ratio and the volumic mass. The other mechanical parameter can be calculated from the previous ones: the bending inertia I j = b j h 3 j 12 , the torsional modulus G j , the torsional iner-

tia I j = b j 2 h j 2 3 16 3 -3.36 h j b j 1 - h 4 j 12b 4
j , the linear mass Λ j = ρ j b j h j and the torsional constant

ρ j J j = ρ j h j b j b 2 j + h 2 j 12
. The formulation of the torsional inertia of a beam of rectangular section comes from [START_REF] Young | Roark's formulas for stress and strain[END_REF]. With these parameters the equation of the model of an orthotropic plate equivalent to a beam grid is the following one [START_REF] Young | Roark's formulas for stress and strain[END_REF].

E x I x l y ∂ 4 u ∂x 4 + E y I y l x ∂ 4 u ∂y 4 + G x I x l y + G y I y l x ∂ 4 u ∂x 2 ∂y 2 +ω 2 ρ y J y l x ∂ 2 u ∂x 2 + ρ x J x l y ∂ 2 u ∂y 2 = ω 2 Λ x l y + Λ y l x u (30)
The homogenised model of an orthotropic plate equivalent to a bi-directional ribbed plate is obtained by enriching this model [START_REF] Young | Roark's formulas for stress and strain[END_REF] with the internal dynamics from the internal plates. They are supposed all the same, of side length l x , l y and width d. Their mechanical properties are the followings: The Young modulus E p , the Poisson ratio ν p , the volume mass ρ p , the surface mass Λ ′ p = ρ p d, the bending inertia I p = d 3 12 and the modified module

E ′ p = E p 1 -ν 2 p
. The 4 th -order partial differential equation is then (31).

E x I x l y ∂ 4 u ∂x 4 + E y I y l x ∂ 4 u ∂y 4 + G x I x l y + G y I y l x ∂ 4 u ∂x 2 ∂y 2 +ω 2 ρ y J y l x ∂ 2 u ∂x 2 + ρ x J x l y ∂ 2 u ∂y 2 + K bc ωx ∂u ∂x + K bc ωy ∂u ∂y (31) = ω 2 Λ x l y + Λ y l x + Λ ′ p ϕ bc ω u
This equation contains the flexural rigidity terms of the stiffeners 4 and and an added mass term Λ ′ p ϕ bc ω that gives the effect of the dynamics of the internal plates. The 1 st -order term appears when the internal plates are non symmetric.

E x I x l y ∂ 4 u ∂x
E y I y l x , a term of torsional rigidity G x I x l y + G y I y l x ∂ 4 u ∂x 2 ∂y 2 , a rotational inertia term ω 2 ρ y J y l x ∂ 2 u ∂x 2 + ρ x J x l y ∂ 2 u
Here the boundary conditions are supposed square (l x = l y ) and clamped in every direction so the term becomes null K bc ω = 0. ϕ bc ω is calculated as the mean of the out of plane displacement on an internal plate with a unitary displacement of the beam grid. This term is function of the frequency and depends on the boundary conditions of the internal plates. With plates clamped at all their sides there is no analytic formula so it is approximated by the case of a circular plate with a unitary clamp motion which gives: Moreover a damping of 5% is added to the Young modulus as an imaginary part. From the values of Tables 1 and2, the value of the coefficients of equation ( 33) can be determined.

D x

D y D xy a x a y 575 + 28.75i J 575 + 28.75i J 747.5 + 37.375i J 4.5 10 -5 kg 4.5 10 -5 kg The equation is solved on a square domain of side 1m, with simply supported boundary conditions and a unitary load at the center. The frequencies go to 10000Hz. The INCOME method is used on a grid of 600 by 600 points. The parameters can be identified using the wavenumbers. For example D x or a 0 (ω). On the dispersion relation it can be noticed that even if the trend seems to follow a square root function, the term a 0 (ω) brought some resonances. The internal resonance frequencies can be interpreted as the frequencies of the modes with non zeros mean of an internal plate. They can be calculated as the zeros of the denominator of (32). The first ones are at 275 Hz, 1070 Hz, 2398 Hz, 4256 Hz, 6648 Hz and 9570 Hz, it agrees with the dispersion curve. The frequency dependent term a 0 and another independent D x are identified and compared to their analytic values. For the identification, on figures 17a and 17b the values calculated are close to the expected values. However near the internal resonance frequencies where

( a )

 a Schema of the beam (b) Displacement for 3 frequencies : Analytic at 10 Hz (-), 50 Hz (-), 100 Hz (-), numerical resolution at 10 Hz (•), 50 Hz (×) and 100 Hz (+) (c) Frequency response function of the displacement at the center: analytic solution (-), numerical values (•)

Figure 1 :

 1 Figure 1: Comparison between analytic solution and numerical resolution: (b) displacement function of the position at 3 frequencies, (c) frequency response function

Figure 2 :

 2 Figure 2: Schema of the membrane

( a )

 a Numerical solution field at 5Hz (b) Analytic solution field at 5Hz from formula (15) (c) Sectional view of the solution at frequency 5Hz at the red line of figure 3b from numerical resolution (•) and analytic formula (15) (-) (d) Frequency response function of the displacement at the point (L/4, L/4): : analytic solution (-), numerical values (•)

Figure 3 :Figure 4 :

 34 Figure 3: Comparison between analytic and numerical solutions: solution fields at 5 Hz, (a) numeric, (b) analytic, (c) displacement function of the position, (d) frequency response function

( a )

 a Numerical solution field at 500Hz (b) Analytic solution field at 500Hz from formula (18) (c) Sectional view of the solution at frequency 500Hz at the red line of figure 5b from numerical resolution (•) and analytic formula (18) (-) (d) Frequency response function of the displacement at the point (L/4, L/4): : analytic solution (-), numerical values (•)

Figure 5 :

 5 Figure 5: Comparison between analytic and numerical solutions: solution fields at 500 Hz, (a) numeric, (b) analytic, (c) displacement function of the position, (d) frequency response function

Figure 6 :

 6 Figure 6: Dispersion relation, analytic and from INCOME for the 2 nd -order equation with ax = 1, ay = 1, axy = 0 and a0 = 1 : analytic dispersion relation (-), INCOME wavenumbers at angles 0 rad (•), π 4 rad (×) and π 2 rad (+).

  (a): L=1m on a grid of 400 by 400 points, (b): L=5m on a grid of 400 by 400 points and (c): L=1m on a grid of 60 by 60 points.

Figure 7 :Figure 8 :

 78 Figure 7: Identified values of a0 (a) and ax (b) calculated using the identification method on the solution fields at different frequencies, exact values (-) and identified values (•)

Figure 9 :Figure 10 :Figure 11 :

 91011 Figure 9: Dispersion relation, analytic and from INCOME for 2 nd -order equations in the directions 0, π 4 and π 2 : analytic dispersion relation at 0 rad (-), π 4 rad (-), π 2 rad (-), INCOME wavenumbers at 0 rad (•), π 4 rad (×) and π 2 rad (+). (a): ax = 1, ay = 0.5, axy = 0, a0 = 1 and (b): ax = 1, ay = 1, axy = 0.5, a0 = 1

3. 4 . 1 .

 41 Isotropic plateThe first equation is a case of the Kirchhoff Love plate equation that models the flexion in an isotropic plate.D x = 1, D y = 1, D xy = 2, a x = 0, a y = 0 ans a 0 = 1. In this case, D x = D y = D xy 2 = D.The solution fields are presented below. The equation can be simplified to make a bilaplacian appear.

Figure 12 :

 12 Figure 12: Dispersion relation, analytic and from INCOME for the 4 th -order equations in the directions 0, π 4 and π 2 : analytic dispersion relation at 0 rad (-), π 4 rad (-), π 2 rad (-), INCOME wavenumbers at 0 rad (•), π 4 rad (×) and π 2 rad (+). (a): Dx = 1, Dy = 1, Dxy = 2, ax = 0, ay = 0, a0 = 1 and (b): Dx = 1, Dy = 0.5, Dxy = 2, ax = 0, ay = 1 and a0 = 1

Figure 13 :Figure 14 :

 1314 Figure 13: Identified values of a0 (a) and Dx (b) calculated using the identification method on the solution fields at different frequencies, exact values (-) and identified values (•)

Figure 15 :

 15 Figure 15: Example of bi-directional ribbed plate of 5 by 5 cells

1

 1 (δr)J 1 (δr) I 1 (δr)J 0 (δr) + I 0 (δr)J 1 (δr) (32) with• J 0 , J 1 Bessel functions of the first kind• I 0 , I 1 modified Bessel functions of the first kind• δ 4 = ω 2 Λ p E ′ p I p = (2πf ) 2 Λ p (1 -ν 2 p ) E p I p • r = 0.532l 255The equation (31) can still be written in the following form.D x ∂ 4 u ∂x 4 + D y ∂ 4 u ∂y 4 + D xy ∂ 4 u ∂x 2 ∂y 2 + ω 2 a x ∂ 2 u ∂x 2 + ω 2 a y ∂ 2 u ∂y 2 -ω 2 a 0 (ω)u = 0 (33)But the term a 0 is frequency dependent and a x , a y are replaced by ω 2 a x and ω 2 a y . a 0 (ω) = Λ ′ + Λ ′ p ϕ CCCC ω (δ(ω)r). The values of the aluminium are chosen for the stiffeners and those of the Plexiglas for the internal plates. The stiffeners are supposed to be the same in the x and y directions.Parameter : h x , h y b x , b y E x , E y ν x , ν y , ρ x , ρ y

Figure 17 :

 17 Figure 17: Identified values of a0 (a) and Dx (b) calculated using the identification method on the solution fields for frequencies below 10 kHz, exact values (-) and identified values (•)

Figure 16 :

 16 Figure 16: Dispersion relation, analytic and from INCOME of the equation (33) in the directions 0, π 4 and π 2 : analytic dispersion relation at 0 rad (-), π 4 rad (-), π 2 rad (-), INCOME wavenumbers at 0 rad (•), π 4 rad (×) and π 2 rad (+)

Table 1 :

 1 Geometrical and mechanical properties values of the beam grid

	Parameter : l x , l y	d	E p	ν p	ρ p
	Value :	0.1m 0.001m 3GPa 0.3 1200 kg/m 3

Table 2 :

 2 Geometrical and mechanical properties values of the internal plates

Table 3 :

 3 Value of the parameters of equation (33) /m 2 1.2 kg/m 2 4.3571 -0.21786i s 2 /m 4 0.0532 m

	Λ ′	Λ ′ p	Λ ′ p (1 -ν 2 p ) E p I p	r
	5.4 kg

Table 4 :

 4 Value of the parameters in the calculation of a 0

ea is a mass coefficient and da a damping coefficient. They won't be used here as they deal with temporal derivatives

the identification of D x are less precise. The identification method works also for this more original model.

Acknowledgments

This work was supported by the LabEx CeLyA (Centre Lyonnais d' Acoustique, ANR-10-LABX-0060) of Université de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR).

Conclusion

In this paper, a methodology to implement any linear PDE using the COMSOL solver has been presented. In particular, it is possible to implement higher than 2 nd -order equation and frequency dependent coefficients. The validity of this method has been verified using test cases with known analytic solutions. Then a coefficient identification method for PDE based on an original approach has been described. This method is adapted for dimension 2 PDE in structural mechanics, like plates and membranes and provides a frequency dependent identification. The approach is based on the comparison between analytic and reference k-spaces which are extracted from displacement fields using the methods IWC and INCOME.

The identification methodology and in particular the k-space extraction method and analytic dispersion relation definition have been described and tested. Three cases of 2 nd -order equations, deflection of isotropic and orthotropic membranes, two of 4 th -order equations, flexural behavior in isotropic and orthotropic plates and an original equation that model the behavior of a bi-directional ribbed panel with some original frequency dependent terms have been used to tested the identification process. Those test cases have been numerically calculated using this implementation methodology of PDE resolution for finite element software. In all these cases, the coefficients are correctly retrieved. The main limitation of the method seems to be produced by the wavenumber extraction procedure. It is limited in low frequency by the domain size and in high frequency by the spatial sample rate.

This method has been developed in order to extract information on structures using displacement fields from experimental data. It can be used to characterise an unknown material plate. But in particular, it can be used in the modeling of meta-material structures where the PDE coefficients are unknown and can't be easily estimated. This way equivalent parameters can be calculated and then drastically reduce computation time of simulations of the structure.

In future works, computational time could be still reduced by improving the extremum research, both in the wavenumber calculation methods and in the identification. Moreover the process can implement the identification of several parameters at the same time using multi objective optimisation.

CRediT author statement

Appendix A. Short overview of k-space methods

In this appendix, k-space methods for 2 dimensional domains are presented with more details. In the case of 2-dimensional structures like plate or panel, the wavenumbers depend on the propagation direction. Some methods exist to get the k-space from a displacement field. The field is supposed known on a regular grid. The Fourier transform, the Inhomogeneous Wave Correlation method (IWC) ans INverse COnvolution MEthod (INCOME) are described.

Appendix A.1. Fourier transform

The Fourier is an integral transform that transforms a function of the time to a function, of the frequency or a function of the space to a function of the wavenumbers [START_REF] Heckbert | Fourier transforms and the fast fourier transform (fft) algorithm[END_REF]. The Fourier transform can be used to get the wavenumber However, it has some drawbacks. The imaginary part of the wavenumbers is neglected. It is not very adapted to finite domains. And the spatial sample step has to be below enough the studied wavelength (Nyquist-Shannon Theorem).

Appendix A.2. Inhomogeneous wave correlation method

The Inhomogeneous wave correlation method, IWC, [START_REF] Ichchou | Identification of effective sandwich structural properties via an inverse wave approach[END_REF][START_REF] Berthaut | K-space identification of apparent structural behaviour[END_REF] is based on the maximisation of a correlation function using an inhomogeneous wave ô. In two dimension, this wave is parameterized by the angle θ and the wavenumber at this angle k(θ)(1 + iγ(θ)).

ôk,γ,θ (x, y) = e -ik(θ)(1+iγ(θ))(x cos θ+y sin θ) (A.1)

Naming u the displacement field and û its Fourier transform, the correlation function to maximize at each angle needed is the following (A.2).

In practice, at a given frequency, the displacement field û is only known on a finite number of points so the integrals are calculated by numerical methods like the trapezoidal method. Moreover the term S |û| 2 dx dy is not calculated because it is independent of θ, k and γ. The maximum is obtained by evaluation. At the end the method gives for each θ a value k(θ) and a value for γ(θ). Then the complex wavenumber is k(1 + jγ).

Appendix A.3. INverse COnvolution MEthod

The INverse COnvolution MEthod, INCOME, [START_REF] Boukadia | An inverse convolution method for wavenumber extraction (income): Formulations and applications[END_REF][START_REF] Boukadia | Vibroacoustics of periodic structures: model order reduction, characterization, optimization[END_REF] is an other method to get the k-space from a displacement field. u is the discrete displacement field at a given frequency known on a grid. The aim of the method is to find a convolution kernel S such that :

Or at least minimize ∥u * S∥ 2 . Then the k-space is deduced from S. The form of S is determined by the dimension of the field and the hypothesis on the equation. For example the following kernel can be used for a field solution of the equation of the flexural behavior of an orthotropic plate. After finding the coefficient of the convolution kernel, the wavenumber in a direction θ will come from the solution of the minimization of the module of a polynomial equation. For the example of an orthotropic plate, with a grid of mesh size L x L y , defining λ x (k, θ) = e ik cos θLx and λ y (k, θ) = e ik sin θLy the equation to minimize is:

INCOME is less time consuming than the other methods so the number of frequencies studied and the number of angle of measurements can be increased. However the form of its convolution kernel must be chosen function of the form of the model equation.