Paul Feliot

Julien Bect

Emmanuel Vazquez

PGMO -COPI'14

Keywords: Optimization, Kriging, Gaussian Process, Subset Simulation, Sequential Monte Carlo Mathematics Subject Classification (2000) 90C56, 90C59, 62K99, 62L05, 62P30

This talk addresses the problem of derivative-free global optimization of a real-valued function under multiple inequality constraints. Both the objective function and the constraint functions are assumed to be smooth, nonlinear, expensive-to-evaluate black-box functions. As a consequence, the number of evaluations that can be used to carry out the optimization is very limited. We focus in this work on the case of strongly constrained problems, where finding a feasible design, using such a limited budget of simulations, is a challenge in itself. The method that we propose to overcome this difficulty has its roots in the recent literature on Gaussian process-based methods for reliability analysis-in particular, the Bayesian Subset Simulation (BSS) algorithm of Li, Bect and Vazquez-and multi-objective optimization. More specifically, we consider a decreasing sequence of nested subsets of the design space, which is defined and explored sequentially using a combination of Sequential Monte Carlo (SMC) techniques and sequential Bayesian design of experiments. The proposed method obtains promising result on challenging test cases from the literature.

This talk addresses the problem of derivative-free global optimization of a real-valued function under multiple inequality constraints:

Minimize f (x), Subject to x ∈ X and c(x) ≤ 0,
where f is the function to be minimized, X ⊂ R d is the design space and c = (c 1 , . . . , c q ) is the vector of constraint functions, c j : X → R, 1 ≤ j ≤ q.

Both the objective function f and the constraint functions c j are assumed to be smooth, nonlinear, expensive-to-evaluate black-box functions. More specifically, it is assumed that the values of f (x)

and c(x), for a given x ∈ X, are provided simultaneously by a single call to some time-consuming computer program-a setup that typically applies to industrial design problems, where numerical simulations are used to mimic the actual physical behavior of the system to be designed. Such simulations may for instance require fluid dynamic, heat transfer or mechanical deformation computations and can take from a few hours to several days to compute.

The number of runs that can be afforded to carry out the optimization is therefore very limited. We focus in this work on the case of strongly constrained problems, where finding a feasible design, using such a limited budget of simulations, is a challenge in itself.

Proposed method

Global optimization methods have been investigated intensively for the last decades. When expensiveto-evaluate functions are involved, cheap-to-evaluate approximations of the objective and constraints functions-often referred to as surrogate models or meta-models-are classically relied upon. We adopt here a Bayesian approach, which provides not-only natural surrogate models for the expensive-to-evaluate functions that we have to deal with, but also an elegant framework to help design efficient optimization algorithms. More precisely, the objective function f and the constraint functions c j are modeled as (independent) Gaussian processes, following a now classical approach that has been made popular by We focus in this talk on the case of strongly constrained problems, where the volume of the feasible space is small compared to the size of the design space. As a consequence, locating even a single feasible point becomes difficult and most existing Bayesian optimization methods, which require at least one feasible point to begin with, fail to be applicable. The method that we propose to overcome this difficulty has its roots in the literature on Bayesian sequential design of experiments for reliability analysis (see, e.g., [START_REF] Bect | Sequential design of computer experiments for the estimation of a probability of failure[END_REF][START_REF] Li | Bayesian Subset Simulation: a kriging-based subset simulation algorithm for the estimation of small probabilities of failure[END_REF]-estimating a failure region or a feasible set are very similar problems-and multi-objective optimization (see, e.g., [START_REF] Emmerich | Single-and multi-objective evolutionary optimization assisted by Gaussian random field metamodels[END_REF][START_REF] Wagner | On expected-improvement criteria for modelbased multi-objective optimization[END_REF]Picheny, 2014b). More specifically, we consider a decreasing sequence of nested subsets of the design space, which is defined and explored sequentially using a combination of Sequential Monte Carlo (SMC) techniques and sequential Bayesian design of experiments, in the spirit of [START_REF] Li | Bayesian Subset Simulation: a kriging-based subset simulation algorithm for the estimation of small probabilities of failure[END_REF]; [START_REF] Benassi | Bayesian optimization using sequential Monte Carlo[END_REF]; [START_REF] Li | Sequential Design of Experiments to Estimate a Probability of Failure[END_REF]; [START_REF] Benassi | Nouvel algorithme d'optimisation bayésien utilisant une approche Monte-Carlo séquentielle[END_REF].

Results and future work

We are able to report good results on challenging test cases from the literature. Future work will include the extension of our method to multi-objective problems and the optimization of various aspects of our algorithm (Sequential Monte Carlo algorithm, sampling criterion...). Simulation failures should also be taken into account as they are inherent to complex industrial simulation codes. The performance of the method will be evaluated on a real-life industrial problem provided by Safran, with the contribution of Cenaero (optimization of the performances of a turbo-machine fan blade under aerodynamic, mechanic and acoustic constraints).

  [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF] for unconstrained optimization problems. Subsequent developments for constrained optimization problems have been proposed by, among others,[START_REF] Schonlau | Global versus local search in constrained optimization of computer models[END_REF];[START_REF] Sasena | Exploration of metamodeling sampling criteria for constrained global optimization[END_REF];[START_REF] Gramacy | Optimization under unknown constraints[END_REF];[START_REF] Parr | Infill sampling criteria for surrogate-based optimization with constraint handling[END_REF];Picheny (2014a).

Acknowledgements This research work has been carried out in the frame of the Technological Research Institute Sys-temX, and therefore granted with public funds within the scope of the French Program "Investissements d'Avenir".


R M Gelet

B Loret

N Khalili

The significance of local thermal non-equilibrium in simulations of enhanced geothermal recovery

. Hydraulic diffusion, thermal diffusion, forced convection and deformation are considered in a two-phase framework, the solid phase being made by impermeable solid blocks separated by saturated fractures. Each of the two phases is endowed with its own temperature. The resulting system of equations is used to address a generic HDR reservoir subjected to temperature and pressure gradients. A change of time profile of the outlet fluid temperature is observed as the fracture spacing increases, switching from a single-step pattern to a double-step pattern, a feature which is viewed as characteristic of established LTNE. A dimensionless number is proposed to delineate between local thermal equilibrium (LTE) and non-equilibrium. This number embodies local physical properties of the mixture, elements of the geometry of the reservoir and the production flow rate. All the above properties being fixed, the resulting fracture spacing threshold between LTNE and LTE is found to decrease with increasing porosity. The thermally induced effective stress is tensile near the injection well, illustrating the thermal contraction of the rock, while the pressure contribution of the fracture fluid is negligible during the late period.

INTRODUCTION

The production of geothermal energy is achieved by cold water injection, in fractured igneous rocks (originally with low matrix permeability), at sites where the vertical temperature gradient is favorable. Of crucial importance to the economic viability of these HDR reservoirs is the knowledge of thermal output evolution, induced thermal stress and fluid loss, at various time scales of the circulation tests [START_REF] Armstead | Heat mining[END_REF]. A key factor influencing geothermal energy recovery is the difference in characteristic times between diffusion in the solid phase and convection in the fluid phase. The thermally induced effective stress which results from these two contributions may lead to permeability change and fluid loss. Also, thermal shrinkage and pressure drop, across the body of the reservoir, occur at various time scales and the understanding of their interaction requires the simulation of a continuum mixture in LTNE [START_REF] Aifantis | Further comments on the problem of heat extraction from hot dry rocks[END_REF].

Coupled thermo-hydro-mechanical analyses for geothermal systems are scarce [START_REF] Hayashi | Numerical models of HDR geothermal reservoirs -a review of current thinking and progress[END_REF][START_REF] O'sullivan | State of the art of geothermal reservoir simulation[END_REF]) and focus mainly on partially coupled systems in an effort to implement a precise fracture network through a system of discrete discontinuities [START_REF] Duteaux | Insight from modelling discrete fractures using GEOCRACK[END_REF][START_REF] Bruel | Impact of induced thermal stresses during circulation tests in an engineered fractured geothermal reservoir: Example of the Soultz-Sous-Forêts european hot fractured rock geothermal project, rhine graben, france[END_REF], to couple free and forced convection [START_REF] Bataillé | Modeling the coupling between free and forced convection in a vertical permeable slot: implications for the heat production of an enhanced geothermal system[END_REF], or to characterize joint closure with a stress dependent law [START_REF] Kohl | FRACture" -a simulation code for forced fluid flow and transport in fractured, porous rock[END_REF][START_REF] Bower | A numerical model for thermo-hydro-mechanical coupling in fractured rock[END_REF] or stress-dissolution/precipitation effects [START_REF] Taron | Numerical simulation of thermal-hydrologic-mechanical-chemical processes in deformable, fractured porous media[END_REF]. However, before [START_REF] Gelet | Thermo-hydro-mechanical study of deformable porous media with double porosity in local thermal nonequilibrium[END_REF] no numerical analysis was proposed to simulate reservoir behavior with a continuum mixture in LTNE.

Extending previous studies, [START_REF] Loret | A three-phase model for unsaturated soils[END_REF] and [START_REF] Khalili | An elasto-plastic model for nonisothermal analysis of flow and deformation in unsaturated porous media formulation[END_REF] proposed comprehensive fully coupled thermo-hydro-mechanical formulations based on the mixture theory for unsaturated porous media. Both contributions focus on the importance of an appropriate definition of the effective stress. In the present contribution, the work of [START_REF] Khalili | An elasto-plastic model for nonisothermal analysis of flow and deformation in unsaturated porous media formulation[END_REF] is extended to describe thermal recovery from a geothermal reservoir viewed as a deformable fractured medium in LTNE [START_REF] Gelet | Thermal recovery from a fractured medium in local thermal non-equilibrium[END_REF].

THERMOPOROELASTIC MODEL IN LTNE

Hydraulic diffusion, thermal diffusion, forced convection and deformation are considered in a twophase framework, the solid phase being made by impermeable solid blocks separated by saturated fractures. Local thermal non-equilibrium between the phases holds at all points of the simulated zone. Each phase k is endowed with its own temperature T k and its own material properties. LTNE emanates from the contrast between the rapid convection of heat by the moving fluid in the fractures and the slow diffusion of heat through the solid matrix. Indeed, the characteristic time associated with convection is several orders of magnitude smaller than the characteristic times of diffusion in both solid skeleton and fluid.

Governing equations

Within a single porosity conceptual framework [START_REF] Gelet | Thermal recovery from a fractured medium in local thermal non-equilibrium[END_REF], differential equations describing the problem at hand may be expressed as,

µ DS ∇ 2 u + (λ DS + µ DS ) ∇(div u) -ξ f ∇p f - c T c ∇T s + ρ g = 0 , (1)
-div k f µ f (∇p f -ρ f g) + a f f ∂p f ∂t +ξ f div ∂u ∂t + a f Ts ∂T s ∂t + a f T f ∂T f ∂t = 0 , (2)
-div (n s Λ s ∇T s ) + a TsTs ∂T s ∂t +T s a f Ts ∂p f ∂t + T s c T c div ∂u ∂t +κ sf (T s -T f ) = 0 , (3) -div(n f Λ f ∇T f )+a T f T f ∂T f ∂t +T f a f T f ∂p f ∂t + ρ f J f • ∇H f +κ sf (T f -T s ) = 0 , (4)
in which

ξ f = 1 -c s /c a f f = n f c f H + (ξ f -n f )c s a f Ts = (n f -ξ f )c T a f T f = -n f c f T a TsTs = n s ρ s C (v) s a T f T f = n f ρ f C (p) f ∇H f = (1 -T f c f T ) ∇p f ρ f + C (p) f ∇T f . ( 5
)
where u is the displacement of the solid, p f is the pressure of the fracture fluid, T s is the temperature of the solid and T f that of the fracture fluid. Throughout, the indices s and f refer to the solid and to the fluid, respectively. Equation (1) uses the Lamé constants of the drained solid λ DS and µ DS , the body force vector g, the overall density ρ, the compressibility of the mixture c, the compressibility of the solid grains c s , and the volumetric thermal expansion coefficient of the solid phase c T . Equation (2) requires the macroscopic porosity n f , the macroscopic intrinsic permeability k f , and the dynamic viscosity µ f . The hydraulic compressibility c f H and the thermal compressibility c f T are defined from the intrinsic density ρ f of the fluid. Equation (3) involves the solid heat capacity at constant volume and fluid pressure C

(v)
s , the solid density ρ s and porosity n s , and the solid thermal conductivity Λ s . The specific inter-phase heat transfer parameter parameter κ sf is discussed in Section 2.2. Finally, equation (4) involves the fluid heat capacity at constant fluid pressure C (p) f , the apparent volume flux J f (Darcy's law), and the fluid thermal conductivity Λ f . Equations (1) to (4) are direct consequences of the field equations (balance of momentum for the mixture, balance of mass for the fluid, balance of energy for the solid and balance of energy for the fluid). Additional information on the parameters of the model can be found in previous works [START_REF] Gelet | Thermo-hydro-mechanical study of deformable porous media with double porosity in local thermal nonequilibrium[END_REF][START_REF] Gelet | Thermal recovery from a fractured medium in local thermal non-equilibrium[END_REF].

The approach is amendable to a double porosity framework: it can be extended to zoom on the solid blocks and include their innate pore pressure, temperature and the fluid and energy exchanges with the fractures [START_REF] Gelet | A thermo-hydromechanical model in local thermal non-equilibrium for fractured HDR reservoirs with double porosity[END_REF].

The weak form of the governing equations, the discretization and time-integration procedures to solve the coupled equations through a finite element method are detailed in [START_REF] Gelet | Thermo-hydro-mechanical study of deformable porous media with double porosity in local thermal nonequilibrium[END_REF]. The primary variables are the displacement vector, the pressure of the fluid, the temperature of the solid and the temperature of the fluid. The resulting system of equations is used to address a generic HDR reservoir subjected to temperature and pressure gradients.

The specific inter-phase heat transfer parameter κ sf

The coefficient of specific inter-phase heat transfer κ sf [W/m 3 .K] controls the rate at which the two phase system (solid-fluid) reaches thermal equilibrium. The higher this coefficient, the faster thermal equilibrium is reached. This coefficient can be defined by idealizing the solid blocks to squares of length B bordered by fractures of width 2b (Fig. 1).

In geothermal applications, the fracture width is much smaller than the fracture spacing, namely 2 b ≪ B. Hence, the specific coefficient of heat transfer between the solid and the fluid phase can be linearly related to B -2 [START_REF] Gelet | Thermal recovery from a fractured medium in local thermal non-equilibrium[END_REF],

κ sf = 8 B 2 n s Λ s × 2 Λ f n s Λ s + 2 Λ f . (6)
The thermal conductivities of the two phases are involved. For a geothermal reservoir, typical values entail n s ≈ 1, Λ s = 2.71 W/m.K and Λ f = 0.6 W/m.K, and the two phases are seen to contribute to the overall conductivity.

Since κ sf depends on the solid-fluid specific surface, the experimental measurement of this coefficient proves to be difficult. One way around is to calibrate the specific inter-phase heat transfer coefficient, with help of data provided from a specific HDR reservoir as proposed in [START_REF] Gelet | A thermo-hydromechanical model in local thermal non-equilibrium for fractured HDR reservoirs with double porosity[END_REF][START_REF] Gelet | Thermal recovery from a fractured medium in local thermal non-equilibrium[END_REF]. Instead, here, a sensitivity analysis is performed to determine the threshold value that separates local thermal equilibrium from the local thermal non-equilibrium.

HDR RESERVOIR ANALYSIS

The fully coupled thermo-hydro-mechanical model is used in the subsequent section to investigate the thermal drawdown of a generic HDR reservoir in LTNE. Information on the geometry and the boundary conditions are provided first.

Geometry

An artificially fractured reservoir with horizontal injection and production wells is considered (Fig. 1). The simulations assume a plane strain analysis, in the xz plane, and symmetry with respect to z-axis. The finite element mesh includes half of the reservoir and a portion of the surrounding low permeability rock formation. The HDR fractured reservoir is idealized by a single porosity saturated medium.

z x 0 0 Solid Fluid X R X F X W B k f 2b Z R
The injection and the production wells are located at the bottom and at the top of the reservoir, respectively. The length of the wells is a key parameter that governs the flow path and the heat transfer between the rock and the fluid. Here, the wells are assumed to penetrate horizontally into the entire reservoir, that is X W = X R , and the velocity of the fluid v f is vertical and uniform between the two wells. The main purpose of this paper is to focus on the influence of LTNE and alternative setups (flow paths) are described elsewhere [START_REF] Gelet | Thermal recovery from a fractured medium in local thermal non-equilibrium[END_REF]).

Boundary conditions

The thermal, hydraulic and mechanical boundary conditions are sketched in Fig. 2.

The injection temperature T inj is applied to all the phases along the injection well. The outlet temperature T f,out along the producing well is an output of the simulations. The remaining horizontal and vertical boundaries are thermally insulated from the surroundings.

The injection and production pressures, p f,inj and p f,out respectively, are both specified along the injection and production wells. The remaining boundaries are hydraulically impermeable.

The vertical stress σ v remains constant along the top boundary z = Z R , and so does the horizontal stress σ H along the lateral boundary x = X F . The displacements on the other boundaries are specified, u x (x = 0, z) = 0 and u z (x, z = 0) = 0.

At the contact between the reservoir and the surrounding rock x = X R , the reservoir is assumed to be thermally insulated from the rock formation. Alternative setups are studied in [START_REF] Gelet | Thermal recovery from a fractured medium in local thermal non-equilibrium[END_REF]. (1) or ( 2)

T s = T f = T inj X F X W z x Z R T s = T f = T 0 X R - 0 0 Impermeable boundary - Impermeable boundary X F z x Z R p p = p f = p inj p p = p f = p out Formation X W Impermeable boundary Impermeable boundary X R - 0 0 X F z x Z R No vertical displacement No horizontal displacement Formation X R Thermal BC Hydraulic BC Mechanical BC σ v σ H

NUMERICAL RESULTS

The resulting system of equations is used to address a generic HDR reservoir subjected to temperature and pressure gradients with a focus on the LTNE assumption. First, thermal drawdown curves (time profiles of the fluid temperature at the producing well) are observed to show a single-step pattern for small fracture spacings B and double-step patterns for large fracture spacings (Sect. 4.1). The threshold values between LTNE and LTE are sought in terms of the fracture spacing B T for an average steady state velocity v ∞ and several fluid porosities in Sect. 4.2. Next, a dimensionless parameter is proposed in Sect. 4.3 to delineate between local thermal equilibrium (LTE) and non-equilibrium. Finally, the reservoir response in terms of overall effective stress is presented in Sect. 4.4. Throughout, temperatures are presented in a dimensionless format. The dimensionless fluid temperature T D built from the outlet fluid temperature T f,out , the initial temperature T 0 and the injection temperature T inj ,

T D = T 0 -T f,out T 0 -T inj , (7)
ranges between 0 (T f,out = T 0 ) at early time and 1 (T f,out = T inj ) ultimately.

The double step pattern

A change of time profile is observed as the fracture spacing B increases (Fig. 3, center). Time profiles for small B are clearly dominated by convection and display a characteristic single-step pattern. Conversely, the time profiles for large B display three stages and typical double-step patterns.

(1) In the first stage, the strong increase of the dimensionless temperature corresponds to the abrupt propagation of the injection temperature dominated by convection. This effect is mostly attributed to the difference in characteristic times between forced convection in the fracture network and diffusion of heat in the solid matrix. Accordingly for each fracture porosity n f , the magnitude of the dimensionless temperature T D at the end of this first stage is increasing as a function of B. Thus, the smaller B (the larger the specific solid-fluid heat transfer κ sf ), the higher the 'instantaneous damping' of the cooled fluid front by the solid phase. On the other hand, the duration of the first stage is independent of B, and it is about equal to the characteristic time for a particle of fluid to flow the length Z R at the steady velocity v ∞ ,

t stage 1 c = Z R v ∞ ≈ 2 × τ S1 = 13 days . (8)
(2) During the second stage, the dimensionless temperature remains constant and displays a plateau. Heat transfer between the solid and the fluid phases takes place partly at constant fluid temperature.

(3) The third stage is characterized by a second strong increase of the temperature T D . The time required for the outlet temperature T f,out to reach 95% of the injection temperature T inj is higher for a model in LTNE than in LTE. In other words, the response of the mixture is delayed by the transfer of heat from the solid phase to the fluid phase. This type of behavior is characteristic of the response of phases in nonequilibrium: as an example, the consolidation time of a dual porosity medium is delayed compared with a single porosity medium [START_REF] Khalili | Consolidation of fissured clays[END_REF]).

Influence of the porosity n f

The influence of the fracture porosity, in the range of 0.001 < n f < 0.01, is mainly visible during the first and the second stage of the double-step pattern (Fig. 3). ( 1) For each fracture spacing B, the magnitude of T D at the end of this first stage is increasing as a function of the fluid volume ratio n f , i.e. the larger the volume of cooled fluid, the smaller the temperature outlet. Since n s ≈ 1, the rate of heat transfer can be considered as constant, eqn (6), as opposed to the overall amount of heat supply in the system. In addition, the duration of the first stage is independent of n f in agreement with eqn (8). 7), as a function of time for three porosities n f and three fracture spacings B. All results are for Z R = 230 m, v ∞ = 2×10 -4 m/s and thermal properties from [START_REF] Gelet | Thermal recovery from a fractured medium in local thermal non-equilibrium[END_REF]. In LTNE, time profiles display a double-step pattern, whereas in LTE, time profiles display a continuous pattern. The thresholds between LTE and LTNE are associated with the fracture spacings B T = 6 m (left), B T = 2.5 m (middle) and B T = 2 m (right). The corresponding dimensionless values of η D are equal respectively to 11, 13 and 10. Therefore, LTE is associated with η D larger than, say 13, while values of η D smaller than 13 require a LTNE analysis. (2) Time profiles in Fig. 3 indicate that increasing the porosity of the fracture network n f reduces the time length of the second stage and, hence, speeds up thermal depletion. An increase of fluid porosity modifies the phase and the heat partition in the system and hence reduces the overall amount of heat to be transferred by the solid to the fluid. Since the rate of heat transfer is almost independent of the fluid porosity, see eqn (6) in which n s ≈ 1, the heat transfer period requires less time.

In terms of fracture spacing, the threshold between LTE and LTNE decreases with increasing fluid porosity, B T = 6 m (n f = 0.001), B T = 2.5 m (n f = 0.005) and B T = 2 m (n f = 0.01), due to the increasing amount of cooled fluid. Those values are obtained by trial and error and correspond to T s = T f , at the reservoir outlet, over the entire circulation test. For smaller fracture spacings B < B T , the local thermal non-equilibrium responses are indistinguishable from the responses in equilibrium.

A dimensionless LTNE parameter

In the context of geothermal energy recovery, reservoir performances can be expressed in terms of a dimensionless LTNE parameter [START_REF] Gringarten | Theory of heat extraction from fractured hot dry rock[END_REF]. The dimensionless LTNE parameter η D serves to delineate LTE and LTNE and reflects the overwhelming importance of the characteristic lengths of the fracture network, namely the average aperture 2 b and the average spacing B.

The general form of η D [START_REF] Gelet | Thermal recovery from a fractured medium in local thermal non-equilibrium[END_REF]) expresses in terms of the thermal properties of the porous medium (Λ s , Λ f and C (p) f ), the fluid porosity n f , the diffusion-convection length Z R , the fracture spacing B, and the steady state fluid velocity v ∞ . For n s Λ s ≫ n f Λ f and for n s ≈ 1, the dimensionless LTNE parameter writes

η D = 1 n f 16 Λ f Λ s + 2 Λ f Z R α T,f v ∞ B 2 . (9) in which α T,f = n f Λ f /ρ f C (p)
f is the thermal diffusivity of the fluid and B depends on the parameters n f and 2b, B = 4b/n f , for idealized solid square blocks of length B borded by fractures of width 2b (Sect. 2.2 and Fig. 1).

In the previous section, three fracture spacing thresholds B T have been identified, one for each porosity. The corresponding dimensionless values of η D are equal respectively to 11 (n f = 0.001), 13 (n f = 0.005) and 10 (n f = 0.01). Consequently, within the range of interest 0.001 ≤ n f ≤ 0.01, dimensionless values η D ≥ 13 are likely to represent LTE and η D < 13 represents LTNE.

THM coupled behavior

The LTNE model is used to investigate the coupled thermo-hydro-mechanical behavior of the Fenton Hill HDR reservoir (Fig. 4).

The coupled behavior of porous media in LTNE is governed by the difference in characteristic times between the thermal depletion of the fluid phase and that of the solid phase. During the first weeks, the thermal depletion of the fluid has already started to propagate inside the reservoir, whereas the thermal depletion of the solid is confined to the immediate vicinity of the injection well. However, as heat transfer between the two constituents gets completed, the mixture tends to thermal equilibrium. On the other hand, the pressure field reaches steady state within a few hours due to the high hydraulic diffusivity of the fracture network. The large difference between the characteristic times of the fracture fluid response (early/intermediate period) and of the solid temperature response (late period) allows segregation of the respective contributions to the effective stress due to hydraulic and thermal effects. As expected, the latter is larger than the former especially in the horizontal and out of plane directions. The effective stresses are governed by the thermal depletion of the solid phase during the late period. Cooling of the lower part of the reservoir induces characteristic contractive strains and tensile stresses.

CONCLUSIONS

Diffusion and forced convection mechanisms between two phases have been studied in the context of heat extraction from a fractured hot dry rock (HDR) reservoir. A model describing the behavior of poroelastic fractured media has been extended to account for local thermal non-equilibrium (LTNE).

A parametric analysis has been carried out to study the influences of the solid-fluid heat transfer coefficient and of the fluid porosity. As a typical feature of LTNE, the temperature outlet profile displays a double-step pattern, which is interpreted as a characteristic of established LTNE. An attempt to define a dimensionless threshold above which a LTNE analysis is required is exposed (Sect. 4.3).

The model is used to describe the thermo-hydromechanical behavior of a long term circulation test on the Fenton Hill HDR reservoir. It has been found that the thermally induced effective stresses are tensile near the injection well due to the thermal contraction of the solid, suggesting a possible increase of fracture aperture. This work highlights how LTNE versus LTE can strongly impact the mechanical response of the reservoir and possibly affect the evolution of the fracture network. The LTNE model leads to a characteristic double-step reservoir response.

Figure 1 :

 1 Figure 1: Representation of a generic HDR reservoir (not at scale). The permeability k f of the reservoir depends on the average fracture spacing B and on the average fracture aperture 2 b.

Figure 2 :

 2 Figure 2: Thermal, hydraulic and mechanical boundary conditions.

Figure 3 :

 3 Figure 3: Dimensionless temperature outlet T D , eq. (7), as a function of time for three porosities n f and three fracture spacings B. All results are for Z R = 230 m, v ∞ = 2×10 -4 m/s and thermal properties from[START_REF] Gelet | Thermal recovery from a fractured medium in local thermal non-equilibrium[END_REF]. In LTNE, time profiles display a double-step pattern, whereas in LTE, time profiles display a continuous pattern. The thresholds between LTE and LTNE are associated with the fracture spacings B T = 6 m (left), B T = 2.5 m (middle) and B T = 2 m (right). The corresponding dimensionless values of η D are equal respectively to 11, 13 and 10. Therefore, LTE is associated with η D larger than, say 13, while values of η D smaller than 13 require a LTNE analysis.

Figure 4 :

 4 Figure 4: Late period profiles, along the z-axis at x = 60 m, of the induced change in vertical effective stress (left), lateral effective stress (center) and out of plane effective stress (right), with n f = 0.005, k f = 8.0×10 -15 m 2 , κ sf = 33.0 mW/m 3 .K. Tensile stresses are counted positive. The contraction of the solid phase induces the effective lateral stress to be tensile near the injection well (z < 120 m) and compressive near the production well (z > 120 m).


Keywords: marine protected area, socio-economic impact, fishing unit profitability, fishery household income, methodological proposals, village monographs, sample-based survey, multivariate analysis, Chisquare tests, Lorenz curves

The assessment of socio-economic effects of marine protected areas on fishing activity points out some methodological difficulties. As a methodology we propose a selection of two components: the fishing unit profitability and the fishery household income distribution and a comparison between a marine protected area versus an "unprotected zone". The methodology has been applied on a case study: the Chumphon National Park versus the Chumphon Province (Gulf of Thailand). Three phases have been implemented: a bibliographical analysis, the carrying out of village monographs, the carrying out of a sample-based survey. The sample-based survey results are presented as follows : the characterization of categories of the variable "marine protected area" by quantitative variables, a multivariate analysis, the performance of Chi-square tests on the relationship between a "profit per fishing day" variable divided into two categories ("in or around MPA" and "remote from the MPA", the measurement of the concentration of income or operating profit per fishery household and the Lorenz curves.

Methodological difficulties and proposals

The assessment of socio-economic impacts of a marine protected area on fishing activities points out some methodological difficulties, particularly when multifleet-multispecies fisheries are active, explain pro parte a weak research effort on the socioeconomic impact of fishing activities after the implementation of a marine protected area.

To assess this impact, we proposed: a selection of two components of this impact (the fishing unit profitability and the fishery household income distribution ) and a comparison between a marine protected area versus an "unprotected zone"

The institutional context was a research agreement between IRD (Institut de Recherche pour le Développement) and Kasetsart University (Faculty of Economics) under the supervision of Pr Ruangrai Tokrishna. This programme was funded by IRD, ECOST Project and Kasetsart University. The case study was the Chumphon Province and the Mu Ko Chumphon National Park (MKCNP) along the Gulf of Thailand.

A case study : Chumphon National Park versus Chumphon Province

(Gulf of Thailand)

Chumphon Province

-6010 square kilometers, 446000 inhabitants, a population density of 74 inhabitants per square kilometer -8 districts : Mueang Chumphon, Lang Suan, Thung Takao, Sawi, Pathio, Tha Sae, Lamae, Phato -2880 fishery households and 9580 fishery household members

Mu Ko Chumphon National Park (MKCNP)

-established in 1999 under the supervision of the Department of National Parks : 317 square kilometers of which 265 of marine area; 70 kilometers of coastline -four main ecological systems : coral reef, seagrass bed, mangrove forest, limestone forest and tropical rainforest, food-plain mire, mud beach -MKCNP is spreading over a part of fourth Chumphon province districts : Mueang Chumphon, Lang Suan, Thung Takao, Sawi, Pathio, Tha Sae, Lamae, Phato -500 fishery households and 1700 fishery household members (inside the Park and in communes just adjacent to the Park) -a zoning : a strict nature reserve and a general use zone in which the fishing from the residents is tolerated

The method

Bibliographical analysis

A first phase devoted to a bibliographical analysis and interviews with Chumphon Province Department of Fisheries (DOF) officers and Direction of National Parks (DNP) officers focused on the features of fisheries activities and on the legislative and regulatory framework.

Village monographs

A second phase devoted to field investigations which took the shape of village monographs and a sample-based survey -The fishers village monographs were focused on fisheries activities (type of fishing units, seasonal activities, location of fishing grounds) and on fishery household occupational structure, monetary costs and fishery profitability, fishery household income distribution by fishing unit and by extra fishing source: 7 fishers village monographs in Mu Ko Chumphon National Park (MKCNP) or adjacent to the MKCNP and 12 fishers village in the rest of Chumphon Province -The sample-based survey was focused on fishery profitability and income distribution : 126 fishery households forming 225 fishing units were surveyed

A sample based survey

A third phase devoted to the processing of village monographs and to the sample-based survey

The sampling strategy

The sampling unit is the fishery household which forms one or several fishing units The data source is a 2006 census from the Department of Fisheries of Chumphon Province (Gulf of Thailand). Two types of stratification :

-a geographical stratification : the communes (tambon) inside or adjacent to the MPA, called "in and around the MPA", versus the remote communes from the MPA called "remote from the MPA" -a stratification by the main "métier" used in and around the MPA: otter board trawler, anchovy purse seiner, anchovy falling netter, squid falling netter. The sampling rate was 15%: 126 fishery household forming 225 fishing units were drawn.

The survey method and the carrying out of the sample survey

The drawing of the sample was as follows: for each stratum one must have a list of fishery household, for each list one applies a systematic random procedure: the first fishery household is drawn at random, then one fishery household is drawn every five fishery household.

The frequency of the survey: every fishery household has been surveyed twice a year relating to the year 2007

The criterion) five clusters that present several determining features (Fig. 1).

 A performance of Chi-square tests to test the relationship between profit per fishing day (PPFD) variable and MPA variable ("in or around the MPA" versus "remote from the MPA"). Two Chisquare tests have been performed ; the "profit per fishing day" (PPFD) variable was divided into two categories (positive and negative PPFD).  The measurement of the concentration of income or operating profit per fishery household by a single indicator derived from Herfindhal index defined as the sum of the squares of the market shares compared to the sum of incomes or operating profits of all households. The indicator derived from Herfindhal index has two advantages : make possible the comparison of the levels of concentration between samples with different size and mark the indicator between 0 and 1.

 The drawing of Lorenz curves

Main results of village monographs

A fishery household is formed by one or several fishing units defined by a métier which associates a fleet, a fishing gear and target species.

The tolerance regarding the small-scale fishery activities of the MPA residents within the borders or just around Mu Ko Chumphon National Park (MKCNP).

The pointing out of the most representative métiers in and just around the selected marine protected area: the squid falling netters (netter/squid falling net/squid) with 283 households, the otter board trawler (trawler/otter board trawl/trashfish) with 139 households, the anchovy falling netter (natter/anchovy falling net/anchovy) with 60 households, the anchovy purse seine (seiner/anchovy purser seine/anchovy) with 12 households. The identification of main fishing grounds for each commune included in the sample, specially for the communes in and just adjacent to the MPA. This identification confirmed that fishery units from these communes fish in or in the vicinity of the MPA, at the opposite of fishery units from the rest of Chumphon Province. It means the implicit acknowledgement of exclusive access rights for the benefit of residents

Main results of the sample based survey

 from the characterization of the categories of the variable "MPA" by quantitative variables  from the multivariate analysis :

The results are not so evident. But we can point out that the proportion of insiders is higher in the high wage rate cluster (class 1; 7 obs.) and the proportion of outsiders is higher in the low profit rate cluster (class 2; 46 obs.). Or the proportion of outsiders facing negative profit rates is higher than the proportion of insiders.

But, at the overall level of the PCA, the "in or around MPA" variable is not significant in the three other classes. This is confirmed by the position of the MPA categories (inside vs outside): although situated on the left-hand side of the horizontal axis where all the returns and profits variables are linked together (thus the units having the highest returns are rather on this left-hand side of the factorial map), the "INSIDE MPA" category remains close to the centre of the map, hence to the average values of the ten variables. The OUTSIDE MPA category is located on the right-hand side of the map (where the individuals having lower returns are) but it is also quite close to the centre of gravity (average values).  from the performance of Chi-square tests on the relationship between a "profit per fishing day" (PPFD) variable divided into two categories (positive and negative PPFD).

-A first test shows that the outsiders having a low profit are twice more important than the insiders : at the 95% of significance we found a significant relationship between the profit per fishing day level and the MPA variable ("in or around the MPA" or "remote from the MPA")

-A second test shows the greater homogeneity and the lower variability of profits for insiders: four fishing units out of five made up the mid-profit category for the insiders against only two thirds as far as the outsiders are concerned.

 from the measurement of the concentration of income or operating profit per fishery household and from the Lorenz curves

The two Lorenz curves show that the concentration of incomes or operating profits is higher for the outsiders. A such concentration refers to a more non-egalitarian distribution.

In calculating the indicator derived from Herfindhal index on incomes per fishery household (n = 78 IN et n = 51 OUT), we have got a value of 1,72% for the insiders versus 5,81% for the outsiders; it means a higher concentration for outsiders. The non-egalitarian structure of outsiders is more evident in considering the operating profits; if we exclude negative operating profits, the value of the index is 3,45% for the insiders and 14,57% for the outsiders.

Fig 1 :

 1 Fig 1: Average profit per fishing day by métier (bath)

Fig2:

  Fig2: Multivariate analysis on the basis of ten variables

Fig 3 .

 3 Fig 3. Lorenz curve of income per fishery household

Table 0 :

 0 Breakdown of the sample of Chumphon Province fishery household survey

	Métier	In and around the	Remote from the	Total
			MPA		MPA
		(MKCNP)	(outsiders)
		(insiders)
		Census	Sample	Census	Sample	Census	Sample
	Trawler/Otter board trawl/Trashfish	139	22	48	7	187	29
	Seiner/Anchovy purse seine/ Anchovy	12	2	11	1	23	3
	Netter/Anchovy falling net/ Anchovy	60	9	71	11	131	20
	Netter/Squid falling net/Squid	283	42	209	32	492	74
	TOTAL	494	75	339	51	833	126

  carrying out of the sample-based survey : -a preliminary inquiry (February-March 2007) : interviews with fishermen leaders at the level of each commune (tambon) to draw up the sample, to explain about the content of the survey with fishery households -126 fishery households covering 225 fishing units have been surveyed: a first inquiry on July-August 2007 and a second inquiry on February-March 2008 3.3.3. The processing of the sample-based survey  The development of an excel file to enter and store the data  The characterisation of the categories of the variable « in and around the MPA » (insiders) versus « remote from the MPA » (outsiders) by a selection of quantitative variables. We have ranked the quantitative or continuous variables by decreasing order of Test-values significant at 5% for both positive and negative statistics (respectively greater and lower than average values)  A principal component analysis and a clustering with a caracterisation of the categories « in and around the MPA » versus « remote from the MPA » by ten following continuous and nominal variables: average landing price for species 1, average landing price for species 2, average landing price for species 3, profit per fishing day, profit rate, wage rate, catch per fishing day, catch per fishing hour, share of wages in variable costs, profit over crew wages. It optimally gave (Ward

Table 1 .

 1

	IN THE MPA
	Distinctive variables	Mean in the category	Overall mean	Standard deviation in the category	Overall standard deviation	Test-Value	Probability
	HIGHER THAN AVERAGE
	Boat size	14,161	12,574	3,478	3,978	4,64	0,000
	Taxes	1560,530	1196,170	1221,870	1118,730	3,78	0,000
	Horse power	190,842	168,449	90,267	95,660	2,72	0,003
	PROFIT	277511,000	239603,000	195106,000	167480,000	2,63	0,004
	value for species 4	41165,100	25443,600	88270,900	70334,800	2,60	0,005
	Share of wages in VC	0,251	0,226	0,102	0,116	2,47	0,007
	FISHERY INCOME	830071,000	710361,000	585998,000	573477,000	2,42	0,008
	fishing hours	2055,440	1773,620	1554,310	1416,050	2,31	0,010
	cost of labour	38231,600	35184,500	15333,200	15775,500	2,24	0,012
	Depreciation cost	32849,300	28410,200	23497,800	23781,900	2,17	0,015
	value for species 3	61643,500	41798,800	125737,000	108376,000	2,13	0,017
	FIXED COSTS	66917,100	56097,500	62362,400	59623,000	2,11	0,018
	Number of crew	5,351	4,806	2,737	3,012	2,10	0,018
	TOTAL COSTS	552559,000	470757,000	444887,000	462504,000	2,05	0,020
	euro	12279,100	10461,300	9886,380	10277,900	2,05	0,020
	Oil	13402,900	10695,000	18678,400	15387,300	2,04	0,020
	Rehabilitation	34,421	27,757	46,254	38,902	1,99	0,023
	Fuel	233344,000	192189,000	233936,000	241155,000	1,98	0,024
	VARIABLE COSTS	485474,000	414524,000	414062,000	430710,000	1,91	0,028
	MANAGEMENT COSTS	167,784	136,270	230,352	193,454	1,89	0,029
	Handling cost	6180,700	5446,940	5299,900	4520,030	1,89	0,030
	Administration	33,341	27,128	46,166	38,737	1,86	0,031
	Enforcement	50,011	40,693	69,249	58,105	1,86	0,031
	Research	50,011	40,693	69,249	58,105	1,86	0,031
	Fishing days	199,123	185,459	79,861	85,712	1,85	0,032
	Wage rate	22410,100	19512,900	15242,900	18461,700	1,82	0,034
	LOWER THAN AVERAGE
	Average landing price for species 6	27,000	34,294	26,098	26,537	-1,98	0,024
	Profit over crew wages	3,572	5,911	2,979	10,894	-2,49	0,006
	Profit rate	0,355	0,392	0,128	0,155	-2,74	0,003

Table 2 :

 2 Profit per fishing day in two categories "in or around the MPA" and "remote from the MPA"

	Negative PPFD	Positive PPFD	TOTAL
	% row Size %column	% row Size %column	% row Size %column
	IN or AROUND THE MPA	16	112	128
	12,5%		87,5%	100,0%
		41,0%	60,2%	56,9%
	REMOTE FROM THE MPA	23	74	97
	23,7%		76,3%	100,0%
		59,0%	39,8%	43,1%
	TOTAL	39	186	225
	17,3%		82,7%	100,0%
		100,0%	100,0%	100,0%
	KHI2 = 4.09 / 1 DEGREES OF FREEDOM / PROBA ( KHI2 > 4.09 ) = 0.043 / TEST-VALUE = 1.72

Table 3 :

 3 Profit per fishing day in three categories "in or around the MPA" and "remote from the MPA"

	Negative PPFD	Medium PPFD	High PPFD	TOTAL
	% row Size %column % row Size %column % row Size %column % row Size %column
	IN or
	AROUND	16	104	8	128
	THE MPA
	12,5%		81,3%	6,3%	100,0%
		41,0%	61,5%	47,1%	56,9%
	REMOTE
	FROM	23	65	9	97
	THE MPA
	23,7%		67,0%	9,3%	100,0%
		59,0%	38,5%	52,9%	43,1%
	TOTAL	39	169	17	225
	17,3%		75,1%	7,6%	100,0%
		100,0%	100,0%	100,0%	100,0%
	KHI2 = 6.16 / 2 DEGREES OF FREEDOM PROBA ( KHI2 > 6.16 ) = 0.046 / V.TEST = 1.69

Funded by the French Pacific Fund (SPP)

A methodology to assess some socio-economic effects of marine protected areas


Combining ZDES with Immersed Boundary

Conditions technique for the treatment of complex geometries L. Mochel , P.-É. Weiss and S. Deck

ONERA -The French Aerospace Lab, F-92190, Meudon, France

The present paper focuses on a numerical strategy called ZIBC consisting of the zonal use of the Immersed Boundary Conditions combined with the ability of Zonal Detached Eddy Simulation [START_REF] Deck | Zonal-Detached-Eddy Simulation of the flow around a high-lift configuration[END_REF][START_REF] Deck | Recent improvements in the Zonal Detached Eddy Simulation (ZDES) formulation[END_REF] to simulate high Reynolds number separated flows. The motivation of such a strategy lies in the accurate handling of geometrically complex configurations with validated unsteady tools. Indeed, in the framework of scale-resolving unsteady simulations with a classical body-fitted mesh approach, the grid generation can become dramatically time-consuming or even unfeasible due to strict grid requirements. The aforementioned methodology appears to be a way to circumvent those limitations. In particular, it would be of considerable interest within the framework of space launcher afterbody by adding geometrical complexity to simplified configurations. Thus, we propose here to assess such a methodology on a simplified space launcher afterbody. The configuration consists of a cylinder elongated by another cylinder of smaller diameter (i.e. an extension). Immersed Boundary Conditions are used to handle the introduction of the extension into a pre-existing structured curvilinear grid fitting the ZDES grid requirements to treat the blunt body configuration (Fig. 1). The governing equations are solved using a standard body-fitted finite volume technique over the whole grid. A direct forcing source term, based on the work of [START_REF] Mohd-Yusof | Combined immersed-boundary/B-spline methods for simulations of flows in complex geometries[END_REF] is used on the cells inside the emergence to drive the velocity and the turbulence variables to the chosen values. Numerical simulations are performed at a Reynolds number based on the larger diameter of 6 10 2 . 1  and a free stream Mach number of 702 . 0

. The numerical results are compared with the experimental (respectively numerical) data of [START_REF] Deprés | Analysis of unsteadiness in afterbody transonic flows[END_REF][START_REF] Meliga | Unsteady transonic flow behind an axisymmetric afterbody equipped with two boosters[END_REF] (respectively [START_REF] Weiss | On the dynamics of axisymmetric turbulent separating/reattaching flows[END_REF]), corresponding to the classical body-fitted mesh approach. Fig. 2 and Fig. 3 suggest that the emergence is properly taken into account. Furthermore, the pressure coefficient (see Fig. 4) illustrates a very good agreement between both the available experimental and ZDES numerical data (classical approach without IBC) and the current ZIBC simulation.

In the final version of the paper, the influence of the combination of Immersed Boundary Conditions technique with ZDES on the fluctuating pressure field around the immersed body will be evaluated within the framework of high Reynolds number separated flows.

Fig. 1 Fig. 2 Fig. 3 Fig. 4

 1234 Fig. 1 Principle of ZIBC decomposed into classical primary approaches (i.e. ZDES and IBC)

Fifth HRLM Symposium, 19-21 March 2014, A&M University, College Station, USA Session:


Guillaume Guguin

Olivier Allix

Pierre Gosselet

Stéphane Guinard

Nonintrusive coupling of 3D and 2D laminated composite models based on finite element 3D recovery

Keywords: Plate-3D coupling, mixed dimensionality, nonintrusive coupling, laminated composites

In order to simulate the mechanical behavior of large structures assembled from thin composite panels, we propose a coupling technique which substitutes local 3D models for the global plate model in the critical zones where plate modeling is inadequate. The transition from 3D to 2D is based on stress and displacement distributions associated with Saint-Venant problems which are precalculated automatically for a simple 3D cell. The hybrid plate/3D model is obtained after convergence of a series of iterations between a global plate model of the structure and localized 3D models of the critical zones. This technique is nonintrusive because the global calculations can be carried out using commercial software. Evaluation tests show that convergence is fast and that the resulting hybrid model is very close to a full 3D model.

Introduction

This paper deals with nonintrusive techniques for the simulation of structures assembled from large flat composite panels. In many such structures, one dimension is much smaller than the other two, so they can be modeled as plates. However, plate theories alone cannot take into account the 3D stress concentration states which occur near edges, geometric discontinuities or zones subjected to concentrated loads. These phenomena are important in the case of laminated composites because they are often responsible for the initiation and propagation of damage, e.g. in the case of delamination [START_REF] Garg | Delamination -a damage mode in composite structures[END_REF]. Therefore, a natural approach to dealing with large composite structures is to attempt, whenever necessary, to couple plate or shell models with 3D models [START_REF] Düster | Applying the hp-d version of the FEM to locally enhance dimensionally reduced models[END_REF].

The main topic of this paper is the development and evaluation of a type a coupling between a plate model and a 3D model which could be applied directly within the nonintrusive framework proposed in [START_REF] Gendre | Non-intrusive and exact global/local techniques for structural problems with local plasticity[END_REF][START_REF] Gendre | A two-scale approximation of the Schur complement and its use for nonintrusive coupling[END_REF] as well as in [5] for dynamics. Essentially, the nonintrusive coupling technique consists, starting from an initial global calculation, in alternating between local calculations with prescribed displacements at the boundaries and global calculations which include corrective loads in order to reduce the imbalance between the models. Because the nonintrusive coupling technique used in the paper was designed for models with the same dimensionality, the main problem addressed here is the recovery of the 3D quantities associated with the 2D model.

Overall, the question of the link between a plate model and its 3D counterpart is of such practical and theoretical importance that it has given rise to many publications. This introduction focuses more specifically on the references of interest in the literature that we tried to incorporate into our formulation.

The quality of a plate model is usually assessed by comparison with a corresponding 3D reference model. This has led to a first group of works focusing on the development of a posteriori error estimators, initially neglecting edge effects [START_REF] Koiter | On the foundations of the linear theory of thin elastic shells[END_REF][START_REF] Danielson | Improved error estimates in the linear theory of thin elastic shells[END_REF][START_REF] Simmonds | An improved estimate for the error in the classical linear theory of plate bending[END_REF][START_REF] Koiter | Foundations of shell theory[END_REF][START_REF] Ladevèze | Validity of the classical linear shell theory[END_REF], then taking them into account [START_REF] Ladevèze | On the validity of linear shell theories[END_REF]. Thus, by the end of the 1980's, the performance of the classical Kirchhoff-Love and Reissner-Mindlin plate and shell theories in 3D elasticity was relatively well-understood. Another group of works concerned the technique of asymptotic developments with the thickness as the small parameter [START_REF] Friedrichs | A boundary layer theory for elastic bending of plates[END_REF][START_REF] Goldenveizer | The principles of reducing three-dimensional problems of elasticity to two-dimensional problems of the theory of plates and shells[END_REF][START_REF] Ciarlet | A justification of the two-dimensional plate model[END_REF][START_REF] Destuynder | Une méthode asymptotique des plaques minces en élasticité linéaire[END_REF][START_REF] Destuynder | Some theoretical aspects in the modelling of delamination for multilayered plates -Local effects in the analysis of structures[END_REF]6,[START_REF] Ciarlet | Plates and junctions in elastic multi-structures : An asymptotic analysis[END_REF]. In addition, using the same approach as Reissner [START_REF] Reissner | The effect of transverse shear deformation on the bending of elastic plates[END_REF][START_REF] Reissner | On a mixed variational theorem and on shear deformable plate theory[END_REF] in statics and Mindlin in dynamics [START_REF] Reissner | Influence of rotatory inertia and shear in flexural motions of isotropic elastic plates[END_REF], a shear factor was often introduced [START_REF] Whitney | Shear correction factors for orthotropic laminates under static load[END_REF] in order to correct what would have been deduced from a pure kinematic assumption by taking into account some a priori information on the three-dimensional stress state. This mixed nature of the link between plate theory and 3D theory is the reason why the construction of improved theories [START_REF] Lure | Three dimensional problems of the theory of elasticity[END_REF][START_REF] Touratier | An efficient standard plate theory[END_REF][START_REF] Levinson | An accurate simple theory of the statics and dynamics of elastic plates[END_REF][START_REF] Reddy | A simple higher-order theory for laminated composite plates[END_REF][START_REF] Carrera | Theories and finite elements for multilayered, anisotropic, composite plates and shells[END_REF] and associated finite elements [START_REF] Carrera | On the use of the Murakami's zig-zag function in the modeling of layered plates and shells[END_REF][START_REF] Auricchio | A mixed-enhanced finite-element for the analysis of laminated composite plates[END_REF] relies heavily on mixed formulations.

Regardless of the method, the assessment of a plate model involves the comparison between a 3D "Saint-Venant" or long-wavelength solution [START_REF] Ladevèze | On the saint-venant principle in elasticity[END_REF][START_REF] Gregory | On plate theories and Saint-Venant principle[END_REF] and the solution obtained from plate theory. A key aspect of the comparison, especially in the case of composites, concerns the distribution of the stress field [START_REF] Pagano | Stress fields in composite laminates[END_REF][START_REF] Davi | Stress fields in general composite laminates[END_REF]. Therefore, a basic technique which is used, sometimes implicitly, when estimating the quality of plate theories is the recovery of equilibrated 3D solutions by integrating the equilibrium equations through the thickness [START_REF] Rohwer | Calculating 3D stresses in layered composite plates and shells[END_REF][START_REF] Allix | New Trends in Thin Structures: Formulation, Optimization and Coupled Problems, volume 519[END_REF]. The application of this technique to finite element calculations requires the recovery of 2D fields which satisfy the 2D equilibrium equations exactly. The use of such a stress recovery technique [START_REF] Noor | Predictor Corrector procedures for stress and free-vibration analysis of multilayered composite plates and shells[END_REF][START_REF] Daghia | A hybrid stress approach for laminated composite plates within the First-order Shear Deformation Theory[END_REF] makes the implementation within a commercial program difficult. A comparative analysis of various models and techniques to evaluate a priori or a posteriori out-of-plane stresses was presented in [START_REF] Carrera | A priori vs. a posteriori evaluation of transverse stresses in multilayered orthotropic plates[END_REF]. An improved hybrid post-processing procedure was proposed in [START_REF] Cen | A new hybrid enhanced displacement based element for the analysis of laminated composite plates[END_REF] along with a new displacementbased element. This procedure, which works element-wise and, thus, is easier to implement, was extended to first-order theories in [START_REF] Kim | Enhanced first-order theory based on mixed formulation and transverse normal effect[END_REF] to recover both displacements and stresses.

In order for the 3D approximation derived from plate theory to be valid everywhere, it should be completed by adding an approximation of the edge effects. A key issue is that the problem deduced from the plate theory residual at the edges must lead to a localized problem [START_REF] Ladevèze | On the validity of linear shell theories[END_REF][START_REF] Ladevèze | A new version of the Reissner-Mindlin's plate theory for orthotropic homogeneous plates[END_REF][START_REF] Allix | New Trends in Thin Structures: Formulation, Optimization and Coupled Problems, volume 519[END_REF]. Otherwise, even using what one calls refined theories, the plate's interior solution would be meaningless for the orders of magnitude considered. This may be one of the reasons why the Reissner-Mindlin theory is the dominant approach today, even for composite structures. Another reason is, of course, that because it involves only first-order derivatives it is much easier to build finite elements based on this theory than, for example, on Kirchhoff-Love's. However, the Reissner-Mindlin theory includes an artificial edge effect which, if not handled properly, could pollute the recovery of the 3D solution near the edges [START_REF] Rössle | Corner singularities and regularity results for the Reissner/Mindlin plate model[END_REF].

In this context, our approach to nonintrusive coupling must comply with several constraints. The first constraint is that our method involves two overlapping models, a 2D (plate or shell) model which is defined over the whole domain to be analyzed and remains unchanged throughout the analysis, and a refined 3D model in which structural details and complex nonlinearities such as contact or delamination can be taken into account. The second constraint is that we need to use results obtained from commercial finite element codes which, usually, are based on the Reissner-Mindlin model and do not provide recovery techniques. The implementation of such recovery techniques would require rather complex manipulations with simultaneous access to data from several elements, which would annihilate the very purpose of nonintrusive techniques. Moreover, the plate model is usually less than optimal compared to its 3D counterpart in that the plate quantities which would be deduced from the calculation of the full 3D model would be somewhat different from those obtained using the plate model, even far away from the edges. That is the reason why [START_REF] Noor | Predictor Corrector procedures for stress and free-vibration analysis of multilayered composite plates and shells[END_REF][START_REF] Yu | Asymptotically accurate 3-D recovery from Reissner-like composite plate finite elements[END_REF] and others have used an asymptotic formulation up to the second order to define a quasi-optimal version of Reissner-Mindlin's composite plate theory by means of an optimization procedure. Then, based on these accurate plate quantities, a recovery technique can be used to obtain the 3D displacement, strain and stress fields as functions of 2D variables.

The paper is organized as follows: Section 2 shows how relevant 3D quantities can be generated for any stacking sequence in order to recover 3D stress and displacement bases. Then, in Section 3, these bases are used to define hybrid models. In Section 4, a first analysis of the results enables us to identify the most favorable coupling technique, which is tested on a more representative example in Section 5.

Numerical recovery of the Saint-Venant stress and warping

In this section we investigate the association of 3D quantities (traction and warping) with composite plate calculations.

Notations and equations

Let us consider a classical 3D problem defined in domain Ω. Let f be the volume force, g the traction prescribed over ∂ N Ω and u D the prescribed displacement over the complementary part ∂ D Ω. The unit vectors of an orthogonal spatial basis are denoted by (e x , e y , e z ). The linear elastic behavior is given by Hooke's tensor H. The symmetric part of the displacement gradient and the Cauchy stress tensor are designated respectively by ǫ(u) and σ.

We introduce the affine space of the admissible displacements U(Ω):

U(Ω) = u ∈ H 1 (Ω), u = u D on ∂ D Ω
whose associated vector space is U 0 (Ω). The 3D problem is:

         Find u ∈ U(Ω) such that, ∀u * ∈ U 0 (Ω), Ω σ : ǫ(u * )dx = Ω f • u * dx + ∂N Ω g • u * dS σ = H : ǫ(u) (1)
Regarding the plate model, we have

Ω = ω × [-h, h].
Using the notation ˜to designate in-plane quantities (x and y coordinates), the first assumption concerns the shape of the plate's displacement:

v = ṽ + v z e z + z θ ∧ e z (2)
where ṽ and v z are respectively the in-plane and out-of-plane displacements of the mid-surface and θ denotes the rotation of the section, all of which are functions defined in ω. Operator ∧ is the cross product defined by:

  x 1 x 2 x 3   ∧   y 1 y 2 y 3   =   x 2 y 3 -x 3 y 2 x 3 y 1 -x 1 y 3 x 1 y 2 -x 2 y 1  
This leads to the introduction of the classical plate strains:

ǫ(v P ) = ǫ ǫz sym v z,z with      ǫ = γ(ṽ) + z χ( θ) ǫz = ∇(v z ) + θ ∧ e z 2 (3)
γ(ṽ) = v x,x vx,y+vy,x 2 sym v y,y , χ( θ) = θ y,x θy,y-θx,x 2 sym -θ x,y
Using notation for the integration through the thickness, the plate's stress tensors are:

Ñ = σ M = z σ Q = σz (4)
The second assumption of the plate model is that the peeling stress σ zz is negligible compared to the other components of the stress tensor. This assumption simplifies the expression of the plate's behavior as a function of the plane stress Hps and out-of-plane stress B. (In 3D, σz = Bǫ z .)

For a plate in equilibrium, the constitutive relations are: Ñ = Hps : γ(ṽ)

M = z 2 Hps : χ( θ) Q = κ B : ǫz (v z , θ) (5)
where κ is a shear correction factor. (In this paper, we use the correction factor from [START_REF] Batoz | Modélisation des structures par elements finis[END_REF].)

The admissible plate displacement field is:

U P (ω) = (ṽ, θ, v z ) ∈ H 1 (ω) 5 , which satisfies Dirichlet's BCs (6)
and the plate problem becomes1 :

           Find (ṽ, θ, v z ) ∈ U P (ω) such that ∀(ṽ * , θ * , v * z ) ∈ U P 0 (ω) ω Ñ : γ(ṽ * ) + M : χ( θ * ) + Q • ǫz (v * z , θ * ) dx = Ω f • v * dx + ∂N Ω g • v * dS Ñ = Hps : γ(ṽ), M = z 2 Hps : χ( θ), Q = κ B : ǫz (v z , θ) (7)
Remark 1. The Reissner-Mindlin kinematic assumption was chosen because it is the most commonly used in commercial software. This model provides an estimation of the transverse shear forces, but, unfortunately, it leads to perturbations of the generalized forces near the edges [START_REF] Rössle | Corner singularities and regularity results for the Reissner/Mindlin plate model[END_REF].

Remark 2 (Notations). From now on, we will use the following notations: for the kinematics:

plate DOFs V = v x v y θ x θ y v z T plate kinematics V = e x e y -ze y ze x e z so that v = 5 k=1 V k V k
and for the generalized forces:

F = N xx N xy M xx M xy Q x N yy M yy Q y T so that F 1:5 = σ : (V 1:5 ⊗ e x ) F 6:8 = σ : (V 2,4,5 ⊗ e y ) (8)
Although these notations seem biased toward the (e x , e y ) coordinate system (especially toward direction e x since the components related to this normal direction are listed first), the following calculations restore the normal-independent property of the method.

In order to deal with interfaces which are not aligned with the axes, we introduce the following notations: for an interface whose local coordinate system is (n, t, e z ) (n being the normal vector and t the tangent vector), φ denotes the angle such that cos(φ) = n • e x = n x . The kinematic basis of the plate in the local coordinate system is V φ = (n, t, -zn, zt, e z ). The associated DOFs are V φ = (v n , v t , θ n , θ t , v z ) and the stress components are

F φ = N nn N nt M nn M nt Q n N tt M tt Q t T .

Recovery of the Saint-Venant traction fields

Since plate solutions are relevant only for problems with long characteristic lengths of variation, we try to associate one 3D Saint-Venant stress field (τ i ) with each plate stress component.

In order to do that, we consider the simplest 3D problems whose solutions are Saint-Venant fields, which consist in a sufficiently large square plate with the same material, the same stacking sequence and the same finite elements as the local model, subjected to 8 components of generalized forces (3 in tension, 3 in bending and 2 in shearing). This leads to 8 basic problems with well-chosen regular loads, each activating one generalized force as defined in Table 1 and illustrated in Figures 12345. We call these problems "cell problems" because of the analogy of the method with micro-macro approaches. Here, we are less interested in the average response than in the distribution of the 3D mechanical quantities through the thickness in the center of the plate, which represents a numerical "3D recovery" of the plate quantities. j loading main load figure

1 σ xx = -1 over ∂Ω 0 , σ xx = 1 over ∂Ω L N xx 1 2 σ xy = -1 over ∂Ω L σ xy = 1 over ∂Ω 0 σ xy = -1 over ∂Ω a σ xy = 1 over ∂Ω -a N xy 2 3 σ xx = z over ∂Ω 0 , σ xx = z over ∂Ω 0 M xx 3 4 σ xy = -z over ∂Ω 0 σ xy = z over ∂Ω L σ xy = z over ∂Ω -a σ xy = -z over ∂Ω a M xy 4 5 σ xz = 1 over ∂Ω L Q x 5 6
same as Problem 1, but rotated 90 From each calculation j ∈ 1, 8 , one can extract the distribution of the 3D stress tensor (σ j ) through the thickness in the center of the plate, i.e. sufficiently far away from any edge effect. This is illustrated in Figure 6 after FE discretization. From this stress tensor, using Equation ( 8), one can post-process the generalized forces (F j ) of the plate. These test cases were designed so that F j j is dominant in the j th experiment (although coupling among generalized forces exists naturally, if only because of the conservation of moments which involves shear forces).

Because of the linearity of the problem, any Saint-Venant stress state is a combination of the extracted (σ j ). The (τ i ) are constructed such that the contributions of the generalized forces are uncoupled:

8 i=1 F i j τ i = σ j =⇒ τ 1 . . . τ 8 = σ 1 . . . σ 8 F i j -1 (9)
where brackets denote matrices. Then, the (τ i ) are obtained by inversion of the 8 × 8 matrix F i j . In other words, the (τ i ) belong to the space of the Saint-Venant stresses and satisfy the orthogonality relation:

1 = τ i : (V k ⊗ e v ) for v = x and i = k ∈ 1, 5 or for v = y and (i, k) ∈ {(6, 1), (7, 3), (8, 5)} 0 = τ i : (V k ⊗ e v ) in all other cases (10)
Figure 7 shows typical distributions of the nodal forces through the thickness, both for an isotropic plate and for 4 orthotropic plies. (A complete description of the material will be presented in Section 4.) In the isotropic case, we recover the analytical functions of [START_REF] Mccune | Mixed-dimensional coupling in finite element models[END_REF].

Remark 3. The (τ i ) were designed for loads which are aligned with axes (e x , e y ). If necessary, one can substitute a tensor family (τ φ i ) in a rotated direction: given a local coordinate system (n, t), let n = n x e x + n y e y so that t = -n y e x + n x e y . Using linearity, the Saint-Venant solution for a load which is aligned with the local frame is:

τ φ Nnn = n 2 x τ Nxx + n x n y τ Nxy + n 2 y τ Nyy τ φ Nnt = -n x n y τ Nxx + (n 2 x -n 2 y )τ Nxy + n y n x τ Nyy τ φ Ntt = n 2 y τ Nxx -n x n y τ Nxy + n 2 x τ Nyy τ φ Qn = n x τ Qx + n y τ Qy , τ φ Qt = -n y τ Qx + n x τ Qy (11)
Obviously, the same transformation can be used for moments as well as for tensions. Remark 4. The quality of the recovered Saint-Venant stress fields can be improved iteratively: the calculated stress fields are applied as traction boundary conditions to the cell problem, resulting in a problem with even more localized edge effects and a better stress distribution through the thickness in the center of the domain. Figure 8 shows that the stress field converges very rapidly because each correction is an order of magnitude smaller than the previous one.

Recovery of the Saint-Venant warping vectors

Now let us try to define relevant information concerning the displacement field. Since the plate's kinematics (rigid section) is unrealistic and would lead to significant errors [START_REF] Allix | New Trends in Thin Structures: Formulation, Optimization and Coupled Problems, volume 519[END_REF], we propose to correct this kinematics using warping vectors (i.e. 3D displacements which are allowed to vary through the thickness) deduced from the Saint-Venant auxiliary problems.

From the eight numerical experiments, we obtain the stress states (σ j ) and the displacement profiles (u j ). By inverting the matrix of the generalized forces of the plate [F i j ], one obtains the family (τ i ) of the uncoupled stresses. The same matrix can be used to define the displacements U i associated with pure Saint-Venant loads:

U 1 . . . U 8 = u 1 . . . u 8 F i j -1 (12)
Each of these displacements consists of a plate part and a warping part. In order to separate the two contributions, we assume that the plate displacements produce the same work as the 3D displacements in the Saint-Venant tractions. At each point of interface γ, one has:

V k i = τ k : (U i ⊗ e x ) (i ∈ 1, 8 , k ∈ 1, 5 ) w i = U i - 5 k=1 V k V k i ( 13
)
Figure 9 presents some of the warping functions extracted from the numerical 3D recovery. The result is in the form of classical zigzags [START_REF] Carrera | Theories and finite elements for multilayered, anisotropic, composite plates and shells[END_REF] with problem-specific shapes.

Remark 5. The extraction of the five kinematic plate components is carried out using only the first five Saint-Venant stress tensors e x : (τ k • e x ) k=1:5 in the normal direction. Nevertheless, this definition is sufficient for the warping to be well-defined (apart from negligible terms) and orthogonal to all eight Saint-Venant stress tensors for any normal direction. Indeed, from [START_REF] Allix | New Trends in Thin Structures: Formulation, Optimization and Coupled Problems, volume 519[END_REF], we know that for Saint-Venant problems the plate displacement corrected by warping in direction e z results in 3D stresses which are O(h/L)-accurate, 2h being the thickness of the plate; additional warping in the plane leads to O(h 2 /L 2 )-accuracy. This leads to the following orders of magnitude (which, in practice, turn out to be very pessimistic because the numerical measurements are much smaller):

(τ 1 • e x -τ 2 • e y ) • U i = (τ 2 • e x -τ 6 • e y ) • U i = o(h 2 /L 2 )( Ñ j : γi ) (τ 3 • e x -τ 4 • e y ) = (τ 4 • e x -τ 7 • e y ) • U i = o(h 2 /L 2 )( M k : χi ) (τ 5 • e x -τ 8 • e y ) • U i = o(h/L)(Q l • ǫz,i )        ∀i ∈ 1, 8 ∀j ∈ {1, 2, 6} ∀k ∈ {3, 4, 7} ∀l ∈ {5, 8} (14)
These relations imply that the warping vectors are independent of the axes chosen and produce zero work in all Saint-Venant stress fields, even those that are rotated (τ φ k ) because of linearity.

-1 -0.5 0 0.5 1 In this section, we focus on the definition of hybrid models for thin structures using plate theories wherever appropriate, but keeping 3D models when necessary (see Figure 10). In our context of using commercial software, only less-thanoptimal plate theories are available and, thus, a direct plate/3D connection would be error-prone. Besides, the analysis of an industrial structure often begins with a global plate model, so we want to take advantage of this model to the greatest possible extent and supplement it with local 3D modeling only where necessary. By convention, we use lowercase for the plate domain and for the interfaces (ω, γ) and uppercase for their 3D counterparts (Ω = ω × [-h, h], Γ = γ × [-h, h]). We distinguish three zones: the zone of interest I, the buffer zone B and the complementary zone C (which exists only in the plate model).

×10 -6 -0.2 -0.1 0 0.1 0.2 w Q x x thickness z orthotropic isotropic (
Then, the hybrid 3D plate model is defined as the limit of an iterative process, illustrated in Figure 11, which alternates between global plate calculations (domain First, a global plate calculation is used to define the boundary conditions for the local 3D problems along γ C (the zooming step). As explained in Section 3.1, the plate-to-3D recovery uses the stress and warping distributions through the thickness deduced from the cell problems of Section 2.

ω = ω I ∪ ω B ∪ ω C ) and local 3D calculations (domain Ω I ∪ Ω B ).
After the local calculation, the imbalance, a plate sense, between the plate model and the 3D model is evaluated along the inner line γ I (Section 3.2). Then, if necessary, a global plate calculation under a corrective loading applied along γ I is carried out (Section 3.3).

Beside the dimensionality transition problem, our approach differs from the algorithm of [START_REF] Gendre | Non-intrusive and exact global/local techniques for structural problems with local plasticity[END_REF] in that we propose to distinguish interface γ C (which serves to prescribe the boundary conditions of the zone of interest) from interface γ I (which is used to apply the global correction) and introduce a buffer zone Ω B in-between. For the descent steps described in Sections 3.1.3 and 3.1.2, it is possible to choose γ C = γ I , but spurious edge effects may occur near Γ C . The buffer zone Ω B is used to let the evanescent components activated by the boundary conditions become negligible before entering the of interest Ω I itself. In a sense, the plate model and the 3D model are made equivalent in the buffer zone, which makes this method different from other approaches based on a volume connection zone, such as the transition element method [START_REF] Garusi | A hybrid stress-assumed transition element for solid-to-beam and plate-to-beam connections[END_REF][START_REF] Schiermeier | Interface elements in global/local analysispart 3: Shell-to-solid transition[END_REF] or the Arlequin method [4]. Besides, the implementation of the method is simpler thanks to its compatibility with nonintrusive approaches.

The zooming step

In this section, we assume that a plate problem has been solved and that all its kinematic and static quantities are known along line γ C . Let (F i ) i∈ 1,8 and (V k ) k∈ 1,5 denote respectively the 8 static plate components and the 5 displacement components at the interface, and let s be the curvilinear abscissa of interface γ C . Thus, the position of a point of Γ C is given by (s, z) and the classical 3D plate displacement is v

P (s, z) = k V k (z)V k (s).
The objective of the zooming step is to define relevant boundary conditions along Γ C = γ C × [-h, h] in order to solve the 3D problem in the extended zone of interest Ω I ∪ Ω B (see Figure 10). In order to do that, we use the stress and warping distributions (τ i ) and (w i )s calculated for the cell problems of Section 2.

Pure traction descent

At each point of interface Γ C , the following traction distribution is prescribed:

σ(s, z) • n(s) = 8 i=1 F i (s)τ i (z) • n(s) along Γ C (15)
Then, the 3D problem becomes:

Find u ∈ U(Ω I ) such that ∀u * ∈ U 0 (Ω I ) ΩI H : ǫ(u) : ǫ(u * )dx = ΩI f • u * dx + ∂N ΩI g • u * dS + ΓC 8 i=1 F i (s)(τ i (z) • n(s)) • u * dsdz (16)
Of course, this technique works only for domains of interest possessing sufficient Dirichlet conditions. Therefore, it will be evaluated only briefly.

Lagrangian descent

The Lagrangian technique is a kinematic coupling method based on an energy equivalence which is similar to that proposed in [START_REF] Mccune | Mixed-dimensional coupling in finite element models[END_REF][START_REF] Shim | Mixed dimensional coupling in finite element stress analysis[END_REF] for isotropic plates with known analytic solutions in order to calculate the distribution of Saint-Venant stresses (τ φ i ). Here, since we are considering stacks of orthotropic plies, no such formula is available and, therefore, we propose to use the numerically calculated (τ φ i ). The zooming step consists in prescribing the 3D displacement at the interface which develops the same amount of work as the given plate displacement with the Saint-Venant stresses (τ

φ i ). Find u ∈ U(Ω I ), (λ k ) ∈ Λ(γ C ) 5 such that ∀u * ∈ U 0 (Ω I ) ΩI H : ǫ(u) : ǫ(u * )dx = ΩI f • u * dx + ∂ f ΩI g • u * + 5 k=1 γC λ k (s) z u * • τ φ k (z) • n(s) dz ds + 8 i=6 ΓC F φ i (s)(τ φ i (z) • n(s)) • u * dzds γC λ * k (s) z (u -v P ) • τ φ k (z) • n(s) dz ds = 0, ∀k ∈ 1, 5 , ∀λ * k ∈ Λ(γ C ) (17)
where Λ(γ C ) is the space of the Lagrange multipliers (λ k ).

In other words, the 5 kinematic degrees of freedom are made consistent with their 3D counterparts. The plate part of the 3D displacement, like warping in Section 2.3, is defined using averages weighted by (τ φ k • n). The Lagrange multipliers (λ k ) are the magnitudes of the normal forces and moments due to the reaction of the 3D domain to the prescribed displacements:

k ∈ 1, 5 , z σ(s, z) • n(s) • V φ k (z)dz = λ k (s) z τ φ k (z) • n(s) • V φ k (z)dz = λ k (s) (18)
The information concerning the last three stress components 6, 8 , which involve only tangential action, is given by a classical Neumann condition. This Lagrangian coupling has the advantage of being applicable to zones of interest with no Dirichlet condition. It also makes the plate's stress components in the 3D zone of interest directly available, which is useful in certain configurations of the global correction step (typically when no buffer zone is used). Unfortunately, this coupling requires the local orientation of the interface to be taken into account, which makes it hard to implement for curved domains. (In the above equations, φ depends implicitly on s.) Also, it is difficult to choose a proper space for multiplier Λ(γ C ) in the case of nonconforming meshes: the LBB condition must be satisfied (see Section 3.5.2).

Pure displacement descent

Pure displacement descent is a very easy coupling to set up. The principle is to transmit as much information as possible using only Dirichlet conditions at the interface. This eliminates the potential instability problems of Lagrangian coupling and makes the coupling independent of the local orientation of the interface. The prescribed 3D displacement at the interface is given by:

u(s, z) = k=1:5 V k (z)V k (s) + i=1:8 w i (z)F i (s), s ∈ γ C , z ∈ [-h, h] (19)
In other words, the classical plate displacements are enriched by warping vectors weighted by the generalized force components of the associated plate.

Evaluation of the residual

Now let us assume that a 3D problem has been solved through a Lagrangian or pure displacement descent step. This leaves us with a plate solution in ω = ω I ∪ ω B ∪ ω C and a 3D solution in Ω I ∪ Ω B with kinematic continuity along γ C (in a sense which depends on the descent chosen). The residual (L φ k ) is defined as the imbalance (in a plate sense) between the plate and the 3D domains along line γ I : F φ,P k (s) : generalized plate forces along γ I viewed as the border of

ω C ∪ ω B F φ,3D k (s) = σ : (n ⊗ V k ) : 3D generalized forces along γ I viewed as the border of ω I L φ k = F φ,3D k + F φ,P k along γ I , k ∈ 1, 5 (20)
The residual is evaluated in a plate sense because the intention is to reuse it as a corrective load applied to the plate model. A 3D residual could be defined through a recovery of the plate forces, but, besides being unsuitable for the global/local strategy, it would be difficult to evaluate because an error of order h L is expected in the σ xz and σ yz components and an error of order h 2 L 2 is expected in the σ zz component. The fact that the 3D quantities have the correct orders of magnitude is interesting information which contributes to the validation of the hybrid model.

The global correction step

The global correction step consists in solving the global plate problem after having updated the loading by residual (R k ) along interface γ I . The updating takes the following form:

5 k=1 γI L φ k V φ, * k ds (21)
where V φ, * k are the degrees of freedom of the plate test displacement field in the coordinate system of γ I . The iterations converge as long as the plate model of the zone of interest ω I is stiffer than the 3D model Ω I (classical fixed point). Otherwise, relaxation is necessary. Such an algorithm also converges for nonlinear monotonic problems [START_REF] Gendre | Non-intrusive and exact global/local techniques for structural problems with local plasticity[END_REF].

In the linear case, the convergence can be improved, as is classically done in Schwarz methods, by using Krylov acceleration [START_REF] Widlund | Domain decomposition methods -algorithms and theory[END_REF]. In this case, the full plate model (ω) can be interpreted as a preconditioner to the hybrid model described in the following section.

In cases where the zooming step has sufficient Dirichlet conditions (i.e. the Lagrangian and pure displacement cases), the bufferless limit case γ I = γ C is acceptable. In the Neumann case, a nonzero buffer is necessary. In any case, overlapping is interesting with all methods because it reduces the small edge effects generated by the interface loading of the zooming step.

Characterization of the limits

At the convergence of the iterative process, the residual is zero and the hybrid model can be interpreted as follows:

• In the zone of interest Ω I , it is a pure 3D model,

• In the complementary zone ω C , it is a pure plate model,

• Along interface γ I , the stresses are in equilibrium,

• Along interface γ C , the plate's kinematic continuity is satisfied (along with some other relations which depend on the zooming technique chosen),

• In the buffer zone, the 3D problem in Ω B is equivalent to the plate problem in ω B in the following sense:

ωB σ : (V φ k ⊗ e φ ) -F P,φ k V * k dx = 0 ∀k ∈ 1, 5 , ∀φ, ∀V * k σ • (V φ k ⊗ n) -F P,φ k = 0 on γ I , ∀k ∈ 1, 5 u • (τ φ k • n) -V P,φ k = 0 on γ C , ∀k ∈ 1, 5 (22)
In other words, the plate model obtained by integrating the 3D model through the thickness in Ω B must be identical to the plate model in ω B .

One should note that in all cases 3D displacements orthogonal to the Saint-Venant tractions are likely to develop in Ω B , but these should be evanescent and should have negligible effect on Ω I .

Implementation into a finite element code

Preprocessing of the Saint-Venant stresses and warping vectors

This preprocessor requires an auxiliary 3D problem to be built with the same stacking sequence as the original problem, but a simpler geometry. Minimal Dirichlet conditions are prescribed, except for shear-dominated experiments (j = 5 and j = 8). An important property is that because the plate's kinematics (spanned by functions (V k )) is linear through the thickness, it is represented exactly by the 3D FE approximation. This means that the calculation of the plate's generalized stresses can be performed directly on the nodal reactions (which are classical outputs of FE software). Let (F l ) be the nodal reactions through the thickness calculated on a face with normal e x . If z l is the height of the l th node), one has:

N xx = σ xx = l F l x N xy = σ xy = l F l y M xx = zσ xx = l z l F l y M xy = zσ xy = - l z l F l x Q x = σ xz = l F l z (23)

Issues related to the Lagrangian approach

The implementation of the Lagrangian approach requires great care in choosing the finite element approximation spaces.

In the case of matching meshes and interpolations (typically 8-node hexahedra on the 3D side and linear plate elements on the plate side), each plate node faces a column of 3D nodes across the thickness. A simple choice consists in using five multiple-point constraints (MPCs) per plate node and corresponding 3D nodes.

The case of nonmatching meshes or interpolations raises the classical issues of mixed approaches: one must check the LBB-stability of the scheme with a well-chosen discretization of the Lagrange multiplier field along the interface.

For example, for a linear plate element and a bilinear 3D element (20-node hexahedron), one method is to use two plate elements facing one column of 3D elements (so each plate node faces a column of 3D nodes), then use a constant Lagrange multiplier for each column of 3D elements.

Numerical study of the various zooming steps

In order to illustrate the descent step, let us consider the cantilever plate of Figure 12, with L = a = 5 mm, h = 0.2 mm and thickness 2h. The plate consists of four orthotropic plies (-45/45) s with the following mechanical properties: The prescribed displacement at the edge of the plate is u d • e z = 2h. The zone of interest is the clamped side because this is where 3D localization takes place. The 3D model was meshed with quadratic hexahedral HEXA20 elements (20 nodes). A regular mesh was used with 4 elements per ply in the z-direction and 20 elements in the x and y directions. The global model defined over the entire plate was obtained using Dhatt and Batoz homogenization [START_REF] Batoz | Modélisation des structures par elements finis[END_REF] and a discrete shear formulation. The reference is a full 3D model over the whole domain. The calculations were carried out using Python scripts interfaced with Code Aster.

E L = 25 MPa E T = E N = 1 MPa G LT = G LN = 0.5 MPa G T N = 0.2 MPa ν LT = ν T N = ν LN = 0.25
Since the zone of interest is clamped, it was possible to use all three zooming techniques. Figure 13 shows the evolution of some components of the 3D stress along the x axis. The reference is defined in Ω (x ∈ [0, 10]) while the 3D parts of the hybrid models are defined in zone x ∈ [0, 5]. In all three methods, as expected theoretically [START_REF] Allix | New Trends in Thin Structures: Formulation, Optimization and Coupled Problems, volume 519[END_REF], σ xz has approximately the same order of magnitude as σ xx ( h L ) and σ zz has approximately the same order of magnitude as σ xx ( h 2 L 2 ). One can observe that in all cases a perturbation appears near the interface (x = 5). However, since the generalized forces transmitted at the interfaces are correct, this local effect decreases rapidly and vanishes within a 2h band. The 3D solutions are close to the reference solution. Altogether, the perturbations caused by the warping technique seem smaller than those caused by the Neumann technique, which, in turn, are smaller than with Lagrangian coupling. The presence of these undesirable stress concentrations is what triggered the extension of the zone of interest by a 3D buffer zone in order to dampen the artificial edge effects.

Figure 14 shows the evolution of the generalized stresses along the interface for the initial plate calculation and for the 3D calculations. By construction, Neumann zooming produces the same generalized stresses as the plate calculation. Conversely, the other zooming techniques lead to slightly different results, which justifies the iterations between the local and global models. As expected, the differences between the plate model and the zooms are particularly significant near the edges. In this section, we study the application of the iterative coupling technique in order to achieve convergence toward the hybrid model. Based on previous experiments, we chose the warping-based descent algorithm (see 3.1.3) and used a nonzero buffer zone (see 3.3).

The objective was to model the holed composite plate shown in Figure 15. The local model consisted in a 3D representation of the hole and its vicinity with enough extent to capture the 3D edge effect. The global model was a plate, in which we could choose to represent (or not) the hole in the zone ω I corresponding to the 3D model. We used the same orthotropic layers as in the previous section. The dimensions were L = a = 60 mm and h = 1 mm. The radius of the hole was r = 2 mm. The loading was a transverse displacement u d = 0.45 along the right side.

The zone of interest ω I was a square of side 20 mm and the width of the buffer zone ω B was 2 mm. It was meshed with 24 elements per side and 16 elements in the thickness.

Convergence toward the hybrid model

Let us define the relative norm η i of the residual L of Section 3.2 as follows:

η 2 i = γI |( Ñ 3D i + Ñ P i ) • n| 2 ds max γI | Ñ P 0 • n| 2 + γI |( M 3D i + M P i ) • n| 2 ds max γI | M P 0 • n| 2 + γI |( Q3D i + QP i ) • n| 2 ds max γI | QP 0 • n| 2
Figure 16 shows the decrease in η i as a function of the number of iterations for the classical fixed point algorithm and for the conjugate gradient algorithm. The complexity of the plate model in the zone of interest ω I is a parameter of the method which influences the convergence, but not the limit. We considered two cases, depending on whether the hole was represented in the plate model ω I or not. When the hole was included, the plate model gave a better representation of the zone of interest, leading to much less initial imbalance. In both cases, the convergence was fast. In particular, one can observe that two conjugate gradient iterations of the model without a hole led to a similar quality approximation as the classical submodeling technique applied to the holed model (which is more difficult to generate and to mesh).

Validation of the hybrid model

Now let us study the hybrid model obtained after a sufficient number of iterations of the coupling procedure and compare the results of this model with those of a full 3D calculation of Ω in terms of the reliability of the 3D quantities obtained.

The quantities of interest chosen are the stresses in the vicinity of the hole, which would initiate delamination. Figures 17(a,b,e) show the variations of some stress components along the particular line described in Figure 17f. A comparison of the full 3D reference values, the values using submodeling (which is equivalent to the hybrid model with a single descent and no iteration) and the values given by the hybrid model after 5 iterations clearly shows that submodeling is very unreliable while the hybrid model is almost identical to the reference model.

Figures 17(c,d) show the evolution of the generalized forces along the same line. A comparison of the reference 3D calculation, the hybrid model at iterations 0 and 5 and a direct (holed) plate model shows that both the plate and the submodel lead to poor results while the hybrid model comes very close to the reference.

Table 2 quantifies the changes induced by the iterations in the plate quantities in the plate domain: the final values of the deflection, shear stress and moment (i.e. the values corresponding to the hybrid model) at a point of γ I (x = 20 mm, y = 0 mm) are compared with the initial values obtained with the plate model (with or without a hole). Of course, the model without a hole had significant errors, but even the model with a hole was improved by taking into account 3D phenomena. (In particular, the shear stress was underestimated by more than 5%.) -with a hole without a hole error in w (%) -1.93

Conclusion

In this paper, we proposed a new strategy for modeling thin composite structures by coupling a global plate model with a refined 3D model in zones of interest. The objective is to combine the efficiency of plate modeling for most of the structure with the added accuracy of 3D calculations where necessary (near edges, holes or defects). The hybrid plate/3D model is obtained after the convergence of a nonintrusive iterative process which enables one to use standard commercial software and industrial plate models without even having to represent the critical zones. The coupling is achieved by means of a preliminary analysis in which numerical Saint-Venant stress and warping fields which are suitable for the stacking sequence and the chosen discretization are recovered automatically. Several coupling methods are possible, but the warping-based technique is the simplest to implement and seems to lead to the least amount of perturbation at the interface. The convergence can be monitored by looking at the imbalance (in a plate sense) between the 3D domain and the plate domain. In the case of linear 3D domains, it can be accelerated by using a Krylov solver.

Our works in progress concern the application of the method to the simulation of complex assemblies of composite plates. In this case, the global model uses classical simplified 1D connectors to maintain computational efficiency while realistic 3D nonlinear models of bolts, rivets and damageable composites are used in the local models. A particular issue arises from the use of different in-plane discretizations at the interface between the local model and the global model. Two-scale approaches, such as that of [START_REF] Pa Guidault | A two-scale approach with homogenization for the computation of cracked structures[END_REF] in the context of the XFEM, have interesting advantages, such as the possibility to transfer average static plate quantities to the 3D model without introducing local forces and spurious edge effects. [5] Omar Bettinotti, Olivier Allix, and Benoît Malherbe. A coupling strategy for adaptive local refinement in space and time with a fixed global model in explicit dynamics. Computational Mechanics, pages 1-14, 2013.

[6] D. Caillerie. Non homogeneous plate theory and conduction in fibered composites. In E. Sanchez-Palencia and A. Zaoui, editors, Homogenization Techniques for Composite Media, pages 251-270. Springer, 1987.

Figure 1 :

 1 Figure 1: Cell problem, loading N x

Figure 2 :

 2 Figure 2: Cell problem, loading N y

Figure 3 :

 3 Figure 3: Cell problem, loading M x

Figure 4 :Figure 5 :

 45 Figure 4: Cell problem, loading M y

  Shear loading, zz component

Figure 7 :Figure 8 :

 78 Figure 7: Distribution of the nodal forces through the thickness for isotropic and orthotropic materials

Figure 9 :

 9 Figure 9: The extracted warping functions for isotropic and orthotropic materials

Figure 10 :Figure 11 :

 1011 Figure 10: The reference 3D model and the hybrid model

Figure 12 :

 12 Figure 12: The cantilever plate

Figure 13 :

 13 Figure 13: Evolution of stresses along the x-axis

Figure 14 : 3 Figure 15 : 1 Figure 16 :

 14315116 Figure 14: The generalized forces along the interface

Figure 17 :

 17 Figure 17: Stress and generalized force distributions near the hole

Table 1 :

 1 Construction of the Saint-Venant solutions

	•	N yy

Table 2 :

 2 The corrections between the first descent and the fifth iteration

		• 10 -2	-1.97
	error in Q x (%)	5.99	-35.94
	error in M xx (%)	1.18	3.3

for the sake of simplicity, we did not develop the right-hand-side in terms of plate components

Acknowledgement

This work was partially funded by the French National Research Agency as part of project ICARE (ANR-12-MONU-0002-04).