O Dorival

P Rouch

O Allix

A substructured Trefftz method for updating joint models in the medium-frequency range

The response of a structure to medium-frequency vibrations is highly dependent on several phenomena. One of these, the behavior of joints between substructures, is essential because it controls the distribution of the power injected into the structure and governs most of the dissipative effects. This paper introduces a new method for the identification of joint parameters from experimental vibration data, inspired by previous works on finite element model updating using the error in the constitutive relation. Since such numerical techniques are not suitable for the medium-frequency range, in which the finite element mesh must be refined, our work uses the numerical framework of the variational theory of complex rays, which is a Trefftz approach entirely dedicated to the calculation of medium-frequency vibrations at very low cost. The main scope of this paper is the presentation of the formulation and its validation against actual and numerically-simulated experimental results.

Introduction

Predicting the behavior of complex structures is very important for both engineers and scientists. Among the various types of problems, the response to medium-frequency vibration loading poses several specific difficulties. Such problems are frequently encountered, for instance, in the aviation industry and in aerospace engineering.

The first difficulty in the medium-frequency range is due to an insufficient knowledge of some physical phenomena the structure depends on, or to the impossibility of taking their level of complexity into account in reasonable calculations. For example, dissipative phenomena play a very important role in the response of a vibrating structure, not only in the medium-frequency range but also at low and high frequencies. These are due, in part, to internal material dissipation governed by very small-scale phenomena, which do not lend themselves to the use of available microscopic models in structural calculations. Classically, material damping is simulated using viscous or hysteretic damping. But a large part of the dissipation, which can represent up to eighty or ninety percent of the total energy loss, takes place at interfaces between substructures [START_REF] Esteban | Theoretical modeling of wave propagation and energy dissipation in joints[END_REF][START_REF] Nanda | Study of damping in layered and jointed structures with uniform pressure distribution at the interface[END_REF], such as bolted, (spot-)welded or riveted joints. Most of the local phenomena involved are known only very approximately. Moreover, they are nonlinear and highly dependent on the geometric and material properties as well as on the loading conditions at the interface. A classical way of obtaining a reliable response consists in using a simple model whose parameters are derived from experimental data through an identification method. The objective of this paper is to introduce an updating method dedicated to the improvement of joint models used in medium-frequency vibrations. The vast number of possible strategies for doing this falls into two categories: direct methods (also called identification methods) and indirect (or parametric) methods. Depending on the category, different quantities can be used to built a cost function for the comparison of experimental and numerical solutions: input residuals [START_REF] Berger | Méthodes de recalage de structures en dynamique[END_REF][START_REF] Farhat | Updating finite element dynamic models using an element-by-element sensitivity methodology[END_REF], output residuals [START_REF] Piranda | Parametric correction of finite element models by minimization of an output residual: Improvement of the sensitivity method[END_REF][START_REF] Lammens | Application of a FRF based model updating technique for the validation of a finite element model of components of the automotive industry[END_REF] or energy residuals [START_REF] Pascual | A frequency domain correlation technique for model correlation[END_REF][START_REF] Balmès | Review and evaluation of shape expansion methods[END_REF]. Our approach is an updating method based on the error in the constitutive relation [START_REF] Ladevèze | Recalage de modélisations des structures complexes[END_REF], which was originally used to quantify finite element approximation errors.

A second type of difficulty encountered in medium-frequency calculations is due to the small-wavelength phenomena involved in the response of the structure in this frequency range. In the low-frequency range, where the wavelengths are more-or-less comparable to the size of the structure, the local description of the motion requires only a small number of degrees of freedom. Finite element methods, often in conjunction with modal synthesis, are the most common techniques used for the resolution of low-frequency problems, but in order to avoid approximation errors known as the pollution effect [START_REF] Ihlenburg | Dispersion analysis and error estimation of Galerkin finite element methods for Helmholtz equation[END_REF] the cost of these techniques increases dramatically with the frequency. When the frequencies are very high, the description of the structure can be averaged within a frequency range and/or on a spatial patch, and statistical strategies, such as Statistical Energy Analysis (SEA) [START_REF] Lyon | Theory and Application of Statistical Energy Analysis[END_REF], can be used to achieve good simulations at very low cost: typically, one unknown is used per subsystem (its energy). In the medium-frequency range, however, the basic assumptions of the SEA must be relaxed because the structure still has modal behavior. Therefore, medium-frequencies constitute an interesting field for many numerical researchers. Some of the proposed approaches (e.g. the discontinuous enrichment method [START_REF] Farhat | A discontinuous galerkin method with Lagrange multipliers for the solution of Helmholtz problems in the midfrequency regime[END_REF], the generalized finite element method [START_REF] Strouboulis | The generalized finite element method for Helmholtz equation: Theory, computations and open problems[END_REF], Trefftz methods [START_REF] Kita | Trefftz methods: An overview[END_REF] or boundary element methods [START_REF] Bonnet | Boundary Integral Equation Methods in Solids and Fluids[END_REF]) use techniques based on the finite element framework. In the SEA category, one should mention the works on wave intensity analysis [START_REF] Langley | A wave intensity technique for the analysis of high frequency vibrations[END_REF] and on power flow analysis [START_REF] Nefske | Power flow finite element analysis of dynamic systems: Basic theory and application to beams[END_REF]. Our approach uses the Variational Theory of Complex Rays (VTCR) [START_REF] Ladevèze | The variational theory of complex rays for the calculation of medium-frequency vibrations[END_REF], a Trefftz method which enables the calculation of medium-frequency vibrations at relatively low and frequency-independent cost.

This paper deals with a new method for updating joint models in the medium-frequency range. Our identification method is based on the error in the constitutive relation, which was initially developed for the validation of finite element calculations. The transposition of this concept to the improvement of finite element models based on experimental data [START_REF] Ladevèze | Error in the constitutive relation in dynamics[END_REF] produced robust results compared to other parametric updating techniques such as input or output error methods. In the context of medium-frequency vibrations, our objective was, for the reasons described above, to avoid using finite element techniques. Therefore, the original formulation was adapted to the VTCR framework. The principle consists in a domain decomposition approach in which the joint to be identified is separated from the reliable substructures. The reliable information is satisfied exactly, whereas a weak solution of the unreliable information is sought through the minimization of a residual. This residual is used as a cost function in an iterative correction process, leading to an inverse problem.

The paper is structured as follows: in the second section, we introduce the notations and present the available numerical and experimental information sets. The equations are split into two groups depending on their reliable or unreliable status. In the third section, we build an admissibility space in order to satisfy the reliable equations exactly. Since our numerical approach for the calculation of medium-frequency vibrations is not standard, we also review the numerical framework of the VTCR. In the fourth section, we present the details of the inverse problem of the minimization of a cost function constructed from the unreliable equations. Finally, in the last section, we present examples of joint identification in which the experimental data were simulated numerically (with the exception of the last example, which uses actual experimental data). We conclude that the identification process is robust and that the resolution is inexpensive. These encouraging results should now be confirmed on larger structures using actual experimental data, which is what we are currently working on.

Review of the information available for the inverse problem

For the sake of simplicity, let us express the problem in the classical threedimensional framework. Let Ω be a structure composed of two substructures Ω 1 and Ω 2 connected by a joint denoted Γ (see Figure 5). The available information is of two types. On the one hand, through the modeling process of the actual structure, we have the equations of the continuous model, which can be solved either theoretically (if possible) or numerically. On the other hand, tests provide an experimental model, which is imperfect in the sense that it is incomplete and contains measurement noise. Usually, a linear elastic vibration problem at a fixed angular frequency ω is described using complex quantities whose harmonic parts with respect to time reduce to zero in the subsequent equations.

Γ F 2d u 1d Ω 1 Ω 2 ∂ F Ω 2 ∂ u Ω 1 ∂ F Ω 1 ∂ u Ω 2 u 2d F 1d
Fig. 1 The reference problem for a structure Ω

The available information

Governing equations of the continuous model Let u i and σ i denote respectively the displacement field and the stress field of Substructure Ω i . The governing equations to be satisfied consist of substructure equations and interface equations. Displacements u id are prescribed on a part ∂ u Ω i of the substructure boundary ∂Ω i , while the complementary part

∂ F Ω i = ∂Ω i -(∂ u Ω i ∪ Γ ) is subjected to prescribed surface forces F id .
In the absence of body forces, the dynamic equilibrium equation and the constitutive relation between u i and σ i involve the mass density ρ i and the Hooke's tensor K i . Then, over the substructures Ω i , (u i , σ i ) must satisfy:

(u i , σ i ) ∈ [H 1 (Ω i )] 3 ⊗ [L 2 (Ω i )] 6
set of finite energy fields

div σ i = -ρ i ω 2 u i within Ω i σ i = (1 + j η i ) K i : ε(u i ) within Ω i u i = u id over ∂ u Ω i σ i • n i = F id over ∂ F Ω i ( 1
)
where ε is the strain operator classically defined as the symmetric part of the gradient of the displacement fields. Here, structural dissipation is introduced into the constitutive law under assumption of hysteretic damping, and η i designates the structural loss factor (also called structural damping in this paper) which can depend on the angular frequency ω. Note that other models of internal damping could have been used here.

In the case of a perfect interface, Γ is described through the equations of static equilibrium of the strain fields and continuity of the displacements. Here, in order to take more realistic situations into account, the behavior of Interface Γ will be modeled by its own dynamic equilibrium equation and constitutive law:

σ 1 • n 1 + σ 2 • n 2 = -ω 2 µ u m over Γ σ m = -k (u 2 -u 1 ) over Γ (2)
where u m and σ m represent the mean values of displacements and stresses over Γ :

u m = u 1 + u 2 2 over Γ σ m = σ 1 • n 1 -σ 2 • n 2 2 over Γ (3)
In the general three dimensional case µ and k represent respectively the mass and the stiffness per unit area of the interface and can take into account dissipative behavior by including an imaginary part. (However it should be noted that in the applications treated at the end of this paper µ and k are given per unit length since the substructures are two dimensional). Contrary to the substructure case, it is difficult to link one's knowledge of the interface to geometric and material parameters alone because of the complexity of the phenomena localized in the joint's behavior (friction, heating, abrasion. . . ) Therefore, the interface's parameters must be identified using experimental results.

Experimental information

In order to obtain an inverse problem, additional information is sought through an experimental study of the vibrational response of the structure, so that the problem becomes overspecified. Let a denote the experimental value of some quantity a. This complex representation takes into account both the magnitude and the phase. Whatever form the measurement takes, we will assume that the measured quantities are:

-the angular frequency ω, -the excitation forces F d , -the displacements u( X n ) at several locations X n within the structure.

Moreover, let us consider that these measurements are based on assumed locations X n and on assumed measurement orientations θ n of the measurement points. These assumptions may be satisfied only approximately or may be limited to a certain precision.

The measurement process can be viewed as a projection of the real displacement field u(X) onto points X n in the directions θ n . Let Π u ( e Xn, e θn)

denote the corresponding projector:

Π u ( e Xn, e θn) : u -→ u( X n ) • θ n (4)
Finally, one must also take into account the location X F and orientation θ F of the excitation. Let Π F ( e XF, e θF) denote the projector associated with the measurement of the excitation.

Partitioning between reliable and unreliable information

Our method for updating the interface model is based on the modified error method presented in [START_REF] Ladevèze | Recalage de modélisations des structures complexes[END_REF] in the context of identification. This approach consists in dividing the available information into two groups by separating the (relatively) reliable quantities from the unreliable quantities. The distinction between the two groups depends on the problem being treated: thus, the engineer's experience can be of great importance for this step.

Concerning the continuous model, considering one's knowledge of the material and the large domain of validity of simplified models (beams, plates, shells,. . . ), the substructures will be assumed to be modeled properly and their material and geometric parameters to be known precisely. This information will be assumed to be reliable. Conversely, the equations of an interface can be viewed as a global behavior model whose parameters are known only very imperfectly. Therefore, these equations will be considered to be unreliable information.

Similarly, the experimental data are split into two groups. In general, the angular frequency ω is measured very precisely. In this paper, we will assume that the main errors come from the unavoidable noise which affects the magnitude of the measurements. The values of F d and u( X n ) will be considered questionable, whereas the locations X n and orientations θ n of the measurement points as well as the location X F and orientation θ F of the exciting forces will be considered to be known with sufficient precision. Table 1 summarizes the categorization of the different types of information. 1 The reliable and unreliable information for the inverse problem Depending on the type of problem being treated, other choices can be made: in [START_REF] Chouaki | Model updating using the error in the constitutive relation: Real case study[END_REF], the dissipation error enabled the elastic part of the constitutive law to be considered reliable in order to identify the damping phenomena. In [START_REF] Nguyen | Application of the CRE for parameter identification in nonlinear dynamics with corrupted measurements[END_REF], highly corrupted measurement data were used to identify the viscoplastic behavior in dynamics.

S 1 Γ F 2d u 1d Ω 1 Ω 2 ∂ F Ω 2 ∂ u Ω 1 ∂ F Ω 1 ∂ u Ω 2 u 2d F 1d w 2 w 1 S 2
Fig. 2 The substructured problem

The key to the method is to seek a compromise solution for the above equations in the sense that this solution satisfies the reliable equations exactly, but satisfies the other equations of the continuous and experimental models only approximately. Then, this solution enables us to identify which among all the unreliable equations is the most inexactly satisfied and to adjust its parameters, which leads to the inverse problem. The numerical strategy we chose in order to solve this problem at a low and predictable cost is a Trefftz approach, called the variational theory of complex rays, in which all the boundary equations are combined into a single variational form. In order to treat the substructures and their interfaces separately, we will now introduce a substructured version of this technique.

The interface quantities

The principle, which was first introduced in [START_REF] Dorival | A substructured version of the variational theory of complex rays dedicated to the calculation of assemblies with dissipative joints in the medium-frequency range[END_REF] for the calculation of dissipative and heterogeneous joints, is based on a condensation approach in which the interface is treated like an actual substructure with its own equations, parameters, unknowns and discretization. Figure 2 illustrates the proposed substructured approach. Let w i and S i denote respectively the interface's displacement and stress unknowns relative to Substructure Ω i . (One unknown is defined for each side of Interface Γ .)

The continuous problem is composed of two parts. First, the interface equations ( 2) and (3) (i.e. the unreliable equations) yield:

S 1 + S 2 = -ω 2 µ w 1 + w 2 2 over Γ S 1 -S 2 2 = -k (w 2 -w 1 ) over Γ (5)
Then, the substructure equations (i.e. the reliable equations) are completed in order to ensure compatibility between the substructure unknowns (u i , σ i ) and the interface unknowns (w i , S i ). Among several possibilities, we chose a condensation approach in which the interface quantities w i and S i are prescribed on the part ∂ Γ Ω i of the boundary:

div σ i = -ρ i ω 2 u i within Ω i σ i = (1 + j η i ) K i : ε(u i ) within Ω i u i = u id over ∂ u Ω i σ i • n i = F id over ∂ F Ω i u i = w i over ∂ Γ Ω i σ i • n i = S i over ∂ Γ Ω i (6)
With this substructuring technique, we can deal with the reliable substructures in any desired numerical framework without taking the interface equations into account. Now, we will go through the schemes we chose to solve each group of equations.

3 Satisfaction of the reliable information

The reliable equations of the model

The numerical approach used to solve Equations ( 6) is the VTCR [START_REF] Ladevèze | A new computational approach for structure vibrations in the medium frequency range[END_REF]. This method was developed as a very inexpensive approach to medium-frequency vibration problems, which would otherwise be very cumbersome due to the medium-frequency framework and to the iterative procedures classically encountered in identification problems. A key feature of the VTCR which turns out to be very convenient for our strategy is that it allows the local equations of the substructure to be satisfied exactly. This enables us to build shape functions which have a strong physical meaning. The boundary conditions on the substructures are introduced through a variational approach and satisfied by discretized shape functions, resulting in a finite-dimension problem.

Construction of admissible fields

The construction of VTCR shape functions is detailed in [START_REF] Ladevèze | Analysis of mediumfrequency vibrations in a frequency range[END_REF]. Here, we will limit ourselves to a brief description. The displacement/stress pair v i = (u i , σ i ) is chosen in the spaces of the finite-energy fields U i × S i in order to satisfy both the local equilibrium equation and the substructure's constitutive law exactly.

S i ad = (u i , σ i ) ∈ U i × S i such that div σ i = -ρ i ω 2 u i σ i = (1 + j η i ) K i : ε(u i ) (7)
Such solutions are sought as two-scale approximation solutions:

u i (X) = U i (X, k) e k(ω)•X within Ω i ( 8
)
Fig. 3 Examples of VTCR complex rays for a rectangular plate where U i is the large scale of the quantities which is assumed to vary slowly, and the exponential part e k(ω)•X corresponds to the small scale, i.e. the oscillating part of the vibration. Note that the wave number k is introduced in a complex form: the imaginary part of k gives an oscillating term and the real part of k -which is negative -stands for the evanescent part of the wave together with the extenuation due to the structural damping. The reason for not putting this non-oscillating part in U i is that it will be analytically known, whereas U i will remain an unknown in the problem. The use of this approximation in Equations ( 7) leads to the domain of definition of k(ω). For instance, the application to the vibrations of isotropic plates in bending yields:

(1 + j η i ) E i h 2 i 12 (1 -ν 2 i ) ∆∆u i = ρ i ω 2 u i (9)
where E i , ν i and h i denote the Young's modulus, Poisson's ratio and thickness of the plate respectively. Then, the complex wave number k must satisfy the classical dispersion equation:

|k| 4 = 12 (1 -ν 2 i ) ρ i ω 2 (1 + j η i ) E i h 2 i (10)
For plates, these functions describe the vibration solutions in an unbounded plane, a half-plane or a quadrant (see Figure 3). The oscillating part of what we call complex rays is known explicitly. Consequently, U i (X, k) vary very slowly and can be described by only a few unknowns.

Variational form of the boundary conditions

The boundary conditions on ∂Ω i are satisfied weakly for any test function δv i through the following variational equation, which expresses the compatibility of Substructure Ω i with the prescribed interface forces S i :

Compatibility of Substructure Ω i with the interface Find v i ∈ S i ad such that:

A i (δv i , v i ) + B i (δv i , S i ) = L i (δv i ) ∀ δv i ∈ S i ad ( 11
)
where A i and B i are nonsymmetrical bilinear operators defined by:

A i (δv i , v i ) = Re -jω ∂uΩi (δσ i • n i ) • u * i dl + ∂FΩi (σ i • n i ) • δu * i dl + ∂Γ Ωi (σ i • n i ) • δu * i dl (12)
B i (δv i , S i ) = Re jω ∂Γ Ωi S i • δu * i dl (13)
and L i is a linear form given by:

L i (δv i ) = Re -jω ∂uΩi (δσ i • n i ) • u * id dl + ∂FΩi F id • δu * i dl (14)
Re(a) and a * designate respectively the real part and the conjugate of a complex quantity a. This variational formulation can be viewed as the expression of the balance of the virtual powers at the boundary of Substructure Ω i . Let us note that the use of the real part Re(-jω) is not an approximation because the imaginary part would lead to exactly the same system. The variational form [START_REF] Kita | Trefftz methods: An overview[END_REF] is equivalent to the substructure problem (6) under the following conditions which, in practice, are not restrictive [START_REF] Ladevèze | The variational theory of complex rays for the calculation of medium-frequency vibrations[END_REF]:

-the substructure reference problem (6) has a solution, -Hooke's tensor K i is positive definite, -the structural damping ratio η i is strictly positive.

The reliable experimental data

Now, we will define a cost function in order to quantify the satisfaction of the unreliable equations. We will compare the numerical solution and the experimental solution using the projector Π u ( e Xn, e θn) defined in 2.1. In order to ensure that the measurement points and their orientations are reliable, we define the corresponding projectors Π u for the numerical model and choose the same points X n and orientations θ n :

Π u ( e Xn, e θn) : u -→ u( X n ) • θ n (15)
The reliable location X F and orientation θ F of excitation force F are applied by prescribing the corresponding force in the numerical simulation.

Finally, the measured angular frequency ω is enforced precisely by setting exactly the same value in the VTCR simulation. Therefore, both the variational formulation [START_REF] Kita | Trefftz methods: An overview[END_REF] and the construction of admissible complex rays (7) use the experimental value ω = ω of the frequency.

Construction of the cost function from the unreliable information

The unreliable information on the models will be satisfied only approximately through a compromise solution obtained by minimizing a cost function composed of two terms: an interface residual (which quantifies the satisfaction of the interface equations) and a measurement term (which represents a distance between the experimental and numerical solutions).

The interface residual

The interface residual E 2 int , also called the model error, is obtained as an L 2 error in the interface equations ( 5):

E 2 int (w, w) = 1 N 1 Γ S 1 + S 2 + ω 2 µ w 1 + w 2 2 2 dl + 1 N 2 Γ S 1 -S 2 2 + k (w 2 -w 1 ) 2 dl (16)
where N 1 and N 2 are normalisation/weighting parameters used to balance the two terms of the interface error. Let us note that this residual is a quadratic form of the unknowns w i and S i . For the purpose of this paper, N 1 and N 2 were chosen equal to 1.

If the interface error is zero, the solution satisfies the interface equations exactly. From a technical point of view, this residual is polluted by approximation errors (essentially discretization errors) due to the numerical resolution of the problem. Then, the magnitude of the residual must be compared to the discretization error of the VTCR.

The measurement term

The measurement term is constructed as a distance between the numerical and experimental results obtained by Projectors Π and Π respectively:

E 2 mes (v, v) = Π u (u) -Π u ( u) 2 u + Π F (u) -Π F ( u) 2 F = n u( X n ) • θ n -u( X n ) • θ n 2 u + F d -F d 2 F ( 17
)
The norms are chosen such that the two terms are balanced. Finally, the measurement term is adjusted in order to be expressed in the same units as the interface term.

If this error is zero, the numerical solution matches the experimental solution exactly. In practice, this residual is inevitably polluted by measurement noise.

The final direct and inverse problems

The solution of the modified problem is obtained by minimizing the cost function, which is the weighted sum of the two residuals defined previously:

E 2 mod = (1 -r) E 2 int + r E 2 mes ( 18
)
E 2 mod is also called the modified residual because in the absence of the measurement term it would be equivalent to the continuous reference problem. Parameter r can be used to assign greater confidence to the interface error or to the measurement error depending on the precision of the experimental setup and one's knowledge of the interface's behavior.

The final problem consists in minimizing this residual over the unreliable information while satisfying all the reliable equations:

The modified problem Find (u i , σ i , w i , S i ) such that:

• (u i , σ i , w i , S i ) satisfy the admissible equations [START_REF] Kita | Trefftz methods: An overview[END_REF],

• (u i , σ i , w i , S i ) minimize the modified residual [START_REF] Langley | A wave intensity technique for the analysis of high frequency vibrations[END_REF].

(19)

If r = 0, the solution does not depend on the experimental data. If r = 1, the solution does not take into account the proposed interface model. In all other cases, the solution is a compromise among the sets of unreliable equations. An often used value for r is about 0.8. This will be our choice in the rest of the paper.

The calculated solution and the magnitudes of the terms of the residual depend on the assumed values of the interface's mass and stiffness parameters µ and k. The inverse problem consists in finding the optimum values of these interface parameters. Then, the calculated residual becomes a cost function in an optimization scheme:

The inverse problem Find the optimum interface parameters p opt = (µ opt , k opt ) such that:

p opt minimizes the cost function E 2 mod (p) ( 20
)
This step is nonlinear and can be solved by classical optimization schemes, such as Newton methods.

Discretization of the problem

The previous problem must be solved numerically. This section describes the discretization process which leads to a finite-dimension problem.

Discretization of the substructure unknowns (u, σ)

The discretization of the substructure unknowns (u, σ) is obtained by using a subspace of the admissibility space S ad defined in Equations [START_REF] Esteban | Theoretical modeling of wave propagation and energy dissipation in joints[END_REF]. Following Equation ( 8), these admissible solutions are decomposed according to a small scale e k(ω)•X and a large scale U(X, k), which becomes the unknown. First, we assume in this paper that the large scale U varies sufficiently slowly with respect to the spatial coordinates that it can be considered to remain constant with respect to X. Let us note that according to previous papers on the VTCR this assumption gives very accurate results, but other possible descriptions (e.g. a polynomial dependence with respect to X [START_REF] Ladevèze | The variational theory of complex rays for the calculation of medium-frequency vibrations[END_REF]) have also been studied.

Discretization involves making assumptions about the variation of the large scale U with respect to the propagation direction k. For example, in the case of interior complex rays, real part of k follows a curve C int which is a circle of radius k 0 (see Figure 4 left) given by:

k 4 0 = 3 ρ (1 -ν 2 ) E h 2 (21)
For edge and corner complex rays, the corresponding loci of k are denoted C edg and C cor . The continuous description of the displacement field u is given by:

u(X) = k(α)∈ Cint∪C edg ∪Ccor U(k(α)) • e k(α)•X dα (22)
Then, curve C int is split into N u non-overlapping angular sectors denoted

C i = [α i , α i+1 ].
The discretized description of the large scale U is chosen to be constant over each angular sector C i . For example, Figure 4 (right) shows an approximation in which C int is composed of eight equal angular sectors. Then, the displacement can be approximated by the following finite dimension approximation u h : where Φ i is a known shape function. This choice retains for k an integral representation in which no propagation direction is missing. As in the case of finite elements, the representation can be enriched by refining the mesh over C int , i.e. by increasing the number of angular sectors N u , which is very similar to the h-method. Another option would be to choose a more complex representation (e.g. linear or quadratic) for a particular sector, which can be viewed as a VTCR equivalent of the p-method. Let us note that the final representation of u also involves shape functions Φ i relative to edge and corner complex rays which are obtained in a very similar way. Previous studies showed that if such local modes are missing the medium-frequency response cannot be calculated correctly.

u h ∈ S h ad ⇔ u h (X) = Nu i=1 U i • α∈Ci e k(α)•X dα = Nu i=1 U i • Φ i (X)

Discretization of the interface unknowns (w, S)

In order to obtain a finite-dimension formulation, the interface quantities (w, S) must also be discretized. Since (w, S) are equivalent to the substructure quantities (u, σ), it seems natural to choose to represent them using the same type of approximation space. From a theoretical point of view, one could seek w h in the same approximation space S h ad as u h or in a subspace of S h ad . For example, one could omit some modes which are far removed from Interface Γ . According to our personal experience, such omissions could lead to inaccurate results in certain situations. In practice, they would not be worthwhile anyway since the number of complex rays remains quite reasonably small (less than 100 complex rays per substructure unknown, plus the same number for the associated interface quantity). Consequently, in this paper, w h will be sought in the same way as u h , but the interface shape functions will be denoted Ψ i to indicate that they can be different from the substructure shape functions Φ i :

w h ∈ S h ad ⇔ w h (X) = Nw i=1 W i • Ψ i (X) (24)

The associated matrix problems

Applying the previous discretizations ( 23) and ( 24) to the substructure/interface compatibility equations [START_REF] Kita | Trefftz methods: An overview[END_REF], one gets the following linear system which expresses substructure admissibility:

K A k, l U l + K B k, l W l = F k ( 25
)
where the matrix terms come from the operators involved in the variational formulation [START_REF] Kita | Trefftz methods: An overview[END_REF]:

K A k, l = A(Φ k , Φ l ) K B k, l = B(Φ k , Ψ l ) F k = L(Φ k ) (26)
By defining the complete vector of the unknowns Y T = U T W T and by assembling the complete admissibility matrix K adm = K A K B , the previous admissibility system can be rewritten as:

K adm k, l Y l = F k (27)
Similarly, the residual terms E 2 mod , E 2 int and E 2 mes inherit the following form from Equations ( 16), ( 17) and ( 18):

E 2 mod = (1 -r) E 2 int + r E 2 mes = (1 -r) 2 Y T K int Y + r 2 Y T K mes Y = 1 2 Y T K mod Y (28)
with the notations:

K int i, j = 2 E 2 int (Ψ i , Ψ j ) K mes i, j = 2 E 2 mes (Φ i , Φ j ) K mod i, j = (1 -r) K int i, j + r K mes i, j (29)

Matrix form of the modified problem (19)

The solution Y s of the modified problem is obtained by minimizing the previous quadratic residual E 2 mod under the admissibility constraints (27). This is done by introducing Lagrange multipliers λ, leading to the following linear problem:

K mod K adm T K adm 0 • Y s λ s = 0 F (30)
Due to the choice of the VTCR as the numerical framework, the total number of unknowns remain very reasonable. Let us note that the admissibility matrix K adm is independent of the interface parameters and can be calculated once and for all. This is very interesting economically since the inverse problem is part of an iterative scheme to identify the joint parameters.

Matrix form of the inverse problem (20)

Since the interface matrix K mod depends on the assumed mass and stiffness parameters p = (µ, k) of the interface, so does the solution Y s of the modified problem. Then, the inverse problem consists in finding the optimum parameters which minimize the modified residual used as the cost function:

Find the optimum interface parameters p opt = (µ opt , k opt ) such that:

p opt minimize the cost function E 2 mod (p) = 1 2 Y s (p) T K mod (p) Y s (p)
(31) Currently, this nonlinear step is carried out using classical optimization schemes.

5 Application: identification of joint parameters using numerically simulated experimental data

Numerical example: identification of a joint between two plates

The first example consists of two plates connected by a complex interface. For the sake of simplicity, the joint is assumed to have perfectly rigid Kirchhoff shear behavior, so we need to identify only the joint's bending moment behavior. The joint model used here takes into account a joint stiffness k, and a joint dissipation which is introduced in the form of a joint loss factor denoted by ζ. Note that ζ pertains to the joint damping and one should not confuse with internal damping η of substructures. The aim of this section is to show an application of identification of joint's parameter k with the method which has been introduced in the previous part. The influence of assumed joint damping ζ on the identification results will also be studied. The identification of ζ is not described here since it will be adressed in a forthcoming paper. The exact experimental values of the joint's parameters are given in Table 2. For example, let us consider a fixed frequency equal to a resonance frequency 155 Hz which leads to several wavelengths along the characteristic dimension of the structure. The first plate is subjected to prescribed Kirchhoff shear while the second plate is clamped rigidly (see Figure 5). The experimental data were obtained by VTCR simulation. In order to get reliable experimental information, we chose to simulate field measurements (as opposed to point measurements) using the measurement mesh represented on the figure. Indeed, in the medium-frequency range, classical measurement techniques (essentially point measurements) would not be useful because they would give only very local information. This type of measurement usually implies intrusive techniques which would modify the local behavior of the structure. (At high frequencies, the response is very sensitive to additional masses or stiffnesses.) Also, the structure is very sensitive to various parameters (not only joint parameters, but also geometric parameters): therefore, a local representation would provide unreliable information which could be polluted by very local modes likely to appear in this frequency range. On the contrary, field measurements can lead to very rich experimental information using nonintrusive tools such as laser vibrometer. VTCR simulations were carried out using 100 degrees of freedom per Fig. 7 Cost function, interface error and measurement term as a function of the interface stiffness, with exact assumed interface damping and ± 20 % noise plate: 64 interior complex rays, 9 edge complex rays per edge, and 0 corner complex rays. The interface unknowns were discretized in the same way. The final number of DOFs in the linear system (30) was 400. The parameter r which balances the two terms of the cost function was set equal to the classical value 0.8.

The first result represents an ideal case in which the measurements are built without noise and the assumed interface damping is equal to the exact value. Figure 6 shows the evolution of the residual E 2 mod as a function of Modified residual the joint stiffness parameter k. The vertical dotted line corresponds to the actual stiffness k ex . This curve shows that the minimization process easily leads to the exact value of the interface parameter k ex . The interface error E 2 int and the measurement term E 2 mes are also shown. We can observe that the latter is numerically zero for the exact value k ex due to the absence of noise. However, the interface error is minimum for k ex , but is not zero: this minimum value represents the approximation error due to the VTCR numerical framework. The discretization which was chosen yields an order of magnitude of the error between 10 -3 and 10 -2 .

E 2 mod -ζ exact E 2 mod -ζ +10% E 2 mod -ζ -10% E 2 mod -ζ +50% E 2 mod -ζ -50% E 2 mod -ζ +300% E 2 mod -ζ -95% Exact stiffness
The second calculation was performed in the presence of measurement noise. Each measured value was polluted by a uniformly distributed white noise. The corresponding curves are given in Figure 7. Despite the noise, the minimum of the cost function E 2 mod points to the exact value of k ex . This due to two reasons: the cost function represents a compromise between the experiment and the interface model which regularizes the identification process, and the field measurements induce a sort of averaging process which stabilizes the identification against noise. The only difference from the ideal case is that the measurement term is no longer zero, but represents the noise level injected into the measurements.

The last cases were calculated in the presence of errors in the interface damping coefficient ζ. Figures 8 and9 show the evolution of the cost function with the interface stiffness for a slightly dissipative joint and for a highly dissipative joint. One can observe two types of results: if the joint dissipation is small (9% of the total dissipated energy, Figure 8), the minimum of the cost function is not polluted by the error in the joint damping parameter, even if this error is huge. Conversely, if the joint dissipation is large (68% of the total dissipated energy, Figure 9) it becomes impossible to obtain a proper identified stiffness. This can be explained by the fact that on this particular example both the interface stiffness and the interface damping govern the distribution of energy between the plates. In the presence of an error in joint damping, the identified interface stiffness tends to readjust the energy distribution to counterbalance the effect of the erroneous interface damping factor.

E 2 mod -ζ exact E 2 mod -ζ +10% E 2 mod -ζ -10% E 2 mod -ζ +50% E 2 mod -ζ -50% E 2 mod -ζ +300% E 2 mod -ζ -95% Exact stiffness

Numerical example: identification of four joints connecting four plates

The second structure, shown in Figure 10, consists of four identical plates connected by four different joints. The objective of this section is to check whether our method enables one to find out where the main modeling errors are localized. Indeed, an industrial structure consists of many substructures and joints. In such a case, the number of interface parameters which need to be identified can be very large. A possible strategy, which was already validated for the improvement of finite element models [START_REF] Deraemaeker | Reduced bases for model updating in structural dynamics based on constitutive relation error[END_REF], consists in seeking the most poorly modeled zones. Only the elements in these zones are marked for correction, which reduces the number of parameters in the optimization process. After a first correction step, another localization step is carried out and if other elements show significant errors they are marked and corrected, too, and so on.

Table 3 summarizes the value of each joint parameter for the real structure and for the assumed initial structure. In this example, the interface 3 Experimental (numerically simulated) and assumed initial stiffnesses damping coefficients are set to zero and are not sought for correction. With four joints connecting the substructures, the interface residual is defined as the sum of the residuals E 2

Γi of each joint equation:

E 2 int = α 1 E 2 Γ1 + α 2 E 2 Γ3 + α 3 E 2 Γ4 + α 4 E 2 Γ4 (32)
where α i are weighting coefficients which enable us to increase the confidence level in some joints which are known more precisely (e.g. joints which have already been tested). If all the interfaces are equally poorly known, one takes α i = 1 Nint , N int being the number of joints in the structure. The VTCR calculations required 864 degrees of freedom. Each substructure unknown was discretized into 48 interior complex rays, 13 edge complex rays per edge and 1 corner ray per corner. We chose to underdiscretize the interface unknowns by using 0 interior complex rays, 13 edge complex rays per edge and 1 corner ray per corner. The measurements were field measurements produced with a mesh size equal to 0.1 m. The experimental results were perturbed by ± 10 % uniform white noise.

Table 4 summarizes the results of the identification process. The first residual calculation identified Γ 1 and Γ 2 as being modeled very poorly, while the error in the other two joints was small. The first stiffness was corrected: the optimum value decreased (1.36) and the error on Γ 1 was reduced significantly. Then, Γ 2 was selected as the most erroneous joint and was corrected to a nearly exact stiffness (2.94 instead of 3.00). In the third correction step, the first joint was selected again and, this time, the identified stiffness became virtually exact (0.99 instead of 1.00). After another small correction of Γ 2 , the total residual remained practically constant. This was also the case of the subsequent steps, which are not represented in the table because

Initial

Step 1

Step 2

Step 3

Step 4 log ki E 2

Γ i log ki E 2 Γ i log ki E 2 Γ i log ki E 2 Γ i log ki E 2 Γ i (N.

Example: identification of a joint from actual experimental data

In this final section, we present a first application using actual experimental data. The structure is a steel plate folded in the middle. The folding process introduces a local thickness variation along the fold line. The resulting "joint" was modeled by a line spring whose unknown stiffness we wanted to identify. For the purpose of the measurements, the structure was simply supported by two needles along the joint and loaded by a third needle: see Figure 11 (left). We chose this experimental setup in order to avoid any dissipation other than structural (material and interface). The dissipation due to air movement was neglected since the order of magnitude of the displacements was 10 -5 m. The excitation needle was equipped with an impedance head in order to measure both the force and the acceleration at this point.

Throughout the test, we checked that the force was always positive, thus ensuring that contact was maintained. In order not to perturb the structural response, contact-free measurements of the surface vibration velocity field were carried out using an OMETRON laser vibrometer driven by a PULSE system [25]. Figure 11 (left) shows the measurement mesh, which consisted of 18x18 points distributed equally over the vertical plate. No measurement was taken for the horizontal plate, except at the excitation point. Figure 12 shows the displacement field obtained at 550 Hz, which is a resonance peak (see Figure 11, right). One can observe that a few points experience noise due to the measurement process. At this stage, no attempt was made to smooth the measurement data: the erroneous measurements were used in the same way as those of the other points in the measurement term of the updating process.

Due to the folding process, the connection between the two plates was unknown and was initially assumed to be very stiff (k 0 = 10 5 N.m). In Displacement (m)

Fig. 12 The experimental displacement field at 550 Hz this example, since there were no contact phenomena at the interface between the plates, the joint's damping was not taken into account. Table 5 summarizes the reliable and unreliable information for this identification example. Table 5 Reliable and unreliable information for the inverse problem The simulation produced the results shown in Figure 13 (left) with a residual error of about 37 %. As in the previous examples, the joint's equations were subjected to the updating process, which returned an identified stiffness equal to k opt = 10 2.4 N.m and a final error of about 8.8 %. One should note that no measurement regularization was performed. In fact, regularization came from the field measurements (which generated a large amount of data) and from the balance between the two terms of the cost function. The optimization process enabled us to improve the joint's stiffness: Figure 13 (right) shows the response of the updated structure, which is very close to the experimental response. The cost function is shown in Figure 14.

Continuous model Experimental model

If we had wanted to take damping into account, we could not have identified a damping coefficient at this frequency. This is due to the fact that in this situation the energy dissipated in the joint is negligible compared to the power injected into the structure: most of the energy is dissipated by the plate's material. Therefore, even a slight imprecision in the plate's damping coefficient η given by the manufacturer would have completely polluted the identification of the joint's damping. For a successful joint damping identification, the material damping at this resonance frequency should be known very precisely.

Conclusion

This paper presents the theoretical foundations of an identification method dedicated to the modeling of joints in the medium-frequency range. This approach is based on two main concepts. First, the information is treated according to its level of reliability: the reliable equations are satisfied exactly while the unreliable equations lead to a mixed experimental-numerical residual whose minimum is a compromise among these unreliable equations. Second, the formulation was tailored to the variational theory of complex rays, a numerical framework which presents advantages in the medium-frequency range. This was achieved through the introduction of interface unknowns Updated joint stiffness k opt Fig. 14 The modified residual as a function of the joint's stiffness k during the updating process of the folded plate using a domain decomposition method. The problem's discretized formulation, which we presented in detail, leads to a very reasonably small matrix system. This approach was validated on two numerical examples: the identification of a joint in a folded plate and the identification of four different joints connecting four plates. The identification method leads to an accurate simulated stiffness and appears to be very robust in the presence of measurement noise. When other sources of modeling errors are present in the structure (in our case, the joint's damping), the robustness of the method seems to depend on the influence of this erroneous parameter. The last example was based on actual experimental results. The identification method gave promising results for joint's stiffness, even though the use of more complex joint behavior should lead to an even better response. Identification of joint's damping will be adressed in a forthcoming paper.

Fig. 4

 4 Fig. 4 Locus of admissible k for interior VTCR modes (left); discretized representation of the amplitude U(k) (right)

Fig. 5

 5 Fig. 5 Example of two plates connected by a joint to be identified

First

Fig. 6

 6 Fig.6Cost function, interface error and measurement term as functions of the interface stiffness, with exact assumed interface damping in the absence of noise

Fig. 8

 8 Fig. 8 Cost function as a function of the interface stiffness, with different values of the error in assumed interface damping and a slightly dissipative joint (ζ = 0.1, ± 20 % noise)

Fig. 9

 9 Fig. 9 Cost function as a function of the interface stiffness, with different values of the error in assumed interface damping and a highly dissipative joint (ζ = 10, ± 20 % noise)

Fig. 10

 10 Fig. 10 Example of four plates connected by four joints to be identified log k (N.m) Γ1 Γ2 Γ3 Γ4 Experimental 0.5 4.5 2 1.5 Numerical (start) 2.00 2.00 2.00 2.00Table 3 Experimental (numerically simulated) and assumed initial stiffnesses

Fig. 11

 11 Fig. 11 Experimental setup: field measurement mesh (left) and frequency response function obtained for the velocity averaged in the vertical plate (right)

  Reliable Substructure geometry and material Angular frequency e ω = 550 Hz Lx = Ly = 0.2 m h = 0.0007 m Locations of the measurement points e Xn E = 210 Gpa ν = 0.3 Orientations of the measurement points e θn ρ = 7, 800 kg.m -3 η = 0.0001 Applied force e F d Unreliable Unreliable Interface constitutive law (k) Magnitudes of the displacements e u( e Xn)

5 Fig. 13

 513 Fig. 13 Initial (left) and updated (right) numerical models

Table

	Continuous model	Experimental model
	Reliable	Reliable
	Substructure geometry and material Angular frequency e ω
	Substructure equations:	Locations of the measurement points e Xn
	-boundary conditions	Orientations of the measurement points e θn
	-dynamic equilibrium equation	Locations of the prescribed forces e XF
	-constitutive law	Orientations of the prescribed forces e θF
	Unreliable	Unreliable
	Interface equations	Magnitudes of the displacements e u( e Xn)
	-dynamic equation (µ)	Magnitudes of the prescribed forces e F d
	-constitutive law (k)

Table 2

 2 Experimental (numerically simulated) values of the interface parameters

Table 4

 4 Iterative correction steps for the structure with four joints. The values in bold mark the joint selected for correction at each step no correction was made. It is impossible to obtain better parameters than those produced in the fourth step. With this level of error (0.83 %) the main contribution to the cost function is that of the measurement term, which represents the level of measurement noise in the experimental results.

		m)	% (N.m)	% (N.m)	% (N.m)	% (N.m)	%
	Γ1	2.00 8.61	1.36	1.85	1.36 1.14	0.99	0.21	0.99	0.22
	Γ2	2.00	7.85	2.00 11.57	2.94	0.78	2.94 0.85	3.02 0.77
	Γ3	2.00	0.62	2.00	0.86	2.00	0.23	2.00	0.23	2.00	0.23
	Γ4	2.00	0.67	2.00	0.13	2.00	0.15	2.00	0.11	2.00	0.11
	E 2 mod	2.27 %	1.74 %		0.88 %	0.84 %	0.83 %