Elie Hachem
email: elie.hachem@mines-paristech.fr

Stephanie Feghali
email: stephanie.el_feghali@mines-paristech.fr

Thierry Coupez
email: thierry.coupez@mines-paristech.fr

Stabilized Finite Element Methods vs LES modelling for fluid-strucure interaction with anisotropic adaptive meshing

Keywords: Immersed Stress Method, Stabilized Finite Elements, Fluid-Structure Interaction, Anisotropic Mesh Adaptation, Monolithic Approach

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

Fluid Structure Interaction (FSI) is of great relevance in many fields of engineering as well as in the applied sciences and material forming with applications ranging from bioengineering to aerodynamics and from civil engineering to automotive. Often, when interaction effects are essential this comes along with large/small structural deformations and/or with turbulent flows. However, many available approaches may lack robustness especially in such severe situations. The components in all engineering fields are continuously pushed towards higher performance by seeking new developments that must be able to deal with different situations and regimes. In particularly it must be able to treat encountered problems ranging from the mesh adaptation issues to the coupling engines between different codes, and from small/large deformations to low/high Reynolds numbers flows.

Most of the commercial software packages solve FSI problems using an Arbitrary Lagrangian Eulerian (ALE) formulation [START_REF] Hirt | An Arbitrary Lagrangian-Eulerian computing method for all speeds[END_REF][START_REF] Hughes | Lagrangian-Eulerian finite element formulation for incompressible viscous flow[END_REF][START_REF] Donea | An Arbitrary Lagrangian-Eulerian finite element method for transient dynamic fluid-structure interactions[END_REF]. The solid domain is treated with a Lagrangian formulation. The nodes belonging to the interface between the solid and the fluid are moved with the solid. The displacement of the nodes in the fluid domain do not depend on the fluid motion, but only ensures the continuity between the fluid and the solid domain, and a good mesh quality. ALE methods are robust and accurate, and do not need any extra degrees of freedom. However, important problems arise if the deformations, displacements and rotations of the solid becomes very important [START_REF] Benson | An efficient, accurate, simple ALE method for nonlinear finite element programs[END_REF][START_REF] Souli | ALE formulation for fluid structure interaction problems[END_REF][START_REF] Van Loon | Comparison of Various fluid-structure interaction methods for deformable bodies[END_REF].

A higher popularity has been gained recently by partitioned approaches which allow the use specific solver for each domain. The difficulty remains in transfering the information between the codes. The coupling between the two phases can be enforced using different schemes : weakly or strongly coupled version. The former approach manages with just one solution of either field per time step but consequently lack accurate fulfilment of the coupling conditions. The latter requires sub-iterations. The predominant approach consists in solving the problem iteratively, using fixed-point schemes [START_REF] Le Tallec | Fluid structure interaction with large structural displacements[END_REF] or Newton Krylov methods [START_REF] Fernández | A Newton method using exact jacobians for solving fluid-structure coupling[END_REF][START_REF] Gerbeau | A quasi-Newton algorithm based on a reduced model for fluid structure interaction problems in blood flow[END_REF][START_REF] Michler | An interface NewtonKrylov solver for fluidstructure interaction[END_REF][START_REF] Gerbeau | Fluid structure interaction in blood flows on geometries coming from medical imaging[END_REF]. Actually, the fixed-point methods with dynamic relaxation seem to be the most interesting variant [START_REF] Küttler | Fixed-point fluidstructure interaction solvers with dynamic relaxation[END_REF]. This approach allows the use of fluid and solid solvers for each of the two phases. It is accurate and quite efficient but present an inherent instability depending on the ratio of the densities and the geometry of the domain [START_REF] Caussin | Added-mass effect in the design of partitioned algorithms for fluidstructure problems[END_REF]. As a result, the numerical cost increase drastically and coupling algorithms may not converge. For 3D problems, such difficulties become even more severe.

Monolithic approaches have been proposed to overcome these drawbacks. The whole domain (composed by fluid and solid phases) is considered as a single one, meshed by a single grid, and solved with an Eulerian framework. The continuity at the interface is then obtained naturally and there is no need to enforce it, as it was the case in partitioned methods. If the multi-mesh approaches permit the use of classical fluid and solid solvers, monolithic approaches impose the use of an appropriated unique constitutive equation describing both the fluid and the solid domain. Interface tracking, between the two different domains, can be completed by Immersed Boundary (IB) methods [START_REF] Peskin | The immersed boundary method[END_REF] where the interface is convected on a Lagrangian way. Other methods such as the fictitious domain [START_REF] Peskin | The immersed boundary method[END_REF][START_REF] Glowinski | A distributed lagrange multiplier/fictious domain methods for particulate flows[END_REF] treat the coupling between the domains by applying a constraints across the rigid body using a Lagrange multiplier.

Here in this work, a new monolithic method is developped : the Immersed Stress Method. This method can be seen as an extension of the Immersed Volume Method (IVM) [START_REF] Hachem | Stabilized finite element solution to handle complex heat and fluid flows in industrial furnace using the immersed volume method[END_REF] to treat real fluid-structure interactions. The motivation of pursuing such general approach comes from the desire of not solving two equations, e.g. one in the solid and another in the fluid, where in some cases ; one may still need to provide the boundary conditions between the two domains. Recall also that the complexity to ensure such conditions is amplified when simulating turbulent fluid structure interactions. When dealing with a large diversity of shapes, dimensions and physical properties of structures, such simulations become rapidly very costly, time consuming and limited. A complete description and details about the immersed volume method but used for a different context (developed previously and applied to thermal couplings) is given in [START_REF] Hachem | Stabilized Finite Element Method for Heat Transfer and Turbulent Flows inside Industrial Furnaces[END_REF][START_REF] Hachem | Stabilized finite element solution to handle complex heat and fluid flows in industrial furnace using the immersed volume method[END_REF].

Therefore, we retain the use of a monolithic formulation for fluid/solid and coupling it to some additional features for accurate resolution, in particular at the interface. The monolithic approach in here is made of a unique mesh in which the different domains are taking place by the level set function. Consequently, different strcutures are immersed in a larger domain of different material properties so that boundary conditions at the interface can be replaced naturally.

The second important ingredient of the approach is the use of anisotropic mesh adaptation [START_REF] Coupez | Metric construction by length distribution tensor and edge based error for anisotropic adaptive meshing[END_REF][START_REF] Coupez | Génération de maillage et adaptation de maillage par optimisation locale[END_REF][START_REF] Gruau | Coupez, 3D tetrahedral, unstructured and anisotropic mesh generation with adaptation to natural and multidomain metric[END_REF][START_REF] Hachem | Stabilized finite element solution to handle complex heat and fluid flows in industrial furnace using the immersed volume method[END_REF] at the interface between two different materials. The idea is to apply a fast mesh generation algorithm that allows the creation of meshes with extremely anisotropic elements stretched along the interface, which is an important requirement for FSI problems having internal/boundary layers. It is successfully applied for fixed and some moving objects [START_REF] Hachem | Stabilized Finite Element Method for Heat Transfer and Turbulent Flows inside Industrial Furnaces[END_REF][START_REF] Hachem | Stabilized finite element solution to handle complex heat and fluid flows in industrial furnace using the immersed volume method[END_REF].

The last and most important ingredient focuses more on the finite element solver : on modeling the interaction between the fluid (laminar or turbulent) and the structure in question (rigid, elastic, viscoelastic, etc). For FSI simulations of elastic/rigid body immersed for instance in an incompressible fluid, the global behavior is described by the classical Navier-Stokes equations, with an extra stress tensor [START_REF] Feghali | Monolithic stabilized finite element method for rigid body motions in the incompressible Navier-Stokes flow, submitted to[END_REF]. For instantce, we simulate a rigid solid using the Navier-Stokes solver under constrains to impose the nullity of the deformations. This can be done by simply penalizing the strain rate using a very important viscosity in the solid, which can sometimes be sufficient [START_REF] Ritz | A Numerical Continuous Model for the Hydrodynamics of Fluid Particle Systems[END_REF][START_REF] Hachem | Stabilized Finite Element Method for Heat Transfer and Turbulent Flows inside Industrial Furnaces[END_REF]. It is also possible to enforce directly the nullity of the strain by using an Augmented Lagrange Mulptipliers method [START_REF] Glowinski | A fictitious domain approach to the direct numerical simulation of incompressible viscous flow past moving rigid bodies : application to particulate flow[END_REF][START_REF] Coquerelle | A vortex level set method for the two-way coupling of an incompressible fluid with colliding rigid bodies[END_REF][START_REF] Janela | A penalty method for the simulation of fluid-rigid body interaction[END_REF], solved by an iterative Uzawa algorithm. The problem is then solved by adding an extra-stress tensor comming from the presence of the structure in the fluid. Linear or harmonic mixture laws of the mechanical properties characterizing each domain are then applied at the interface.

Finite Element Formulation 2.1 Governing equations

The governing equations are considered to be three-dimensional and unsteady. Using a monolithic approach, a unique constitutive equation will be solved on the whole domain, with a variation of the parameters depending on the phase that should be modelled. Recall that the the concept of the Immersed Stress Method (ISM) is based on solving the single set of equations by differentiating the subdomains and refining the mesh at this interface using the level set method. The ISM allows the immersion of any strucutre using the level-set function, mixes the physical properties (ρ and η) and finally applies the anisotropic mesh adaptation at the vicinity of the interface. The strong form for the whole domain reads then :

                       ρ(∂ t v + v • ∇v) -∇ • (2η f ε ε ε(v) + τ s -p I d ) = 0 ∇ • v = 0 ε ε ε(u) - 3 2E τ s = 0 ∂ t u + v • ∇u = v + Boundary conditions (1)
where τ s is the extra stress tensor reflecting the presence of the immersed structure (elastic/rigid...) in the incompressible fluid and E is the Young modulus. Note that depending on the value of E, the treated structure will inherit the appropriate law (rigid, elastic, ...).

Stabilized finite element method

Based on a mesh K h of Ω into set of N el elements K, the functional spaces for the velocity, the pressure and the stress are approached by the finite dimensional spaces spanned by V h , P h and T h . The variational multiscale method is used to stabilize the Galerkin formulation and allows the use of equal order continuous interpolations for the velocity and the pressure unkowns (see [START_REF] Hachem | Stabilized finite element method for incompressible flows with high Reynolds number[END_REF] for details). A piecewise constant interpolation for stresses can be used. It consists in here of a decomposition for both the velocity and the pressure fields into coarse/resolved scales and fine/unresolved scales. The distinctive feature of the proposed approach resides in the efficient enrichment of the extra constraint. We first solve and then we substitute the fine-scale solution into the large-scale problem providing additional terms, tuned by a local time-dependent stabilizing parameter, that enhance the stability and accuracy of the standard Galerkin formulation for the transient Navier-Stokes equations. Such approach can deal with laminar/turbulent FSI problems and can handle large/small deformations. Additionally, when higher accuracy is needed, by applying a robust and fast mesh adaptation it provides a much more computational efficiency than coupling solvers.

The enrichment of the functional spaces for the velocity, pressure and stress solutions is performed as follows : V h ⊕V , P h ⊕ P and T h ⊕ T . To this end, v, p, τ will be approximated as :

v = v h + v ∈ V h ⊕ V , p = p h + p ∈ P h ⊕ P , τ = τ h + τ ∈ T h ⊕ T (2)

Numerical scheme

Three equation with three primary variables requires larger computational cost. To circumvent this issue an agmentated Lagragian method and Uzawa's algorithm would be used to solve the system without increasing the size of the linear system. In the same iteration, the problem of non linearity, the time integration and the computation of the Lagrange multiplier would be solved. An implicit time scheme with a Newton method for the non-linear term is used.

LevelSet method

The interface between the phases is resolved using a convected level set approach developed in [START_REF] Ville | Convected level set method for the numerical simulation of fluid buckling[END_REF]. This approach enables first to restrict convection resolution to the neighbourhood of the interface and second to replace the reinitialisation steps by an advective reinitialisation. This enables an efficient resolution and accurate computations of flows even with large density and viscosity differences. The level set function is discretized using a stabilized upwind Petrov-Galerkin method and can be coupled to a direct anisotropic mesh adaptation process enhancing the interface representation.

Anisotropic mesh adaptation

Accurate calculation of the velocities, strains and stresses along the fluid-solid interface is critical for a correct modelling of industrial applications. The difficulty arises due to the discontinuity of the properties of the material across the interface. If this latter is not aligned with the element edges, it may intersect the element arbitrarily such that the accuracy of the finite element approach can be compromised. In order to circumvent this issue, the level-set process is thus coupled to an anisotropic mesh adaptation as described in [START_REF] Coupez | Metric construction by length distribution tensor and edge based error for anisotropic adaptive meshing[END_REF]. The idea of this method is to pre-adapt the mesh at the interface. The mesh becomes locally refined, elements are streched, which enables to sharply define the interface and to save a great number of elements compared to classical isotropic refinement. This anisotropic adaptation is performed by constructing a metric map that allows the mesh size to be imposed in the direction of the distance function gradient. Let us briefly described the main principles of this technique. First of all, one has to resort to a so-called metric which is a symmetric positive defined tensor representing a local base that modify the distance computation, such that :

||x|| M = √ T x • M • x , < x, y > M = T x • M • y . ( 3
)
The metric M can be regarded as a tensor whose eigenvalues are related to the mesh sizes, and whose eigenvectors define the directions for which these sizes are applied. For instance, using the identity tensor, one recovers the usual distances and directions of the Euclidean space. In our case the direction of mesh refinement is given by the unit normal to the interface which corresponds to the gradient of the level-set function : x = ∇α/||∇α||. A default mesh size, or background mesh size, h d is imposed far from the interface and it is reduced as the interface comes closer. A likely choice for the mesh size evolution is the following :

h =      h d if |α(x)| > e/2 2h d (m -1) m e |α(x)| + h d m if |α(x)| ≤ e/2 (4)
Eventually, at the interface, the mesh size is reduced by a factor m with respect to the default value h d . Then this size increases until equalling h d for a distance that corresponds to the half of a given thickness e. The unit normal to the interface x and the mesh size h defined above, lead to the following metric :

M = C (x ⊗ x) + 1 h d I with C =    0 if |α(x)| ≥ e/2 1 h 2 - 1 h 2 d if |α(x)| < e/2 ( 5
)
where I is the identity tensor. This metric returns to isotropic far from the interface (with a mesh size equal to h d for all directions) and to anisotropic near the interface ( with a mesh size equal to h i in the direction x and equal to h d in the others). This method can be assisted by a posteriori anisotropic error estimator, the search of the optimal mesh (metric) that minimizes the error estimator. As a result, an optimal metric as a minimum of an error indicator function and for a given number of elements is obtained. In practice, the mesh is generated in several steps using the MTC mesher and remesher developed by [START_REF] Coupez | Génération de maillage et adaptation de maillage par optimisation locale[END_REF]. The proposed mesh generation algorithm works well for 2D or 3D complex shapes. It allows the creation of meshes with extremely anisotropic elements stretched along the interface. The mesh size is then only refined in the direction of the high physical and mechanical properties gradients. This allow both conserving a high precision in the calculus and in the geometry description, in spite of an important decrease of the total number of degrees of freedom. The grid is furthermore only modified in the vicinity of the interface which keeps the computational work devoted to the grid generation low. Note also that the proposed method can easily handle arbitrary complex geometries. As shown in figure 1 which presents a close-up on the interface zone at the end of the anisotropic adaptation process, the mesh has been gradually refined when approaching the interface. Consequently, only additional nodes are locally added in this region, whereas the rest of domain keeps the same background size. Another method to construct metrics, based on an anisotropic a posteriori error estimator, can also be used with CimLib. This error estimator uses the eigenvalues and eigenvectors of the recovered Hessian matrix for a given function and a fixed number of nodes to construct the metric field [START_REF] Mesri | Advanced parallel computing in material forming with CIMLIB[END_REF][START_REF] Coupez | Metric construction by length distribution tensor and edge based error for anisotropic adaptive meshing[END_REF]. Figure 1 presents a dynamic anisotropic mesh adaptation behind a solid obstacle. As shown, the proposed mesh adaptation approach is capable of handling at the same time a fine mesh around the solid body and the inside the vortices. In this section, we present relatively simple 2D and 3D test cases in order to validate the proposed formulation and to check the accuracy and the efficiency of the immersed stress method. All the numerical simulations were carried out by using the C++ CimLib finite element library .

In figure 3, four solids objects with different densities are falling due to gravity in an air-filled channel. The objective of this test, referred as the Tetris benchmark, is to show the capability of the dynamic mesh adaptation as well as the method to handle high discontinuities of the solids and fluid physical properties. As shown in figure 4, the convection dominated flows of the surrounding air, the four rigid solids movement as well as the representation of the objects are all well taken into account using one a single domain with one set of equations. The final obtained mesh reflects the capability of the method to render a well respected geometry in terms of curvature, angles and complexity. Contrary to others techniques, this promising method can provide an alternative to body-fitted mesh for very complex geometry.

Conclusion

In this paper we have presented a stabilized three-field velocity-pressure-stress, referred as Immersed Stress Method, designed for the computation of rigid bodies in an incompressible Navier-Stokes flow. The presence of the solid is taken into account as an extra stress in the Navier-Stokes equation. The sharp discontinuity of the material properties was captured by an anisotropic refined solid-fluid interface. The robustness of the method to compute the flow and heat transfer with large materials properties differences is demonstrated using stabilized finite element formulations. The approach is applied to the numerical simulation of 2D and 3D test cases. Finally, the capability of the model to simulate the fluid-rigid body interaction was demonstrated.

FIGURE 1 -

 1 FIGURE 1 -Zoom on the fluid-solid interface after anisotropic adaptation

FIGURE 3 -

 3 FIGURE 3 -The Immersed Stress Method with dynamic anisotropic mesh adaptation Figure 5 presents the parallel numerical simulation of turbulent heat transfer inside an industrial furnace using the proposed monolithic fluid-structure approach with fixed anisotropic mesh adaptation. The

FIGURE 4 -

 4 FIGURE 4 -Immersion of four rigid bodies : velocity vectors at different time instants

FIGURE 5 -

 5 FIGURE 5 -Turbulent heat transfer inside an industrial furnace


R Fedele
email: fedele@stru.polimi.it

R Santoro

Presqu'île de Giens (Var) Simultaneous assessment of mechanical properties and boundary conditions by Digital Image Correlation

Keywords: Cohesive crack model, Inverse Problems, Digital Image Correlation

This communication represents a further contribution to the study of full-field identification procedures based on Digital Image Correlation. An augmented identification procedure is presented herein, including as unknowns both the (cohesive) mechanical parameters of the interface and the actual boundary conditions to be prescribed along the border of the Region-Of-Interest. In this way, data smoothing and parameter identification procedures, so far usually performed in a sequence, are tackled simultaneously in a coupled framework.

Introduction

During the last years, high-accuracy kinematic measurements provided by Digital Image Correlation at different observation scales, are being employed for an increasingly number of applications in material mechanics, see e.g. [START_REF] Avril | Overview of identification methods of mechanical parameters based on full-field measurements[END_REF][START_REF] Roux | Digital Image Mechanical Identification (DIMI)[END_REF]. In particular, the availability of full-field data allows one to monitor during non trivial tests small sub-domains on the specimen surface, in which the phenomenon of interest turns out to be predominant. In this way, virtual tests can be performed by prescribing local boundary conditions, which integrate or even supersede the overall boundary conditions prescribed to the whole tested sample (e.g., the reaction force). In this respect, the whole sample with its boundary conditions represents no more the unique reference of material mechanics for at least two kind of considerations. i) The overall boundary conditions prescribed to the sample (no matter the static or kinematic ones) are usually regarded as certain, and they are deterministically prescribed in the numerical simulation of the whole specimen response without taking into proper account the gap between the ideal constraints and their actual response. ii) The utopia of elementary tests, explicitly conceived to provide a purely constitutive response, easily and accurately, is declining. Design of complex laboratory experiments, aiming at generating strain and stress field close to those expected under service, is playing more and more a crucial role.

In a previous paper [START_REF] Fedele | Identification of adhesive properties in GLARE laminates by Digital Image Correlation[END_REF], a numerical-experimental methodology has been developed, to estimate the parameters governing the mixed mode response of an adhesive joint inside an industrial manufact. A local approach was proposed therein, focusing on a small sub-domain, referred in what follows a Region-Of-Interest (ROI), for which (Dirichlet) boundary conditions are prescribed. As an undesired drawback, the stress and strain fields inside the ROI turn out to be hardly legible, due to noise on boundary data provided by DIC at the microscale. In fact, because of the elastic-plastic behavior of the metal adherents, even high frequency fluctuations are expected to propagate from the outer boundary into the sub-domain, and not to be rapidly attenuated as it occurs in linear elasticity. The slender geometry of the monitored sub-domain worsen this situation.

In this study, an improvement of the original strategy is presented. The boundary displacements to be prescribed along the ROI border at different measurement instants are included as unknowns in the identification problem (amounting to a few thousands). The forward operator, based on a finite element discretization, is now regarded as a function of both the arguments, i.e. the cohesive parameters of the adhesive joint and the boundary displacements. Also the cost (or objective) function to minimize, quantifying the discrepancy between DIC measurements and model predictions, includes now displacements both along the boundary and inside the sub-domain. In the other terms, the optimality of the augmented problem solution rests on an extended notion of distance, including data belonging to the closure of the ROI, and not only the interior part of it.

This communication is organized as follows. Section 2 outlines the main features of the proposed experimental tests and their simulations. The novel extended identification problem is presented in Section 3, where some preliminary results are presented. Closing remarks and future prospects are gathered in Section 4.

Debonding experiments and their simulation

The main features of non-conventional tests performed on adhesively bonded assemblies, and their simulations by finite elements, are briefly outlined in what follows. For further details, the interested reader is referred to [START_REF] Fedele | Identification of adhesive properties in GLARE laminates by Digital Image Correlation[END_REF]. The experiments reported herein are applied to the characterization of a structural (60-µm thick) adhesive layer, made of BMS 5-101 produced by 3M (commercial denomination AF163 degree 10), joining a (1.4-mm thick) laminate skin and a Z-shaped reinforcement (stringer), both made of GLARE composites for aerospace applications. Tests were designed and performed on rectangular joined specimens (surface area: 297 × 50 mm 2 ) to generate within the adhesive complex stress and strain states up to a complete debonding under mixed mode. A small region on the specimen surface around the adhesive layer (of area approximately equal to 0.6 × 1.3 mm 2 ) was selected as the Region-of-Interest (ROI) and monitored during the tests. The estimated pixel size amounts to 1.93 µm.

Recourse is made to a Digital Image Correlation (DIC) code based on a Galerkin, finite-element discretization of the displacement field, with Q4-(bilinear) shape functions, see e.g. [START_REF] Roux | Digital Image Mechanical Identification (DIMI)[END_REF]. The kinematic measurements at the micro-scale provided by DIC concern two strips belonging to the joined substrates, close to the adhesive layer. The Aluminum substrates (cold-rolled 2024-T3 alloy) are modeled accordingly to Ramberg-Osgood's model. All parameters governing the response of the aluminum alloy sheets are assumed to be known a priori with a sufficient accuracy, whilst the actual mechanical properties of an adhesive layer within the joined assembly are the objective of the present identification procedure.

To describe the joint response under mixed-mode loading conditions, the exponential relationship proposed by Van den Bosch et alii [START_REF] Van Den Bosch | An improved description of the exponential Xu and Needleman cohesive zone law for mixed-mode decohesion[END_REF] is adopted

T n = φ n δ n ∆ n δ n exp - ∆ n δ n exp - ∆ 2 t δ 2 t ( 1a
)
T t = 2 r φ n δ t ∆ t δ t 1 + ∆ n δ n exp - ∆ n δ n exp - ∆ 2 t δ 2 t (1b)
Normal and tangential components of traction and (relative) displacement vectors are denoted by T n and T t , ∆ n and ∆ t , respectively. Symbols δ n , δ t , φ n and r • φ n designate the characteristic lengths and the works-of-separation (for a unit surface) under pure Mode I (peel) and Mode II (shear), respectively, being r = φ t /φ n the (non-dimensional) ratio between the tangential and normal works-of-separation. In what follows the three unknown parameters governing the adhesive layer behavior are gathered in a vector form x = {φ n , δ n , δ t } T , being r = 4 assumed a priori. The plane-stress model adopted to simulate the mechanical response of the considered sub-domain is made of 570 finite elements, and the adhesive layer is discretized by means of 39 cohesive elements with a finite thickness (see also [START_REF] Fedele | Identification of adhesive properties in GLARE laminates by Digital Image Correlation[END_REF]).

The nonlinear relationship between the unknown parameters and the selected measurable quantities y comp Ω , will be referred to as forward operator, and denoted by the symbol:

y comp Ω = H x u ∂Ω (2)
Herein the two arguments of the above operator are emphasized, namely the cohesive mechanical properties of the interface x = {φ n , δ n , δ t } T , and vector u ∂Ω gathering boundary displacements at different measurement instants i = 1, ..., n t , in what follows both included in a unique augmented parameter vector, say X = {x T u T ∂Ω } T . The measurable quantities y comp Ω are the displacements belonging to the whole ROI, i.e. to the closure Ω ≡

• Ω ∪ ∂Ω, including the interior part and the frontier of the monitored subdomain. In the local approach proposed in [START_REF] Fedele | Identification of adhesive properties in GLARE laminates by Digital Image Correlation[END_REF], displacements u i ∂Ω prescribed along the ROI boundary have to be considered as independent boundary conditions, since they cannot be derived deterministically by any other model. Therefore, the augmented parameter space in the present case turns out to be the direct sum of the space of the nonnegative cohesive parameters, and the one spanned by boundary nodal displacements. Also the image space of the forward operator, or data space (the space of measurable quantities y Ω ), can be regarded as the direct sum of two sets previously mentioned, the one gathering nodal displacements belonging to the interior part of the ROI, i.e. u • Ω , and the analogous set relevant to the border, namely u ∂ Ω .

Extended inverse problem

In the original inverse methodology presented in [START_REF] Fedele | Identification of adhesive properties in GLARE laminates by Digital Image Correlation[END_REF], experimental displacements gathered in vector u exp i ∂Ω provided by DIC, were deterministically prescribed along the outer boundary ∂Ω of the monitored subdomain at different measurement instants i (i = 1, ..., n t ). The cohesive parameters of the joint were then estimated as follows (see [START_REF] Fedele | Identification of adhesive properties in GLARE laminates by Digital Image Correlation[END_REF]):

x = arg min x S • u exp • Ω -u comp • Ω (x) 2 W (3)
It is worth emphasizing that, in Eq. ( 3), the forward operator is a function of the joint parameters only, and the measurable quantities located closely to the joint, inside the ROI, were selected therein via the Boolean matrix S.

In this communication an extended inverse methodology is developed, in which both the cohesive parameters x = {φ n , δ n , δ t } T , and the "actual" boundary displacements at different instants, i.e. u i ∂Ω , i = 1, ..., n t ≡ 19, are assumed as unknowns to be estimated. The optimal parameter vector X, is obtained by minimizing a suitable discrepancy norm (objective or cost function). Such a function quantifies the discrepancy between measured quantities (e.g., provided by DIC) over the whole problem domain Ω =

• Ω ∪ ∂ Ω ,and their counterparts computed via Equation [START_REF] Van Den Bosch | An improved description of the exponential Xu and Needleman cohesive zone law for mixed-mode decohesion[END_REF]. Of course, the measurable quantities computed by the forward operator inside the ROI, or prescribed along the outer boundary are expected to be close to the ones provided by DIC, but not exactly coincident with them, due to the presence of noise and model error. The extended identification problem can be straightforwardly formulated in a least-square sense as follows:

X = { x u ∂Ω } = arg min { x u ∂Ω } ω x u ∂Ω ω (X) = u exp • Ω -u comp • Ω x u ∂Ω 2 W + u exp ∂Ω -u ∂Ω 2 W + c x • R x (x) + c u • R u (u ∂Ω ) (4)
In Equation ( 4), R x (x) and R u (u ∂Ω ) denote non-negative (≥ 0) regularization terms, concerning cohesive parameters and boundary displacements, and weighted by non-negative penalty factors c x and c u , respectively. The regularization terms are identically null inside the relevant feasible domains, whilst they increase fastly (and smoothly) moving away from it. Diagonal weighting matrices W scale and make non-dimensional the kinematic residuals. As a meaningful numerical exercise developed in a Matlab environment [START_REF]Matlab 6.5, Optimization Toolbox Manual[END_REF], Figures 1 and2 show the vertical displacement field and the stress component σ 11 distribution predicted by the calibrated finite element model at a generic instant of the test simulation, with boundary data provided by DIC (Eq. 3) in a), or estimated by the present extended procedure (Eq. 4) in b).

Closing remarks

In this communication an extended inverse methodology resting on a local approach has been presented, including as unknowns both mechanical parameters and boundary conditions. As a consequence, the parameter space exhibits an unusually high dimensionality. The large-scale inverse problem above defined turns out to be severely ill-posed, and recourse is made to regularizing provisions in a Tychonoff sense. In particular, a penalty approach is used to preserve the constraint-free minimization scheme, with reference to two regularization terms: i) a Neumann (equilibrium) condition has been considered for the traction distributions along the interface, enforced to vanish where a macroscopic (traction-free) crack is observed; ii) the boundary displacements are constrained to belong to a finite-dimensional subspace, spanned by Cebycëv polynomials. As a possible further development, to drastically reduce the problem dimensionality, sequential decoupling will be investigated, between the identification of boundary displacements and that of cohesive parameters, and, for the former, between boundary data relevant to different measurement instants. The preliminary results achieved by the present approach turn out to be rather satisfactory, in terms of both traction distribution along the interface and strain and stress fields inside the ROI predicted by the FE model.

Figure 1 :

 1 Figure 1: Vertical displacement field [mm], predicted by the calibrated FE model inside the ROI (dimensions in mm) at measurement instant t = 7 during the test simulation, being prescribed along the boundary: in a), the experimental (noisy) displacements provided by DIC through Eq. 3; in b), the "actual" displacements identified by the present inverse methodology Eq. 4.

Figure 2 :

 2 Figure 2: Distribution of stress component σ 11 [MPa] inside the ROI at measurement instant t = 7 during the test simulation. Conventions are the same of Figure 1.


M Cremonesi
email: cremonesi@lmt.ens-cachan.fr

F Frangi

U Perego

Presqu'île de Giens (Var) A Lagrangian finite element method for non-Newtonian free-surface fluid flows and fluid-structure interaction problems

Keywords: Lagrangian fluid, free-surface, fluid-structure interaction, non-Newtonian fluids

A Lagrangian Finite Element method is applied to the discretization of the Navier-Stokes equations to simulate fluid problems with significant evolution of the free-surfaces. The same approach is also used to treat fluid-structure interaction problems. Both Newtonian and non-Newtonian fluid behaviours are considered. Different examples are performed to show the potentialities of the approach.

Introduction

A critical point in the simulation of the fluid flow problems is the treatment of the free surface. Similarly, in fluid-structure interaction problems the efficent determination of interface between solid and fluid, is always crucial. The Arbitrary Lagrangian Eulerian method (ALE), in which the movement of the fluid particles is independent of that of the mesh nodes, is often used to solve this kind of problems. Other possible solution schemes are based on the volume of fluid or level set algorithms. Meshfree and meshless methods are often used for their ability to recover the free-surfaces and the interfaces. In this work an other possibility to overcome the difficulties concerning the tracking of the interfaces is presented. A Lagrangian approach is used to describe both the fluid and the structural phases. In particular, a staggered scheme is considered in which the fluid is treated using a new implementation [START_REF] Cremonesi | A lagrangian finite element approach for the analysis of fuid-structure interaction problems[END_REF] of the Particle Finite Element Method (PFEM), first proposed by Oñate, Idelsohn and coworkes [START_REF] Oñate | The particle finite element method. An overview[END_REF][START_REF] Idelsohn | The particle finite element method : a powerful tool to solve incompressible flows with free-surfaces and breaking waves[END_REF], and the structure using a classical finite element method.

Using a Lagrangian approach for the fluid flow the convective terms in the momentum conservation equation, typical difficulty of the Eulerian framework, disappear. The difficulty is however transferred to the necessity of a frequent remeshing of the fluid domain. In fact, in a Lagrangian framework the position of element nodes is updated as a consequence of the fluid flow and the element distortion can become excessive. To avoid these distortions, a possible remedy consists of a systematic remeshing of the volume of the problem. To this purpose, an efficient Delaunay triangulation has been adopted. Moreover, to define the integration domain and to correctly impose the boundary conditions, an "alpha shape method" [START_REF] Edelsbrunner | Three dimensional alpha shapes[END_REF] is used.

The proposed Lagrangian Finite Element Method is particularly suited also for the solution of fluidstructure interaction problems in the presence of free-surfaces, in conjunction with a classical finite element method for the solid part.

Governing equations

The Particle Finite Element Method (PFEM) [START_REF] Oñate | The particle finite element method. An overview[END_REF][START_REF] Idelsohn | The particle finite element method : a powerful tool to solve incompressible flows with free-surfaces and breaking waves[END_REF] is conceived for the solution of fluid-dynamics problems using a Lagrangian approach based on the formulation of the Navier-Stokes equations in material coordinates. In the Lagrangian framework the equations of motion of a fluid can be written as :

ρ 0 Du Dt = DivΠ + ρ 0 b in Ω 0 × (0, T ) (1)
Div

JF -1 u = 0 in Ω 0 × (0, T ) (2)
where ρ 0 is the density of the fluid, u is the velocity, Π = JσF -T is the first Piola-Kirchhoff stress tensor, σ is the Cauchy stress tensor, F is the deformation gradient, J is the determinant of F and Ω 0 represents the initial (and reference) configuration. The Cauchy stress tensor σ can be related to the velocity u and to the pressure p by the relation :

σ = -pI + τ(u) (3)
where τ(u) represents the deviatoric stress. If a Newtonian constitutive law is considered, the deviatoric stress can be directly expressed as a function of the velocity :

τ = µ gradu + gradu T ( 4
)
where µ is the viscosity of the fluid. If, insted, a non-Newtonian behaviour is considered, the deviatoric stress is a generic nonlinear function of the velocity.

Space discretization

To discretize in space the problem (1)-( 2) an isoparametric finite element discretization is introduced. As explained in the next section, the same order of interpolation is used for both velocity and pressure unknows. The semi-discretized form of the equations ( 1)-( 2) reads :

M(x) DU Dt + K(x)V + D T (x)P = B (5) D(x)U = 0 ( 6
)
where M is the mass matrix, K is the fluid stiffness matrix, D is the discretization of the divergence operator and B is the vector of external forces :

M ab = Ω 0 ρ 0 N a N b IdΩ 0 D ab = Ω 0 JN a Grad(N b )F -1 dΩ 0 K ab = Ω 0 µJ Grad(N a )F -1 ⊗ Grad(N b )F -1 dΩ 0 B a = Ω 0 ρ 0 bN a d Ω 0
where I is the identity matrix and N a the nodal shape function. The matrix operators K and D clearly depend on the current configuration through the deformation gradient F.

Pressure stabilization

In the present approach every time the mesh is regenerated (connecting the points with the Delauany triangulation), the informations related to the element are lost. The same occurs for informations related to nodes not coinciding with triangle vertices. Consequently, to avoid convecting data from the old to the new mesh, only vertx node data can be stored, and only linear shape functions can be used in the space discretization of the problem (1)-(2). To avoid spurious oscillations due to the fact that the inf-sup condition is not satisfied [START_REF] Brezzi | Mixed and Hybrid Finite Element Method[END_REF], a stabilization must be introduced. In particular a pressure-stabilizing/Petrov-Galerkin (PSPG) stabilization is adopted [START_REF] Tezduyar | Incompressible flow computations with stabilized bilinear and linear equal-order-interpolation velocity-pressure elements[END_REF]. Introducing the stabilization in the finite element space discretization previously considered, the semi-discrete stabilized equations read :

M(x) DU Dt + K(x)V + D T (x)P = B (7) C(x) DU Dt + D(x)V + L(x)P = H(x) (8)
where the matrices C, L and H are given by :

C ab = N el ∑ e=1 Ω e 0 τ e pspg Grad(N a )N b dΩ 0 L ab = N el ∑ e=1 Ω e 0 τ e pspg ρ Grad(N a ) • Grad(N b )F -1 dΩ 0 H a = N el ∑ e=1 Ω e 0 τ e pspg Grad(N a ) • bdΩ 0
and the stabilization parameter is defined as :

τ e pspg = z e 2 u (9)
and z e is the "element length" defined to be equal to the diameter of the circle which is area-equivalent to the element e.

Time integration

The time integration is performed with a classical backward Euler scheme. Introducing a partition of the time domain into N time steps of the same length ∆t, and choosing as the reference configuration at time t n , the discretized system reads :

1 ∆t M U n+1 -U n + K(U n+1 )U n+1 + D T (U n+1 )P n+1 = B (10) 1 ∆t C U n+1 -U n + D(U n+1 )U n+1 + L(U n+1 )P n+1 = H (11)
where the matrices K, D and L and the vector H depend nonlinearly on the vector U n+1 .

Description of the method

The PFEM, in the form implemented in the present work, consists of the following conceptual steps (see also [START_REF] Oñate | The particle finite element method. An overview[END_REF] and reference therein). Figure 1 shows a pictorial representation of the steps.

1. Fill the fluid domain with a set of points referred to as "particles" (Figure 1(a)). The accuracy of the numerical solution is clearly dependent on the considered number of particles.

2. Generate a finite element mesh using the particles as nodes. This is achieved using a Delaunay triangulation (Figure 1(b)).

3. Identify the external and internal boundaries to compute the domain integrals and to impose correctly the boundary conditions (Figure 1(c)).

4. Solve the non-linear Lagrangian equations of motion ( 10)-(11) finding velocity and pressure at every node of the mesh.

5. Update the particle positions using the computed values of velocity and pressure (Figure 1(d)).

Check the mesh distortion :

-if the mesh is too distorted go back to step 2 and repeat for the next time step.

-if the mesh is not too distorted go back to step 3 and repeat for the next time step.

Mesh generation and identification of external boundaries

To solve the Lagrangian Navier-Stokes equations in a domain defined by the particles, a finite element mesh is necessary. A possibility to generate a mesh starting from a particle distribution is the Delaunay triangulation. The Delaunay triangulation generates the convex hull of the set of particles, i.e. the convex figure with minimum area that encloses all the particles belonging to the set. -identification of the real shape of the particle distribution ; -identification of the nodes that belong to the boundary. The second phase is almost trivial, but the first one is critical. As introduced, the Delaunay triangulation generates the convex hull of the set of particles. However, the convex hull may be not conformal with the external boundaries. To clarify this problem, a set of points is shown in Figure 2 A possibility to overcome this problem is to correct the generated mesh using the so called α-shape method [START_REF] Edelsbrunner | Three dimensional alpha shapes[END_REF] as proposed in [START_REF] Idelsohn | A lagrangian meshless finite element method applied to fluidstructure interaction problems[END_REF]. The key idea is to remove the unnecessary triangles from the mesh using a criterion based on the mesh distortion. For each triangle e of the mesh, the minimal distance h e between two nodes in the element and the radius R e of the circumcircle of the element are defined. If h is computed as the mean value of all the h e , the shape factor :

α e = R e h ≥ 1 √ 3 , (12)
is an index of the element distortion. It is worth recalling that R e /h e = 1/ √ 3 is the ratio for an equilateral triangle. All the elements that do not satisfy the condition :

α e ≤ ᾱ ( 13
)
are removed from the mesh, where ᾱ ≥ 1 is assumed. Figure 2(c) shows the Delaunay triangulation coupled with the alpha-shape scheme for the example introduced before.

In general, increasing the value of ᾱ, fewer triangles are removed from the original mesh and, for ᾱ → ∞, the original Delaunay triangulation is always recovered.

Once the unnecessary triangles are removed from the Delaunay triangulation, the particles which actually belong to the boundary can be identified.

The α-shape method can also be used for the identification of the fluid particles which separates from the rest of the domain. When a particle on the boundary is a vertex of only one triangle and the triangle is overly distorted, the particle is separated from the domain and all the mass of the triangle is assigned to this particle in order to preserve the total fluid mass. After separation, the motion of these particles is governed by the body force and the initial velocity which they are subjected to. At each new Delaunay triangulation, the separated particles become vertices of triangles. In the case that a separated particle has approached enough the fluid free boundary, the triangle is not eliminated by the α-shape method check and the particle is again incorporated in the fluid mass.

Fluid-structure interaction

The proposed Lagrangian Finite Element Method is particularly suited for the solution of fluidstructure interaction problems in the presence of free-surfaces, in conjunction with a classical finite element method for the solid part. A critical issue of fluid-structure interaction schemes is the identification of the contact interfaces between the solid and fluid domains. The evolution of the interaction surfaces is tracked using a novel algorithm which exploits the features of the Lagrangian approach based on the continuous remeshing introduced for the fluid solution. The proposed algorithm is based on the superposition of a set of fictitious fluid particles (called "ghost particles") to the nodes of the solid domain, which can come in contact with the fluid domain. When the Delaunay triangulation is performed, the alpha-shape criterion selects those parts of the interface where the fluid particles can possibly come into contact with the structure. If the two discretized domains are not in contact, the fluid and the solid analyses are performed separately without any interaction. If instead parts of the two discretized domains are in contact, a coupled analysis is performed with a Dirichlet-Neumann iterative approach. The fundamental steps of the interaction can be listed as follows.

-Superpose a set of ghost particles to the solid surfaces which can possibly come into contact with the fluid domain. -Perform the Delaunay triangulation.

-Apply the α-shape method to remove the unnecessary triangles.

-Check whether the fluid is in contact with the solid part. The fluid element is assumed to be in contact with the solid if at least one of its nodes is a ghost particle.

-If the two discretized domains are not in contact, the fluid and the solid analyses are performed separately without any interaction. -If the two discretized domains are in contact, a coupled analysis is necessary.

Applications

The proposed scheme has been applied to both Newtonian and non-Newtonian fluids. Comparisons with numerical benchmarks and with experimental results are presented to show its potentialities.

5.1 Fluid-structure interaction : filling of a deformable container

The present approach has been used to simulate an example of a fluid-structure interaction problem involving large fluid motions with large displacements of the structure. Under the action of the gravity force a fluid drops down from a funnel-shaped rigid container into a thin elastic container. The elastic properties of the container are given as : Young modulus E = 210 KPa, Poisson ratio ν = 0.3 and density ρ s = 2000 Kg/m 3 . For the fluid a fluid viscosity µ = 0.001 Kg/ms and a density ρ f = 1000 Kg/m 3 are adopted. In the solution an initial fluid mesh of 5749 nodes and 11103 elements and a solid mesh of 2188 nodes and 3365 quadrilateral plane strain elements are adopted and a time step of 0.001 s is used. Figure 3 shows the deformation of the elastic container at different time step.

Non-Newtonian fluids : slump flow test

The non-Newtonian solver has been used to simulate typical tests on fresh concrete, mortar and cement pastes. In general, yield stress and viscosity are parameters used to characterize the workability of this kind of materials. The experimental determination of these parameters is a research topic of recent years. However, in the identification of parameters an efficient numerical tool can be also very useful. To this purpose, the non-Newtonian solver is used to simulate typical tests widely used in the experimental identification of the rheological properties of fresh cement paste and concrete. For the cement pastes and mortars a Bingham constitutive behaviour, which is the simplest model encompassing the concept of yield stress, is considered. Using the Papanastasiou approximation [START_REF] Papanastasiou | Flows of materials with yield[END_REF], the Bingham behaviour can be described as a relationship between deviatoric stress and velocity :

τ = μ ε(u) (14) μ = µ + τ 0 ε 1 -e -n ε (15)
where ε(u) = gradu + gradu T is the symmetric part of the gradient of the velocity. Figure 4 shows the simulation of a typical slump flow test. A rigid cone, open at both ends, is filled in with fresh cement paste or concrete, then the cone is suddenly raised allowing the fluid concrete to flow onto a board. The material will move only if the yield stress is exceeded and will stop when the stress is below the yield stress. Table 1 shows a comparison between the experimental radius and the radius computed with the present approach for different values of density, viscosity and yield stress for cement pastes. A good agreement can be observed between experimental and numerical results. Details about the numerical simulation of these tests can be found in [START_REF] Cremonesi | Simulation of the flow of fresh cementitious suspensions by a lagrangian finite element approach[END_REF].

Deformable landslide simulation

The proposed approach is then used to simulate the wave generated by a deformable mass of granular material sliding in a water channel. This example is particularly interesting because impulse waves in bays, lakes or reservoirs generated by landslides or rock falls may propagate with possibly catastrophic consequences for the surronding enviroment.

For the numerical simulation of the sliding mass a non-Newtonian Bingham fluid behaviour is assumed, while a classical Newtonian behaviour is considered for the water [START_REF] Cremonesi | A lagrangian finite element approach for the simulation of waterwaves induced by landslides[END_REF]. Figure 5 shows snapshots of the analysis at different time steps.

  FIGURE 1 -Steps of Particle Finite Element Method

  (a) and its Delaunay triangulation in Figure 2(b). It is clear that the Delaunay triangulation does not match the real boundaries.

FIGURE 2 -

 2 FIGURE 2 -(a) Distrubution of points, (b) Delaunay triangulation, (c) Delaunay triangulation with αshape.

FIGURE 3 -

 3 FIGURE 3 -Deformation of an elastic container. Snapshots at different time step.

FIGURE 4 -

 4 FIGURE 4 -Mini-slump-flow test. Snapshots at different time steps for a cement paste.

FIGURE 5 -

 5 FIGURE 5 -Landslide without initial velocity. Snapshots at different time steps.

TABLE 1 -

 1 Mini-slump-flow test : comparisons for varying density, viscosity and yield stress for cement pastes.


N Valoroso
email: nunziante.valoroso@uniparthenope.it

Giens (Var) On the use of drilling rotations for DIC-based identification of interface parameters in adhesive structural joints

Keywords: Adhesive joints, cohesive-zone models, identification, digital image correlation

The problem of parameter identification via inverse analysis based on full-field kinematic data obtained via Digital Image Correlation (DIC) is discussed with special attention on resolving the surface motion of adhesive joints. The response up to failure of a DCB-like structure is simulated using either classical continuum and interface elements or membrane elements with drilling dofs and interface elements with drilling rotations along the bondline. Identification results are presented that show the differences between two such cases in terms of performances and parameter estimates.

Introduction

In recent years adhesive bonding technologies have become widespread and almost commonplace in industry. Whether it be electronic components [START_REF] Li | Recent advances of conductive adhesives as a lead-free alternative in electronic packaging : Materials, processing, reliability and applications[END_REF] or furniture, food packaging materials [START_REF] Choia | Development of biodegradable hot-melt adhesive based on polye-caprolactone and soy protein isolate for food packaging system[END_REF] or floor coverings, shoes or high-tech solutions [START_REF] Savage | Failure prevention in bonded joints on primary load bearing structures[END_REF], there is almost no product on the market today that is produced without adhesives. This virtually unlimited range of applications calls for an intensive research work in order to fully exploit the potential of adhesive joints.

As for most structural components consisting of the assembly of individual elements, failure of adhesive joints due to damage growth at bonded interfaces, leading to fracture development by sliding and separation, is one of the most important failure modes and for its simulation the cohesive zone concept initially proposed by Barenblatt [4] has become increasingly popular in the last decades. One of the reasons of its success is probably the flexibility of the cohesive approach to fracture, which provides an effective phenomenological description of the complex microscopic processes leading to the progressive decay of cohesive forces and the formation of traction-free surfaces.

For situations where interface positions are a priori known, cohesive models are used in conjunction with zero-thickness interface elements. These have to adequately resolve the process zone and be compatible with the surrounding continuum elements in order to obtain meaningful answers. Moreover, the availability of trustworthy estimates of material parameters that characterize the traction-relative displacement relationship and the general coherence between the data reduction schemes and the FE model in which the cohesive law is used are not secondary aspects, since the assumptions made for computing the material parameters from experiments have a direct impact on the results of computations.

In this study the fundamental problem of parameter identification via inverse analysis based on fullfield kinematic data obtained via Digital Image Correlation (DIC) is discussed with special attention on resolving the surface motion of adhesive joints up to failure. This is not an easy task since computational failure analysis would require elements that are highly flexibile and accurate, and with reduced sensitivity to mesh distortion, whilst for identification purposes elments with low number of nodes are quite desirable since they result in low connectivity of the structural stiffness matrix and reduced computational cost. Obviously, finding the optimal trade-off between accuracy and computational cost is critical.

The evaluation of mode-I parameters governing the damage mechanics-based formulation developed in [START_REF] Valoroso | A damage-mechanics-based approach for modelling decohesion in adhesively bonded assemblies[END_REF] is addressed as the solution of a nonlinear programming problem. A least-squares norm is used as objective function that quantifies the distance between experimental data and the analogous quantities computed via finite elements as a function of the unknown parameters. The data set concern the deformation process of of a suitable region of interest extracted from the tested Double Cantilever Beam (DCB) specimen [START_REF] Valoroso | Characterization of a cohesive-zone model describing damage and decohesion at bonded interfaces. Sensitivity analysis and mode-I parameter identification[END_REF].

For the DCB problem at hand we suggest to simulate the response up to failure using interface elements with drilling rotations located along the bondline [START_REF] Valoroso | A novel interface element with rotational dofs[END_REF]. The interest of such elements is motivated by the fact that thay can be coupled with membrane/shell elements with in-plane rotations, which are well-known to exhibit excellent accuracy and convergence features both for regular and distorted meshes, see e.g. [START_REF] Hughes | Numerical assessment of some membrane elements with drilling degrees of freedom[END_REF]. The interface elements with drilling dofs have been implemented in a customized version of the finite element code FEAP [START_REF] Taylor | FEAP, User and Programmer Manuals[END_REF] along with a suitable procedure for the computation of gradients of displacements and forces with respect to the material parameters, i.e. parameter sensitivities. These are indeed a key ingredient for identification purposes since, one one hand, thay may allow to determine the location of measurements with the highest information content during experiments and, on the other side, can also provide the search direction in a gradient-based minimization algorithm.

The outline of the paper is as follows. In section 2 the cohesive model is briefly introduced; the interface element formulation is then recalled in section 3. The DCB test and the inverse identification problem are discussed in section 4 ; the results of identification based on truly experimental results in terms of full-field displacements data obtained via DIC are finally presented in section 5.

Cohesive damage model

The formulation of a cohesive-like model basically requires the definition of a relationship between interface tractions and displacement discontinuities and the introduction of a criterion for damage to grow and the process zone to move. In particular, we recall in the following the model proposed by Valoroso and Champaney [START_REF] Valoroso | A damage-mechanics-based approach for modelling decohesion in adhesively bonded assemblies[END_REF] and consider the simplest two-parameter version for the mode I case, that is governed by the following equations :

t = (1 -D)k [[u]] + + k -[[u]] - Y = 1 2 k [[u]] 2 + φ = Y -Y * ≤ 0 Ẏ * = Ḋ ∂F ∂D φ ≤ 0; Ḋ ≥ 0; Ḋ φ = 0 (1)
where t and Y respectively denote the interface traction and the damage-driving force, D ∈ [0, 1] is the scalar damage variable, [[u]] is the displacement jump in the direction normal to the interface while k and k -are the undamaged interface stiffnesses in tension and compression, respectively. The impenetrability constraint is introduced in penalty form via the stiffness coefficient k -and by distinguishing between the positive and negative part of the displacement jump. Figure 1 -Cohesive relationship to be identified, (see also [START_REF] Valoroso | A damage-mechanics-based approach for modelling decohesion in adhesively bonded assemblies[END_REF])

The critical damage-driving force Y * represents a non-decreasing energy threshold whose value is determined by a monotonically increasing positive function F given as :

Y * =      t 0 Ẏ * dt = F(D) if D ∈ [0, 1[ max τ∈[0,T ] Y (τ) if D = 1 (2)
Typical forms of F are that of power laws or exponential functions, and their explicit expressions are constructed in a way to ensure that the energy dissipated in the formation of a new unit traction-free surface equals the critical strain energy release rate G c . In particular, the two-parameter version of the exponential traction-separation relationship, see also Figure 1, is obtained using :

F(D) = -G c log(1 -D) (3)
We recall that in the above model the computation of the damage state, is completely explicit ; in particular, for damage loading ( Ḋ > 0) at each time τ the damage variable is computed as :

D(τ) = min 1, max (σ≤ τ) {F -1 (Y * (σ))} (4)
see also [START_REF] Valoroso | A damage-mechanics-based approach for modelling decohesion in adhesively bonded assemblies[END_REF] for a detailed account.

In the following the vector collecting the two material parameters to be identified will be denoted as :

x = k G c
(5)

Interface element formulation

A consistent variational framework for problems including in-plane rotational degrees of freedom has been presented by Hughes and Brezzi in [START_REF] Hughes | On drilling degrees of freedom[END_REF]. In particular, in their formulation the stress tensor is not a priori assumed to be symmetric, whereby the role it plays in the theory is complementary to that of the (infinitesimal) rotation field.

The stiffness matrix of the drilling interface element can be obtained starting from the following degenerated Hughes-Brezzi functional [START_REF] Valoroso | A novel interface element with rotational dofs[END_REF] :

Π γ (u, ϕ, τ ) = 1 2 Ω D[[u]] • [[u]]d Ω + Ω τ • (u -ϕ) dΩ - 1 2 γ -1 Ω τ 2 dΩ - Ω f • u dΩ ( 6
)
where u is the displacement vector, ϕ is the drilling rotation, τ is a Lagrange multiplier that plays the same role of the skew-symmetric part of the stress in the continuum membrane element [START_REF] Ibrahimbegovic | A robust quadrilateral membrane finite element with drilling degrees of freedom[END_REF] and γ is a regularizing parameter to be chosen in accordance with the ellipticity condition.

Introducing the standard linear continuous finite element interpolations for displacements and rotations and adding quadratic modes for linking displacements to rotations, the discrete equations obtained by taking the variations of the above functional can be written in matrix form as [START_REF] Valoroso | A novel interface element with rotational dofs[END_REF] : T -

   K e u g e g e,
Ω e γ    u I τ I = f 0 (7)
The variable τ is in general assumed element-wise discontinuous, so that it can be eliminated using static condensation; hence, the stiffness matrix to be assembled at the global level is :

K e = K e u + γ Ω e (g e ⊗ g e ) (8)
The individual terms of the stiffness matrix are evaluated using a three-point quadrature scheme. as shown in [START_REF] Valoroso | A novel interface element with rotational dofs[END_REF], either Newton-Cotes or Gauss quadrature yield the same identical result.

The DCB test and the identification problem

The Double Cantilever Beam is the standard test for obtaining the mode-I fracture toughness G c of adhesives. Different procedures exist for performing the experiment and for data reduction ; in particular, corrected beam theories and compliance calibration methods are the most popular ones. However, apart from requiring the knowledge of the crack length evolution during the test, the mentioned data reduction schemes can exhibit possible incoherencies with the numerical model to be used for computations. For this reason, in the author's opinion the use of a Finite Element model to compute the material parameters should be preferred. The geometry of the test considered is shown in Figure 2. The response of the DCB during the debonding tests is simulated using the interface model discussed in Section 4. This is taken as the constitutive law for interface elements, which have been implemented as a part of a customized version of the FE code FEAP [START_REF] Taylor | FEAP, User and Programmer Manuals[END_REF]. In the numerical simulations the right-end of the structure (intact part) is free whilst on the left side two supports are prescribed and an increasing vertical displacements is imposed at the end of the upper arm. The experimental raw load-deflection curve is shown in Figure 3 in terms of reaction force P versus the relative displacement δ. The present study aims to effective identification of the material parameters of the interface model presented in section 2 via inverse analysis using as data set kinematic full-field data obtained from a DIC procedure [START_REF] Hild | Digital image correlation : from measurement to identification of elastic propertiesa review[END_REF][START_REF] Besnard | Finite-element displacement fields analysis from digital images : application to Portevin-Le Chatelier bands[END_REF]. In a deterministic framework, the optimal material parameters x can be obtained as the solution of a nonlinear programming problem where a suitable cost function is minimized. Such a function can be conveniently defined to represent, at all considered instants, the distance between measured quantities, i.e. displacements inside the monitored sub-domain, and those computed via the mathematical model as a function of the unknown parameters x.

The kinematic data set of the identification process (the measured kinematic quantities) is constituted by a time-space sampling of suitable displacements. In particular, it consists of a collection of n u -dimensional vectors u meas i each referring to a particular instant t i , i = 1 , ... , n t , of the loading process. The values of the same displacements that are instead computed via the FE model will be denoted as

u comp i
. With this notation in hand the identification problem is set as follows :

x = arg min x∈S ω(x) (9)
for the cost function :

ω(x) = n t ∑ i=1 u meas i -u comp i (x) 2 W -1 u ( 10
)
S being is the constraint set, that in the present case is defined as the set of x that are component-wise positive :

S = x ∈ ℜ | x j > 0, j = 1, 2 (11)
In the previous equation W -1 u denotes a weighting operator that scales the terms entering the definition of the residual in order to render them all of the same magnitude.

The minimization of the objective function ( 10) is performed through either using the Simplex method and/or a gradient-based descent algorithm available in the Matlab package [START_REF]Matlab 7.8 -Optimization Toolbox 4[END_REF], for which an interface has been developed to communicate with the FE code.

Results and discussion

Symmetric double cantilever beam (DCB) specimens were made using two 8 mm thick, 200.0 mm wide and 20.0 mm deep Al 2024 T351 aluminum coupons, bonded with a layer of acrylic adhesive and separed by an initial crack of length a = 50 mm that is used as the starting defect. The tested DCB assembly was connected to the testing device using two 30mm-long by 10mm-thick fixtures after having coated the surface with sprayed black and white paints to create the random texture necessary for measuring kinematic fields measurement using DIC.

The experiment was monitored within a ROI of about 70 × 10mm 2 using a digital camera. The displacement field has been estimated using a correlation code relying upon a FE discretization that uses bilinear shape functions [START_REF] Besnard | Finite-element displacement fields analysis from digital images : application to Portevin-Le Chatelier bands[END_REF] ; the element size in the FE mesh used for image correlation (typically of the order of some micrometers) being too small to be used for identification purposes, a suitable interpolation has been adopted in order to map the measured diplacements onto the FE mesh to be used for identification.

Figure 4 shows at two different instants the kinematic measurements provided by DIC mapped onto a FE mesh that uses a discretization of 1 × 1mm 2 continuum elements. In the following are analyzed the results of identification exercises based on such data. Computations are first carried out using a FE mesh of 4400 Enhanced Assumed Strain quadrilaterals and 150 interface elements. The result of identification is the vector :

x 1 = k G c = 2583.875 N/mm 3 0.198390N/mm (12)
Figure 5 shows the surface plot of the cost function of this series of identification exercises for a quite wide range of material parameters. The shape of the surface confirms what has been already noted in [START_REF] Valoroso | Characterization of a cohesive-zone model describing damage and decohesion at bonded interfaces. Sensitivity analysis and mode-I parameter identification[END_REF], i.e. that a good estimation of the fracture energy G c is likely to be obtained in almost all cases, while a good estimate of the interface stiffness is more difficult. Figure 6 shows the surface plot of the cost function of this series of identification exercises. Though the shape of the surface remains the same as before, the convergence path takes place almost for a constant value of the interface stiffness. However, the good performances of plane elements with drilling dofs for this problem are evident since the solution ( 13) is obtained for a much coarser mesh. In the third series of computations we tried to enrich the information content of the available measures by using an information about the rotation field tant can be obtained from the available measurements of the displacement field. In particular, the drilling rotations computed using the interpolated nodal displacements on the adopted FE mesh and a finite different scheme. The result is shown in Figure 7 in terms of additional data (the infinitesimal rotation field) and in Figure 8 in terms of cost function and convergence path of the minimization algorithms. In this case it is clearly seen that during the identification process the sequence generated by the minimization algorithms lies mostly along the steepest descent path. However, the solution vector The reason of this apparent contradiction is probably the fact that in this case the quality of used data (i.e. their information content) is not of the same for displacements and rotations, whereby adding rotations into the data set has the effect of producing noise rather than improving precision. Improvements in the estimation of material parameters can however be achieved by properly using the results of the sensitivity analysis. The sensitivity information can indeed be very helpful in the development of inverse analysis procedures since it can be taken as the basis of a selection criterion for the choice of the measurable quantities with the highest information content to be included in the cost function.

Figure 2 -

 2 Figure 2 -DCB test specimen.

Figure 3 -

 3 Figure 3 -DCB test. Experimental load-deflection curve.

Figure 4 -

 4 Figure 4 -DCB test. Interpolated displacement field computed via DIC at two different instants.

Figure 5 -

 5 Figure 5 -Cost function surface. Classical continuum and interface elements

Figure 6 -

 6 Figure 6 -Cost function surface. Continuum and interface elements with drilling dofs

  x,y); U 0 =2.1046 mm y [mm]

Figure 7 -

 7 Figure 7 -DCB test. Rotation field computed using interpolated displacements obtained from mapping of DIC displacements.

  not very reliable since since the value of G c is too close to the physical limit value 0.00.

Figure 8 -

 8 Figure 8 -Cost function surface. Continuum and interface elements with drilling dofs. Rotations included in data set

Acknowledgements

This work has been carried out as part of the project Analysis of structural adhesive joints. Modelling and experimental testing via DIC within the MiSE-ICE-CRUI 2008 research program. The financial support of the Italian Ministry of Economic Development (MISE) and the Italian Institute for Foreign Trade (ICE) are gratefully acknowledged.