John Michopoulos

A Iliopoulos

B Staber

J Guilleminot

C Soize

Stochastic modeling and identification of an hyperelastic constitutive model for laminated composites

Keywords: Composite, identification, hyperelasticity, random field, stochastic modeling, uncertainty quantification

come

Stochastic modeling and identification of an hyperelastic constitutive model for laminated composites

Notation

The notations listed in Tabs. 1-4 will be used throughout this paper. B and B (r) Reference configuration of the composite and ply r x Generic point in the reference configuration ϕ Deformation map [F ϕ ] Deformation gradient associated with ϕ, [F ϕ (x)] = ∇ϕ(x) [F ] Any deformation gradient in M 3 + Table 3: Notation related to deterministic constitutive equations.

w and w (r) Strain energy function in the composite and ply r [P ]

First Piola-Kirchhoff stress tensor (f (r) ,

f (r) ⊥ )
Local basis of ply r µ (r) i , 1 i 5 Material parameter defining w (r) [M (r) ]

Structural tensor in ply r

Introduction

Characterization of the constitutive behavior of laminated polymer matrix composites (PMCs) via the solution of deterministic inverse problems, has been pioneered and pursued by the U.S. Naval Research Laboratory (US-NRL) for at least the last five decades [START_REF] Mast | Characterization of strain-induced damage in composites based on the dissipated energy density part iii. general material constitutive relation[END_REF][START_REF] Michopoulos | Towards the robotic characterization of the constitutive response of composite materials[END_REF][START_REF] Michopoulos | Data-driven design optimization for composite material characterization[END_REF][START_REF] Michopoulos | Advances on the constitutive characterization of composites via multiaxial robotic testing and design optimization[END_REF]. These approaches utilize dense data sets generated by custom made multi-degree of freedom robotic testing machines capable of applying multiaxial loading in order to activate rich sets of strain states in the tested materials. In addition, these approaches utilized formalisms of the recoverable and irrecoverable parts of the strain-energy density to construct proper objective functions in conjunction with the multiaxial data for determining the material properties of the bulk lamina participating into these laminates by minimizing these objective functions via optimization methodologies. In the context of the infinitesimal strains approximation these properties (including the elastic moduli) essentially appear as coefficients of monomials of strain components that make up the respective strain energy density functionals [START_REF] Michopoulos | Data-driven design optimization for composite material characterization[END_REF][START_REF] Michopoulos | Advances on the constitutive characterization of composites via multiaxial robotic testing and design optimization[END_REF]. An operationally useful assumption utilized in these studies was that the elastic properties of each lamina were assumed to remain constant throughout the domain of the plies consisting the laminates. However, the mechanical response of actual composites is very prone to uncertainties, both in the linear and the nonlinear regime and prior to final failure. This variability is essentially induced by complex material processing conditions during manufacturing (generating, for instance, topological imperfections usually in the form of local defects, such as voids, fiber bundling, fiber-matrix debonds, resin rich pockets etc.) and strongly impacts the performance of structural parts. One possible way to assess such uncertainties is to develop stochastic constitutive models that enrich physical models with a probabilistic description of all relevant uncertainties. While parametric models involving a few random elastic parameters have been used for many years, information-theoretic modeling approaches in which the elasticity tensor can be treated in full as a random matrix were proposed only recently. Such an approach was first introduced in [START_REF] Soize | Non-gaussian positive-definite matrix-valued random fields for elliptic stochastic partial derivative operators[END_REF] for triclinic materials, and was latter extended to all symmetry classes in [START_REF] Guilleminot | Stochastic model and generator for random fields with symmetry properties: Application to the mesoscopic modeling of elastic random media[END_REF][START_REF] Staber | Stochastic modeling and generation of random fields of elasticity tensors: a unified information-theoretic approach[END_REF]. Based on this formalism, statistical dependences between elastic moduli were also derived in [START_REF] Guilleminot | On the statistical dependence for the components of random elasticity tensors exhibiting material symmetry properties[END_REF]. In contrast, the case of nonlinear constitutive laws has received much less attention to date. The proper randomization of Ogden's strain energy functions was first addressed in [START_REF] Staber | Stochastic modeling of a class of stored energy functions for incompressible hyperelastic materials with uncertainties[END_REF][START_REF] Staber | Stochastic modeling of the ogden class of stored energy functions for hyperelastic materials: the compressible case[END_REF] for homogeneous materials (see also [START_REF] Mihai | Stochastic isotropic hyperelastic materials: constitutive calibration and model selection[END_REF]), and the case of spatially dependent behaviors was treated in [START_REF] Staber | A random field model for anisotropic strain energy functions and its application for uncertainty quantification in vascular mechanics[END_REF]. Studies reporting on the identification of such models can be found in, e.g., [START_REF] Staber | Stochastic hyperelastic constitutive laws and identification procedure for soft biological tissues with intrinsic variability[END_REF][START_REF] Brewick | Uncertainty quantification for constitutive model calibration of brain tissue[END_REF] for homogeneous biological tissues. In this context, the identification of stochastic models for anisotropic, heterogeneous nonlinear materials remain widely unexplored.

The present work is a first step towards the construction of an uncertainty quantification (UQ) framework that will ultimately support simple and robust decision making with respect to the design, certification and qualification of composite materials (and structures at large) under the presence of uncertainty. As the initial step in the development of this approach, the main goal of this paper is to introduce the formulation of a probabilistic modeling approach that accounts for the variability exhibited by the hyperelastic strain energy W and W (r) Random strain energy function in the composite and ply r W Regularized random strain energy function in the composite

{{G (r) i (x), x ∈ B (r) }} 5 i=1
Random fields of material parameters defining W (r) [M (r) ] Random structural tensor in ply r {{Ξ Correlation lengths along f (r) and f

(r) ⊥ , defining R δ
Coefficient of variation for the random material parameters function defining the (recoverable) mechanical response of a laminated composite. Two specific contributions are presented. First, an appropriate formulation defining the random field of anisotropic strain energy function is proposed. This stochastic model is consistent with the anisotropy exhibited by the composite, as well as with all mathematical requirements arising from the theory of nonlinear elasticity. Second, we address the identification of model parameters using experimental, three-dimensional full-field data acquired at the US-NRL [START_REF] Michopoulos | Data-driven design optimization for composite material characterization[END_REF][START_REF] Michopoulos | Overview of constitutive response characterization for composite materials via data-driven design optimization[END_REF][START_REF] Michopoulos | Advances on the constitutive characterization of composites via multiaxial robotic testing and design optimization[END_REF] by using custom developed non-contact strain-field measurement methods [START_REF] Iliopoulos | Direct strain imaging for full field measurements[END_REF][START_REF] Iliopoulos | Direct strain tensor approximation for full-field strain measurement methods[END_REF][START_REF] Iliopoulos | Meshless methods for full field displacement and strain measurement[END_REF]. A two-step strategy is followed where the mean model is calibrated through a linearization of the nonlinear model, while the parameters related to statistical fluctuations are identified by combining a principal component analysis with the maximum likelihood method. The rest of this paper is organized as follows. The modeling framework is presented in Section 3. Selected theoretical results related to the modeling of anisotropic stored energy functions are recalled, with the aim to precisely identify constraints on material parameters, in accordance with existence theorems in nonlinear elasticity. The probabilistic model is then introduced for a prototypical hyperelastic constitutive law, relevant to transversely isotropic materials. Forward simulations are also conducted to briefly illustrate the effect of model parameters on sample paths. Section 4 is focused on the identification of the non-Gaussian random field model. The experimental database, associated with a carbon-epoxy laminated composite, is first introduced. The calibration methodology is next detailed and finally deployed for the identification of the anisotropic strain energy function random field.

Stochastic modeling of the hyperelastic strain energy function

Background on hyperelastic constitutive models for transversely isotropic heterogeneous materials

In this section, we recall the kinematics and constitutive equations of deterministic hyperelasticity (see, e.g., [START_REF] Ciarlet | I: Three-Dimensional Elasticity[END_REF][START_REF] Ogden | Non-Linear Elastic Deformations[END_REF]). Since the laminate under consideration is made of unidirectional transversely isotropic plies, the presentation is restricted to the case of materials exhibiting such material symmetry.

Let (f (r) , f

⊥ ) be the basis of unit vectors defining the local frame of ply r (see Fig. 1) and let (e 1 , e 2 , e 3 ) be the canonical basis of R 3 . The plane of the plies is normal to e 3 . The local basis vectors are given by

f (r) = sin(θ)e 1 + cos(θ)e 2 , f (r)
⊥ = -cos(θ)e 1 + sin(θ)e 2 (1)
for odd plies, and read as

f (r) = -sin(θ)e 1 + cos(θ)e 2 , f (r)
⊥ = -cos(θ)e 1 -sin(θ)e 2 (2)
for even plies, with θ ∈ [0, π/2]. In this paper, the composite is defined by a stacking sequence of [+θ, -θ] 16 (meaning that the first ply is defined by an angle +θ, the second one by an angle -θ, and this two-ply structure is stacked up 16 times; hence a total of 32 plies); see Section 4.1 for a more complete description of the material. Let B be the reference configuration occupied by the composite and let ∂B denote its boundary. Each ply r has a reference configuration denoted by B (r) , such that B = ∪ 32 r=1 B (r) with B (r) ∩ B (s) = {∅} for r = s, 1 r, s 32. In the present case, all domains occupied by the plies are assumed to exhibit the same geometrical characteristics after material processing. The material undergoes a deformation map ϕ : B → R 3 that is assumed to be injective and orientation-preserving. The associated deformation gradient is denoted by

[F ϕ (x)] = ∇ϕ(x) , (3)
such that det([F ϕ (x)]) > 0 for all x in B. In the following, any deformation gradient evaluated at a fixed location x in B will be denoted by [F ] ∈ M 3 + . The material is assumed to be hyperelastic, meaning that there exists a stored energy function w : B × M 3 + → R such that the first Piola-Kirchhoff stress tensor can be obtained by

[P (x, [F ])] ij = ∂w(x, [F ]) ∂[F ] ij , 1 i, j 3 , (4)
for all [F ] ∈ M 3 + . The stored energy function w is further assumed to be frame-invariant, that is, it satisfies

w(x, [F ]) = w(x, [Q][F ]) , ∀[Q] ∈ SO(3) (5)
for all x in B. The stored energy function w can thus be decomposed as

w(x, [F ]) = 32 r=1 χ (r) (x)w (r) (x, [F ]) , ∀[F ] ∈ M 3 + , (6)
where x → χ (r) (x) is the characteristic function of the r-th ply (that is, r) , and χ (r) (x) = 0 otherwise) and w (r) : B (r) × M 3 + → R is the stored energy function in B (r) . In order to represent the material symmetry, the stored energy function w (r) must satisfy the following invariance property:

χ (r) (x) = 1 if x ∈ B (
w (r) (x, [F ]) = w (r) (x, [F ][Q]) , ∀[Q] ∈ G (r) , (7)
for all x in B (r) , where G (r) ⊂ SO(3) is the symmetry group associated with the preferred direction f (r) , i.e.,

G (r) = {[Q] ∈ SO(3) : [Q]f (r) = f (r) } . (8)
In this work, and without loss of generality, the anisotropic stored energy function proposed in [START_REF] Schröder | Anisotropic polyconvex energies on the basis of crystallographic motivated structural tensors[END_REF] (for all symmetry classes) is selected to model each ply (see Appendix A). The stored energy function w (r) is specifically decomposed as

w (r) (x, [F ]) = w (r) iso (x, [F ]) + w (r) ani (x, [F ]) -ψ * (x) , (9)
in which the isotropic and anisotropic contributions, as well as the additional term ψ * , are defined as follows.

(I.) The isotropic potential w (r) iso corresponds to a compressible Mooney-Rivlin model and is given by

w (r) iso (x, [F ]) = µ (r) 1 (x) [F ] 2 F + µ (r) 2 (x) Cof([F ]) 2 F + µ (r) 3 (x) det([F ]) 2 -d r (x) log(det([F ])) , (10)
for all [F ] ∈ M3 + , where {µ i (x)} 3 i=1 are model parameters, with µ 3 (x) > 0, and d r (x) = 2µ

(r) r) . The fields x → µ (r) 1 (x) and x → µ (r) 2 (x) are assumed to be uniformly bounded from below by strictly positive constants, i.e., there exists (τ

1 (x) + 4µ (r) 2 (x) + 2µ (r) 3 (x) ∀x ∈ B (
(r) 1 , τ (r) 2 ) such that ∀x ∈ B (r) , µ (r) 1 (x) τ (r) 1 > 0 , µ (r) 2 (x) τ (r) 2 > 0 . (11)
(II.) The anisotropic contribution w (r)

ani is given by

w (r) ani (x, [F ]) = Tr([F ] T [F ][M (r) (x)]) 1+β4 (1 + β 4 )Tr([M (r) (x)]) β4 + Tr(Cof([F ] T [F ])[M (r) (x)]) 1+β5 (1 + β 5 )Tr([M (r) (x)]) β5 -2Tr([M (r) (x)]) det([F ]) , (12)
where β 4 0, β 5 0, and [M (r) (x)] is a symmetric and positive-definite second-order tensor that is referred to as to a structural tensor at point x. This tensor can be decomposed as

[M (r) (x)] = µ (r) 4 (x)f (r) ⊗ f (r) + µ (r) 5 (x)([I 3 ] -f (r) ⊗ f (r) ) , ∀x ∈ B (r) , (13)
in which x → µ r) , is implied by the positive-definiteness of [M (r) (x)]). It satisfies, by construction, the invariance property

(r) 5 (x) > 0 for all x ∈ B (
[M (r) (x)] = [Q][M (r) (x)][Q] T , ∀[Q] ∈ G (r) , (14)
where G (r) is given by Eq. ( 8).

(III.) The term ψ * is defined as

ψ * (x) = µ (r) 3 (x)δ (r) 0 (x) 2 -d r (x) log(δ (r) 0 (x)) -2Tr([M (r) (x)])δ (r) 0 (x) , (15)
where the mapping x → δ (r) 0 (x) is defined as

δ (r) 0 (x) = 2Tr([M (r) (x)]) + (Tr([M (r) (x)])) 2 + 2µ (r) 3 (x)d r (x) 2µ (r) 3 (x) . (16)
In this setting, it can be shown that the stored energy function w (r) satisfies the objectivity constraint and all the mathematical requirements involved in the existence theorems of nonlinear elasticity [START_REF] Ball | Convexity conditions and existence theorems in nonlinear elasticity[END_REF][START_REF] Müller | On a new class of elastic deformations not allowing for cavitation[END_REF] (see also [START_REF] Ciarlet | I: Three-Dimensional Elasticity[END_REF]); see Appendix B. Remark 3.1. Using Eq. ( 13), the anisotropic invariants involved in the right-hand side of Eq. ( 12) can be expressed as

Tr([F ] T [F ][M (r) (x)]) = µ (r) 4 (x) [F ]f (r) 2 + µ (r) 5 (x)( [F ]f (r) ⊥ 2 + [F ]e 3 2 ) ( 17
)
and

Tr(Cof([F ] T [F ])[M (r) (x)]) = µ (r) 4 (x) Cof([F ])f (r) 2 + µ (r) 5 (x)( Cof([F ])f (r) ⊥

Construction of the random field model

In this section, the proper randomization of the spatially dependent, anisotropic stored energy function introduced in Section 3.1 is addressed, following the methodology proposed in [START_REF] Staber | A random field model for anisotropic strain energy functions and its application for uncertainty quantification in vascular mechanics[END_REF] for hyperelastic materials. The approach consists of substituting the fields of deterministic material parameters by non-Gaussian random fields, the probabilistic models of which must be constructed. In this work, the exponents β 4 and β 5 are taken deterministic and could be randomized as well, depending on the application.

Let {W (r) (x, •), x ∈ B (r) } be the stored energy function random field, defined on a probability space (Ω, T , P) and corresponding to the probabilistic modeling of the anisotropic stored energy function w

(r) , 1 r 32. Let {G (r) (x) = (G (r) 1 (x), . . . , G (r)
5 (x)), x ∈ B (r) } be the random field corresponding to the probabilistic modeling of the vector-valued model parameter field x → (µ

(r) 1 (x), . . . , µ (r) 5 (x)), associated with {W (r) (x, •), x ∈ B (r) }.
Some mathematical preliminaries regarding the functional form of the stochastic potential W (r) (x, •), x being fixed in B (r) , are first discussed in Section 3.2.1. The construction of the probabilistic model for {G (r) (x), x ∈ B (r) } is then tackled in Section 3.2.2. Note at this point that this model can be used per se, or can be integrated as an optimal prior model in a Bayesian approach (see Chapter 26 in [START_REF]Handbook of Uncertainty Quantification[END_REF] for methodological aspects). In the following, it is assumed that {G (r) (x), x ∈ B (r) } is a second-order random field, and that the mean function of x → G (r) (x) is independent of x and r, with

E{G (r) (x)} = g , (19)
where E denotes the mathematical expectation and g = (g 1 , . . . , g 5 ) is such that g j > 0 for 1 j 5.

Regularized form of the stochastic stored energy function

In order to ensure uniform boundedness from below, a regularized form W (r) of the random stored energy function is introduced as

W (r) (x, [F ]) = 1 1 + (W (r) 1 (x, [F ]) + w 1 ([F ])) + W (r) 2 (x, [F ]) + W (r) vol (x, det([F ])) -W (r) vol (x, D (r) 0 ) , ( 20
)
where 0 < 1 is an arbitrary small parameter. The stochastic potential W (r)

1 (x, •) corresponds to the stochastic modeling of the first two terms in the isotropic contribution given by Eq. ( 10) and writes

W (r) 1 (x, [F ]) = G (r) 1 (x) [F ] 2 F + G (r) 2 (x) Cof([F ]) 2 F , ∀[F ] ∈ M 3 + . (21)
By construction, the mean function x → E{W (r) 1 (x, •)} is independent of x and is denoted by

E{W (r) 1 (x, •)} = w 1 (•), with w 1 ([F ]) = g 1 [F ] 2 F + g 2 Cof([F ]) 2 F . (22)
By substituting Eq. ( 21) and Eq. ( 22) into Eq. ( 20), a regularized, isotropic strain energy function W 1, can be introduced as

W 1, (x, [F ]) = G 1 (x) [F ] 2 F + G 2 (x) Cof([F ]) 2 F , (23)
with

G 1 (x) = 1 1 + (G 1 (x) + g 1 ) (24)
and

G 2 (x) = 1 1 + (G 2 (x) + g 2 ) . (25)
Following Eq. ( 12), the random contribution

W (r)
2 takes the form

W (r) 2 (x, [F ]) = Tr([F ] T [F ][M (r) (x)]) 1+β4 (1 + β 4 )(Tr([M (r) (x)])) β4 + Tr(Cof([F ] T [F ])[M (r) (x)]) 1+β5 (1 + β 5 )(Tr([M (r) (x)])) β5 , (26)
where [M (r) (x)] is now a random matrix given by

[M (r) (x)] = G (r) 4 (x)f (r) ⊗ f (r) + G (r) 5 (x)([I 3 ] -f (r) ⊗ f (r) ) . (27)
Similarly, from Eq. ( 15), the random volumetric contribution

W (r)
vol is given by

W (r) vol (x, δ) = G (r) 3 (x)δ 2 -D r (x) log(δ) -2Tr([M (r) (x)])δ , (28)
where the random variable D r (x) reads as D r (x) = 2G

(r)

1 (x) + 4G (r) 2 (x) + 2G (r)
3 (x). Finally, the expression of the random variable D (r) 0 involved in the very last term of the right-hand side of Eq. ( 20) can be deduced from Eq. ( 16), i.e.,

D (r) 0 (x) = 2Tr([M (r) (x)]) + (Tr([M (r) (x)])) 2 + 2G (r) 3 (x)D r (x) 2G (r) 3 (x) . ( 29
)
In order to ensure the well-posedness of the stochastic boundary value problem, it is further assumed that the random variables

G (r) 1 (x), . . . , G (r)
5 (x) are such that

G (r) j (x) > 0 , 1 j 5 , (30)
for all x in B (r) , almost surely. By substituting Eqs. ( 21)-( 30) into the stochastic stored energy function given by Eq. ( 20), it can be deduced that

W (r) (x, [F ]) c w 1 ([F ]) , (31)
with c = /(1 + ). Eq. [START_REF] Iliopoulos | Loading subspace selection for multidimensional characterization tests via computational experiments[END_REF] shows that w 1 is a deterministic lower bound for the random strain energy function, which then satisfies uniform growth conditions (see Appendix B).

Definition of the probabilistic model for the set of coefficients

Given the inequality constraints raised by the well-posedness of the stochastic boundary value problem, the random fields {{G (r) i (x), x ∈ B (r) }} 5 i=1 cannot be Gaussian and must be modeled as non-Gaussian random fields. Following the methodology proposed in [START_REF] Soize | Non-gaussian positive-definite matrix-valued random fields for elliptic stochastic partial derivative operators[END_REF] for anisotropic linear elastic materials (see also [START_REF] Guilleminot | Stochastic model and generator for random fields with symmetry properties: Application to the mesoscopic modeling of elastic random media[END_REF][START_REF] Staber | Stochastic modeling and generation of random fields of elasticity tensors: a unified information-theoretic approach[END_REF]) and advanced in [START_REF] Staber | A random field model for anisotropic strain energy functions and its application for uncertainty quantification in vascular mechanics[END_REF] for spatially dependent anisotropic stored energy functions, each random field {G (r) i (x), x ∈ B (r) } is then expressed through a nonlinear transformation of a Gaussian random field, that is:

G (r) i (x) = H i Ξ (r) i (x) , ∀x ∈ B (r) , 1 i 5 , (32)
where the mapping H i is constructed below and such that all considered constraints (e.g., Eq. ( 30)) hold, and

{Ξ (r) i (x), x ∈ B (r)
} denotes the restriction to B (r) of a real-valued, homogeneous and normalized Gaussian random field. Furthermore, and since the constraints on material parameters do not contain cross information between them (see Section 3.1), the information-theoretical framework (see below) implies that

{G (r) i (x), x ∈ B (r) } and {G (r) j (x), x ∈ B (r) } (and similarly, {Ξ (r) i (x), x ∈ B (r) } and {Ξ (r) j (x), x ∈ B (r) }), 1 i = j
5, are statistically independent random fields. In order to make the identification of the probabilistic model tractable, it is assumed below that (i ) the covariance functions defining the underlying Gaussian random fields exhibit the same algebraic form, regardless of the material parameter; (ii ) the hyperparameters defining the random fields do not change from one ply to another. Note that this latter assumption is consistent with the very structure of the material (which is manufactured by stacking up independent samples of unidirectional plies), provided that the impact of the curing process (which may generate spatial fluctuations in the mean fields for thick composite parts, for instance) can reasonably be neglected.

Let η → R(η) denote the selected covariance function of the Gaussian fields. For a given ply r, this function is given by R

(η) = ρ 1 ( η, f (r) ) × ρ 2 ( η, f (r) ⊥ ) , (33)
where •, • denotes the Euclidean inner product and the univariate covariance functions are defined as

ρ k (d k ) = 2L k πd k 2 sin 2 πd k 2L k , 1 k 2 . ( 34
)
In the equation above, the parameters L 1 and L 2 correspond to the correlation lengths along the directions f (r) and f

⊥ , respectively. The definition of the mappings H i , 1 i 5, is then achieved by invoking the maximum entropy principle [START_REF] Shannon | A mathematical theory of communication[END_REF][START_REF] Jaynes | Information theory and statistical mechanics i[END_REF][START_REF] Jaynes | Information theory and stastitical mechanics ii[END_REF] to construct the probability distributions of G (r) i (x), with 1 i 5 and x fixed in B (r) [START_REF] Staber | A random field model for anisotropic strain energy functions and its application for uncertainty quantification in vascular mechanics[END_REF]; see [START_REF] Staber | Stochastic modeling of a class of stored energy functions for incompressible hyperelastic materials with uncertainties[END_REF][START_REF] Staber | Stochastic modeling of the ogden class of stored energy functions for hyperelastic materials: the compressible case[END_REF] for the case of homogeneous hyperelastic materials, as well as [START_REF] Staber | Stochastic hyperelastic constitutive laws and identification procedure for soft biological tissues with intrinsic variability[END_REF] for a discussion regarding model identification. The knowledge of these probability distributions will be used, below, to derive analytical expressions for the transformations {H i } 5

i=1 . From a probabilistic standpoint, this approach completely defines the full systems of marginal distributions of these fields. Under the constraints

E{G (r) i (x)} = g i (35)
and E{log(G

(r) i (x))} = χ (r) i , |χ (r) i | < +∞ , (36)
arising from Eq. ( 19) for the former and from the fact that G

i (x) > 0 for the latter (see [START_REF] Soize | A nonparametric model of random uncertainties for reduced matrix models in structural dynamics[END_REF]), it can be shown that the random variable G (r) i (x) (x being fixed) follows a Gamma distribution. Consequently, the random field {G (r) i (x), x ∈ B (r) } can simply be defined, in the present case, as

G (r) i (x) = F -1 Γ(ki,θi) Φ(Ξ (r) i (x)) , 1 i 5 , (37)
where F -1 Γ(ki,θi) is the inverse cumulative distribution function of the Gamma law with parameters (k i , θ i ) and Φ is the cumulative distribution function of the normal distribution. Eq. (37) explicitly defines the mapping H i , introduced in Eq. [START_REF] Dacorogna | Direct methods in the calculus of variations[END_REF], and upon expressing the parameters of the Gamma distribution in terms of second-order properties, it can be written as

G (r) i (x) = F -1 Γ(δ -2 ,δ 2 g i ) Φ(Ξ (r) i (x)) , 1 i 5 , (38)
in which it was assumed that all random fields exhibit the same level of statistical fluctuations, measured by the coefficient of variation δ.

As pointed out in [START_REF] Staber | A random field model for anisotropic strain energy functions and its application for uncertainty quantification in vascular mechanics[END_REF], it should be noticed that the above stochastic modeling step encompasses the case of isotropic materials and can readily be generalized to other types of strain energy functions.

Simulations of the random field model

In this section, forward simulations are performed in order to qualitatively illustrate the effect of model parameters on sample characteristics. The random field is generated on the geometry depicted in Fig. 2, which corresponds to the configuration that will be considered in the identification of the random field model in Section 4. Here, a single-ply laminate is considered, and simulations are conducted for the first material parameter random field {G 1 (x), x ∈ B (1) } is carried out using a spectral expansion [START_REF] Shinozuka | Simulation of stochastic processes by spectral representation[END_REF][START_REF] Poirion | Numerical methods and mathematical aspects for simulation of homogeneous and non homogeneous gaussian vector fields[END_REF]. Realizations of the Gaussian field for different sets of correlation lengths (L 1 , L 2 ) are shown in Figs. 34.

It can be seen that the samples exhibit patterns that are characteristic of an isotropic correlation structure for L 1 = L 2 . In contrast, the realizations of the Gaussian field show pronounced anisotropy along either f ⊥ or f , depending on the ordering of the spatial correlation lengths L 1 and L 2 . Realizations of the non-Gaussian field {G

Description of the experimental setup and model identification

In this section, the identification of the random field model for the anisotropic strain energy function is pursued. The experimental database is first described in detail in Section 4.1. The proposed methodology enabling the identification of the hyperparameters is then discussed and applied in Section 4.2.

Description of the experimental database

In this application, the unidirectional ply corresponds to a 3501-6 epoxy resin matrix, reinforced with AS4 carbon fibers. Four ply angles θ were considered and two specimens were tested in tension for each angle (hence, a total of 8 specimens tested). As indicated previously, the layup stacking sequence of each specimen was [θ, -θ] 16 . A displacement distribution u d was applied on the top and fixed boundary conditions were prescribed at the bottom of the composite samples; see Fig. 7. This loading condition represents one of 72 loading paths defined to homogeneously sample a 4-dimensional loading space spanned by the tension and the 3-rotation loading bases. This loading path is linear and starts from the origin (0,0,0,0) and corresponds to the case of pure tension. The particular loading subspace where this path is defined is one of the 15 4-dimensional subspaces that can be selected out of the full 6-dimensional loading space that NRL66.3 can apply loads in and it was selected by a process described in [START_REF] Iliopoulos | Loading subspace selection for multidimensional characterization tests via computational experiments[END_REF]. The loading paths sampling this 4-dimensional subspace are linear combinations of the loading space bases and are applied (one loading path per specimen) on the upper edge of the specimen via the US-NRL's NRL66.3 6-DoF robotic loading frame as described in [START_REF] Michopoulos | Data-driven design optimization for composite material characterization[END_REF][START_REF] Michopoulos | Advances on the constitutive characterization of composites via multiaxial robotic testing and design optimization[END_REF]. Note that the geometry presents two notches that were specifically introduced to disturb the strain field and ensure that some areas of the material will exhibit material nonlinearity in their response [START_REF] Michopoulos | Advances on the constitutive characterization of composites via multiaxial robotic testing and design optimization[END_REF]. Snapshots of the displacement and strain fields were acquired using the digital imaging based method described in [START_REF] Iliopoulos | Meshless methods for full field displacement and strain measurement[END_REF] on the front and rear faces of the composite coupons, until failure was reached (see [START_REF] Michopoulos | Advances on the constitutive characterization of composites via multiaxial robotic testing and design optimization[END_REF] for a complete description of the experimental setup). The measurement was performed on the front and rear faces motivated by the need to have access to the largest data set possible from each test performed. It should be noted that marking the specimens from both sides at the same x, y is particularly challenging and some misalignment of the optical metrology patterns is always expected. For this reason, the present analysis does not make any assumption about alignment or not of the front and back patterns. Instead, it allows for calculating the displacement and strain fields through vector and tensor interpolation respectively, at the locations required for further analysis. A summary of the experimental parameters is given in Tab. 5, where n exp denotes the number of observation points on the back and front faces and T exp is the number of loading steps until the maximum displacement max s∈Texp u d (s) imposed on the top of the specimen has been reached. For further developments related to the statistical inverse problem, the experimental database was postprocessed as follows. The set of observation points with the lower number of entries (n exp = 2357) was selected as the reference set. For every observation point in this set, the nearest neighbors belonging to the seven remaining sets were then determined using the k-nearest neighbors algorithm. An additional set of observation points was then obtained by solving an optimization problem such that this new set constitutes an accurate approximation of all the sets of observation points for all the experimental realizations. The resulting set contains n exp = 2355 observation points.

Similarly, the sets of loading steps were approximated by a single set T exp such that the maximum prescribed displacement on the top is given by 0.1mm with 10 equally spaced steps, regardless of the ply angle and realization. In this setting, the failure of the composite material is not taken into account and may be addressed in the future. The experimentally obtained boundary conditions in terms of the applied displacement evolution and the corresponding reaction forces evolution for all four values of the lamination angle are presented in Fig. 8 while their associated approximations are shown in Fig. 9.

The available realizations of the component x → [E exp (s * ; x, θ)] 22 of the strain field are shown in Figs. 10111213, where [E exp (s; x, θ)] is the experimental Green-Lagrange deformation tensor at loading step s and point x, for a composite defined by a stacking sequence of [+θ, -θ] 16 , and s * is the final loading step such that the prescribed displacement is u d = 0.1 mm.

The variability in the reconstructed strain field is clearly observed and will be exploited in the next section of identify the parameters of the random field model.

Identification methodology

This section is focused on the identification of the random field model for the anisotropic strain energy function, presented in Section 3.2. This model depends on the mean vector g, the deterministic exponents (β 4 , β 5 ), the coefficient of variation δ, and the spatial correlation lengths (L 1 , L 2 ). In this work, the identification problem is attacked through a two-step strategy. The mean model and the deterministic exponents are first determined in Section 4.2.1, by imposing a fit of the linearized elasticity tensor. The hyperparameter controlling the statistical fluctuations and the correlation lengths are then identified in Section 4.2.2, using the maximum likelihood method.

It is assumed below that the level of experimental noise is negligible as compared to that of intrinsic material variability, so that it can be neglected in the calibration stages. In addition, we note that this noise should also be modeled and identified, and that its consideration is more relevant in a Bayesian approach to model calibration.

Identification of the mean model and deterministic exponents

In order to identify the mean vector g and the exponents (β 4 , β 5 ), the mean elasticity tensor at small strains is fitted with known values available in [START_REF] Michopoulos | Overview of constitutive response characterization for composite materials via data-driven design optimization[END_REF][START_REF] Michopoulos | Data-driven design optimization for composite material characterization[END_REF]. From the definition of the stochastic stored energy function, it can be deduced that the associated random elasticity tensor at small strains reads as

A (r) (x) = 1 1 + A (r) iso (x) + E{ A (r) iso (x) } + A (r) sym (x) , (39)
where the isotropic contribution writes and the anisotropic term is given by

A (r) iso (x) = (4G (r) 2 (x) + 4G (r) 3 (x))[I 3 ] ⊗ [I 3 ] + (4G (r) 1 (x) + 4G (r) 2 (x))[I 3 ] [I 3 ] ( 40
)
A (r) sym (x) = 4(β 4 + β 5 ) [M (r) (x)] ⊗ [M (r) (x)] Tr([M (r) (x)]) + 4Tr([M (r) (x)])(1 + β 3 + β 5 )[I 3 ] ⊗ [I 3 ] -4(1 + β 5 )([I 3 ] ⊗ [M (r) (x)] + [M (r) (x)] ⊗ [I 3 ]) + 4([I 3 ] [M (r) (x)] + [M (r) (x)] [I 3 ]) . (41)
Let h = (g, β 4 , β 5 ) be the sought-for vector of mean values, and denote by A(h) the nominal value of the non-regularized linearized elasticity tensor, obtained by summing up the right-hand sides in Eqs. ( 40)-(41):

A(h) = (4g 2 + 4g 3 )[I 3 ] ⊗ [I 3 ] + (4g 1 + 4g 2 )[I 3 ] [I 3 ] + 4(β 4 + β 5 ) [M (r) ] ⊗ [M (r) ] Tr([M (r) ]) + 4Tr([M (r) ])(1 + β 3 + β 5 )[I 3 ] ⊗ [I 3 ] -4(1 + β 5 )([I 3 ] ⊗ [M (r) ] + [M (r) ] ⊗ [I 3 ]) + 4([I 3 ] [M (r) ] + [M (r) ] [I 3 ]) . (42)
In the equation above, the mean value [M (r) ] of the structural tensor [M (r) (x)] is given by (see Eq. ( 27))

[M (r) ] = g 4 f (r) ⊗ f (r) + g 5 ([I 3 ] -f (r) ⊗ f (r) ) . (43)
For given values of the mean model and deterministic exponents, Eq. (42) allows for computing the mean engineering constants, namely the longitudinal and transversal moduli (E L , E T , G L ), and Poisson's ratios (V L , V T ). The vector of mean parameters and deterministic exponents is then defined as: where z(q) = (E L (q), E T (q), V L (q), V T (q), G L (q)) is the vector gathering the mean engineering constants associated with the tensor A(q) , obtained by substituting q for h in Eq. ( 42), and z is the vector of target mean values defined in [START_REF] Michopoulos | Overview of constitutive response characterization for composite materials via data-driven design optimization[END_REF][START_REF] Michopoulos | Data-driven design optimization for composite material characterization[END_REF] as z = (139.1400, 10.3620, 0.2860, 0.3080, 6.6478) (units for elastic moduli: GPa). Following the constraints on material parameters introduced in Section 3, the admissible set H h is defined as H h = (0, +∞) 5 × [0, +∞) 2 . The optimization problem defined by Eq. ( 44) is solved with Matlab function fmincon, for a series of 100 independent guesses obtained from Latin Hypercube Sampling (LHS). The final values of the cost function with respect to the index of the (vector-valued) starting point are shown in the left panel of Fig. 14, in which filled bullets are used to highlight starting points for which the final value of the cost function is below 10 -10 . The convergence of the optimization algorithm is shown in the right panel of Fig. 14 for the optimal starting point. Similarly, the identified values of the mean parameters g 1 , . . . , g 5 and exponents β 4 , β 5 are shown for all starting points in Figs. [START_REF] Michopoulos | Overview of constitutive response characterization for composite materials via data-driven design optimization[END_REF], 16 and 17. The filled bullets represent the cases for which the cost function is below 10 -10 , and it can be seen that the identification turns out to be robust with respect to the initial guess of the solution. Note that the point associated with sample index 57 (see the right panel in Fig. 15, as well as the left panel in Fig. 16 and the right panel in Fig. 17) is not a robust candidate, since the corresponding final value of the cost function is the largest among all other combinations (see Fig. 14,left).

h = min q ∈ H h 5 i=1 (z i (q) -z i ) 2 z 2 i , (44)
The mean parameters are finally identified as follows: Figure 15: Identified optimal values of the mean parameters g 1 and g 2 , obtained for the 100 starting points. Filled bullets correspond to cases for which the final value of the cost function is below 10 -10 . In the case of g 2 , the limits of the y-axis have been reduced for readability.

Identification of the fluctuation parameter and correlation lengths

Let ξ = (L 1 , L 2 , δ) be the vector-valued parameter to be identified, gathering the fluctuation parameter and the correlation length, and denote by H ξ = (0, +∞) 2 × (0, 1) its admissible set. For any angle θ (which specifies the laminate configuration), ζ ∈ H ξ and the value of h identified in Section 4.2.1, we introduce the R ñexp -valued random variable Y(ζ;θ) (recall that ñexp is the number of observation points) defined component-wise as

Y i (ζ; θ) = log   s ∈ T exp [E comp (s; x (i) , ζ, θ)] 2 F   , 1 i ñexp , (47)
where [E comp (s; x (i) , ζ, θ)] is the Green-Lagrange deformation tensor associated with the computational model at loading step s ∈ T exp , taken at location x (i) on the observed boundaries. A few remarks regarding the definition of the random variable Y(ζ;θ) are relevant at this point. First, the summation in the righthand side of Eq. (47) corresponds, for a given component Y i (ζ; θ) (that is, for a given observed point), to the natural discretization of the L 2 norm along the loading paths. Second, it should be noticed that the logarithm is introduced to make the random variable more sensitive to small values of the norm. Finally, the introduction of a vector-valued variable allows spatial information to be preserved and will allow us to formulate the identification problem in reduced coordinates, using a statistical reduction (see below). Let Figure 16: Identified optimal values of the mean parameters g 3 , g 4 and g 5 (from left to right), obtained for the 100 starting points. Filled bullets correspond to cases for which the final value of the cost function is below 10 -10 . In the case of g 3 , the limits of the y-axis have been reduced for readability.

Y (m) = m Y + m k=1 λ k η k ϕ (k) , (48)
where {(λ k , ϕ (k) )} k are the sets of associated eigenvalues and eigenvectors of [C Y ], and η k is the real-valued random variable such that

η k = 1 √ λ k Y -m Y , ϕ (k) , (49)
The order m of the reduced model is deduced from a standard convergence analysis, looking at the convergence of the error function m → E (m) defined as

E (m) = 1 - m α=1 λ α ñexp α=1 λ α . ( 50
)
Note at this stage that m is implicitly depending on both θ and ζ. For a given fiber orientation θ, let Y exp (θ) be the vector-valued random variable with components

Y exp i (θ) = log   s ∈ T exp [E exp (s; x (i) , θ)] 2 F   , 1 i ñexp , (51)
where we recall that [E exp (s; x (i) , θ)] is the experimental version of the random Green-Lagrange deformation tensor at loading step s and location x (i) . Note that Y(ζ;θ) and Y exp (θ) are not defined on the same probability space. Let y exp (θ, ω exp ), ∈ {1, 2}, denote the experimental realizations of Y exp (θ). Consider next the reduced coordinates defined as the projection of the experimental samples onto the basis obtained from the computational model (for fixed values of the hyperparameters, that is, for a fixed value of ζ):

η exp k (ζ; θ, ω exp ) def = 1 λ k (ζ; θ) y exp (θ, ω exp ) -m Y (ζ; θ), ϕ (k) (ζ; θ) , 1 k m(θ, ζ) , (52)
where the set of eigenpairs {(λ k , ϕ (k) )} k was introduced below Eq. ( 48) and parameter dependence is now recalled. The vector-valued parameter ξ is then defined as

ξ = arg max ζ ∈ H ξ L(ζ) , (53)
where L : H ξ → R denotes the log-likelihood function given by

L(ζ) = 4 j=1 2 =1 m(θj ,ξ) k=1 log fη k (θj ,ζ) η exp k ζ; θ j , ω exp . ( 54
)
In Eq. ( 54), fη k (θj ,ζ) is the probability density function of random variable η k (θ j , ζ) (see Eq. ( 49)) estimated with a kernel density estimator, using Monte Carlo simulations (with 400 independent realizations) and the computational model parametrized by ζ (and the value of h identified in the previous section) for the configuration defined by θ j (i.e., with a stacking sequence [+θ j , -θ j ] 16 ). The definition of the angles {θ j } 4 j=1 is provided in Tab. 5. From a computational standpoint, the mesh of the composite domain contains 139 200 8-node hexahedral finite elements (395 721 degrees of freedom). The nonlinear problem is solved with an in-house finite element implementation, using a standard Newton-Raphson method with a load stepping scheme (used libraries: boost, ParMETIS, Trilinos). Linear systems are solved with the preconditioned conjugate gradient method.

Given the computational burden associated with this identification step, the non-convex optimization problem given by Eq. ( 53) is solved by evaluating the cost function on a grid with equally spaced points. Other strategies could be pursued to solve this problem and in particular, to decrease the computational cost associated with the propagation of uncertainties through nonlinear models. A detailed discussion about stochastic solvers is outside the scope of this work, and the reader is referred to the general review [START_REF]Handbook of Uncertainty Quantification[END_REF]. Here, the grid is specifically defined by imposing reasonable physics-guided constraints, in particular for the correlation lengths (which are, for instance, assumed to be smaller than the smallest characteristic length of the specimens). The graph of ζ → L(ζ) is shown in Fig. 18, and the optimal values are found as ξ = (11.18, 11.18, 0.1) .

(55)

The fact that the correlation lengths along f (r) (which is the direction of the fibers) and f

⊥ are found to be equal constitutes an unexpected result. However, this may be a consequence of the relative "isotropization" of the bulk response of the composite with respect to the loading conditions, achieved by the stacking sequence (recall that only boundary measurements are made). The level of statistical fluctuations is rather contained (δ = 0.1), and is consistent with typical values obtained from multiscale simulations with non-separated scales on similar materials (see, e.g., Section 4.2 in [START_REF] Staber | Stochastic modeling of the ogden class of stored energy functions for hyperelastic materials: the compressible case[END_REF]). Realizations of the random fields associated with the engineering constants for transversely isotropic materials are shown in Figs. [START_REF] Ciarlet | I: Three-Dimensional Elasticity[END_REF], 20 and 21, hence showing the substantial variability of the mechanical properties over the considered ply. It should be noticed that the fluctuation patterns and ranges in these figures are induced by the values of the hyperparameters identified using the database. In particular, the parameters L 1 and L 2 control the frequency of the spatial oscillations along the local principal directions in the ply, while the range of fluctuations is controlled by the parameter δ. Given the very limited number of experimental realizations, quality assessment is finally achieved by controlling that most experimental realizations for some components of the Green-Lagrange deformation tensor at observable points are included in the support of the associated probability density functions estimated with the identified stochastic model (note at this point that the experimental realizations of the deformation fields cannot be compared with numerical ones, since the probability spaces are different). The plots for the components 11 and 22 of the deformation tensor at a given point x obs , located near a notch points for the component 22 exhibit a non-vanishing probability level, whereas the realizations for component 11 fall outside the estimated support of the probability density function. A possible explanation for this result is the small magnitude of the uniaxial deformation along e 1 , which makes experimental data more sensitive to measurement noise. Moreover, it should be noticed that this constitutes a severe validation test for the proposed model and identification strategy, since (i) no statistical treatment (e.g., local averaging) is performed on the experimental data, and (ii) potential measurement noise is not accounted for. In this context, the use of a Bayesian approach would be a very natural extension that is left for future works.

Conclusions and Plans

In this work, we have addressed the construction and inverse identification of a random field model describing an anisotropic stored energy function relevant to transversely isotropic materials. The stochastic model relies in part on information theory, which is used to properly integrate constraints related to the well-posedness of the nonlinear boundary value problem. In order to account for potential data paucity at the identification stage, the formulation involves a limited number of parameters that control mean and dispersion values, as well as spatial correlation lengths along appropriate directions. The calibration strategy was tackled as a two-step problem where the linearized model is first considered to identify the mean values. In a second step, principal component analysis is combined with the maximum likelihood method to estimate the dispersion parameter and the correlation lengths. The framework was finally deployed to identify the hyperelastic constitutive model of a carbon-epoxy laminate. This research effort is a first step towards the construction of a UQ framework that will support the development of an approach providing simple and robust decision making with regard to design, certification and qualification of composite materials and structures under the presence of uncertainty. Two envisioned improvements of the framework are the inclusion of additional multiaxial loading paths and the extension of the approach to include the irrecoverable part of the strain energy density. Another potential extension of this work that would be applicable for the high performance high strength high stiffness PMCs, would be to utilize an infinitesimal strain approximation for the strain energy density and the typical Hooke's tensor as the container of the anisotropic elastic properties.

Table 1 :i=1 a i b i • 2 3 +

 123 Notation related to vectors and matrices. a, a and [A] Deterministic scalar, vector and matrix A, A and [A] Random scalar, vector and matrix •, • Euclidean inner product, a, b = n Euclidean norm, a 2 = a, a a ⊗ b Outer product of vectors a and b [I 3 ] Identity matrix of rank 3 Tr(•) Trace of a matrix, Tr([A]) = n i=1 A ii • F Frobenius norm of [A], [A] F = Tr([A] T [A]) M Set of (3 × 3) matrices with strictly positive determinants Preprint submitted to Computer Methods in Applied Mechanics and Engineering June 3, 2019

i

  (x), x ∈ B (r) }} 5 i=1 Normalized Gaussian random fields, associated with ply r η → R(η) Covariance function defining the Gaussian fields L 1 and L 2

Figure 1 :

 1 Figure 1: Definition of the global and local basis vectors.

  ) and x → µ (r) 5 (x) are fields of model parameters (note that µ (r) 4 (x) > 0 and µ

  ), x ∈ B(1) }, with g 1 = 1. The numerical simulation of the associated Gaussian field {Ξ

  ), x ∈ B(1) } are then shown in Figs.5-6, for different values of the correlation lengths and coefficient of variation δ. It can be observed that the field {G (r) j (x), x ∈ B (r) } inherits the aforementioned patterns associated with the underlying Gaussian field.

Figure 2 :

 2 Figure 2: Representation of the geometry under consideration, with L = 50mm, = 25mm and e = 4.1mm.

Figure 3 :

 3 Figure 3: One realization of the Gaussian field {Ξ

  ), x ∈ B (r) } for an even ply, with L 1 = L 2 = 2 (left) and L 1 = L 2 = 10 (right).

Figure 4 :

 4 Figure 4: One realization of the Gaussian field {Ξ

  ), x ∈ B (r) } for an even ply, with (L 1 , L 2 ) = (2, 10) (left) and (L 1 , L 2 ) = (10, 2) (right).

Figure 5 :

 5 Figure 5: One realization of {G (r) 1 (x), x ∈ B (r) } for (L 1 , L 2 , δ) = (2, 2, 0.1) (left) and (L 1 , L 2 , δ) = (2, 2, 0.2) (right).

Figure 6 :

 6 Figure 6: One realization of {G

  ), x ∈ B (r) } for (L 1 , L 2 , δ) = (2, 10, 0.1) (left) and (L 1 , L 2 , δ) = (2, 10, 0.2) (right).

Figure 7 :

 7 Figure 7: Representation of the boundary conditions applied to the composite specimens.

θ θ 1

 1 = 15 • θ 2 = 30 • θ 3 = 60 • θ 4 = 75 • Realization ω exp ω s∈Texp u d (s) [mm] 0.123 0.116 0.158 0.150 0.485 0.745 0.715 0.614

Figure 8 :Figure 9 :

 89 Figure 8: Left: experimental vertical displacement prescribed on the top of the specimens. Right: experimental reaction force the top the specimens with respect to the prescribed displacement.

Figure 10 :

 10 Figure 10: One experimental realization of the strain field x → [E exp (s * ; x, θ)] 22 , for θ = 15 • , on the front (left) and back (right) surfaces.

Figure 11 :

 11 Figure 11: One experimental realization of the strain field x → [E exp (s * ; x, θ)] 22 , for θ = 30 • , on the front (left) and back (right) surfaces.

Figure 12 :

 12 Figure 12: One experimental realization of the strain field x → [E exp (s * ; x, θ)] 22 , for θ = 60 • , on the front (left) and back (right) surfaces.

Figure 13 :

 13 Figure 13: One experimental realization of the strain field x → [E exp (s * ; x, θ)] 22 , for θ = 75 • , on the front (left) and back (right) surfaces.

g 1 =Figure 14 :

 114 Figure 14: Left: final values of the cost function for 100 starting points obtained by Latin Hypercube Sampling. Filled bullets correspond to the cases for which the cost function is below 10 -10 . Right: convergence of the optimization algorithm for the 82th starting guess.

Figure 17 :

 17 Figure17: Identified optimal values of the exponents β 3 and β 4 , obtained for the 100 starting points. Filled bullets correspond to cases for which the final value of the cost function is below 10 -10 . In the case of β 5 , the limits of the y-axis have been reduced for readability.

Figure 18 :

 18 Figure 18: Graph of the log-likelihood function ζ → L(ζ). Note that marker size is proportional to the value of the cost function.

Figure 19 :

 19 Figure 19: One realization of the longitudinal and transversal Poisson coefficients (from left to right) in the first ply (with θ = 60 • , L 1 = L 2 = 11.18, δ = 0.1).

Figure 20 :

 20 Figure 20: One realization of the longitudinal and transversal Young moduli (from left to right) in the first ply (with θ = 60 • , L 1 = L 2 = 11.18, δ = 0.1).

Figure 21 :Figure 22 :

 2122 Figure 21: One realization of the longitudinal shear modulus in the first ply (with θ = 60 • , L 1 = L 2 = 11.18, δ = 0.1).

Table 2 :

 2 Notation related to kinematics.

Table 4 :

 4 Notation related to stochastic constitutive equations.

Table 5 :

 5 Summary of experimental parameters.

+ Cof([F ])e

2 ) .[START_REF] Iliopoulos | Meshless methods for full field displacement and strain measurement[END_REF] For any unit vector a, the terms [F ]a and Cof([F ])a represent the stretch along the direction defined by a and the deformation of an area with unit normal vector a, respectively. This interpretation serves as a physical basis to justify the use of the retained anisotropic contribution.

Acknowledgments

Authors BS and JG gratefully acknowledge partial support from the National Science Foundation under grant CMMI-1726403. Authors JM and AI acknowledge support for this work by the Office of Naval Research through the Naval Research Laboratory's core funding.

Appendix A. Definition of the anisotropic strain energy function

A generic form for anisotropic strain energy functions is proposed in [START_REF] Schröder | Anisotropic polyconvex energies on the basis of crystallographic motivated structural tensors[END_REF] as (see Eq. (3.7) in [START_REF] Schröder | Anisotropic polyconvex energies on the basis of crystallographic motivated structural tensors[END_REF]): ψ(I 1 , I 2 , I 3 , J 4 , J 5 ) = ψ iso (I 1 , I 2 , I 3 ) + ψ ani (I 1 , I 2 , I 3 , J 4 , J 5 ) .

(A. [START_REF] Mast | Characterization of strain-induced damage in composites based on the dissipated energy density part iii. general material constitutive relation[END_REF] where {I i } 3 i=1 and {J i } 5 i=4 are appropriate invariants, defined below. The isotropic contribution is specifically given by (see Eq. (3.12) in [START_REF] Schröder | Anisotropic polyconvex energies on the basis of crystallographic motivated structural tensors[END_REF]):

where

F and I 3 = (det([F ])) 2 ; {α i } 2 i=1 and {δ i } 2 i=1 are model parameters. Various forms for the anisotropic contribution were proposed in [START_REF] Schröder | Anisotropic polyconvex energies on the basis of crystallographic motivated structural tensors[END_REF]. In this paper, the form given in Tab. 1 is selected, with m = n = 1, ξ 11 = 1 and γ 11 = -1/2:

where

Using the following correspondence between the notations in [START_REF] Schröder | Anisotropic polyconvex energies on the basis of crystallographic motivated structural tensors[END_REF] and those considered in this paper:

where spatial dependency has been introduced. For homogeneous material properties, rigorous proofs of polyconvexity and coercivity can be found in [START_REF] Schröder | Anisotropic polyconvex energies on the basis of crystallographic motivated structural tensors[END_REF]. In order to accommodate the case of heterogeneous material parameters, the latter are assumed to be uniformly bounded from below (see Eq. ( 11)) in this paper, and an additional term ψ * (x) was added to bound the volumetric contribution from below. This leads to the final decomposition

of the final energy function, given by Eq. ( 9) (for a given ply).

Appendix B. Mathematical properties of the stored energy function

Under the constraints stated in Section 3.2.1, it can be shown that the stored energy function w (r) satisfies the following properties.

i. Strong polyconvexity. There exist a constant > 0 and a convex function f (r) (see below) such that

for all x ∈ B (r) and [F ] ∈ M 3 + . This generalized convexity condition is stronger than the usual polyconvexity condition [START_REF] Ball | Convexity conditions and existence theorems in nonlinear elasticity[END_REF][START_REF] Dacorogna | Direct methods in the calculus of variations[END_REF], due to the presence of the constant . In the present case, the latter can be chosen as = τ (r) 1 and the function f (r) can be defined as

2)

The convexity of the function

[H], δ) then follows from standard arguments [START_REF] Schröder | Invariant formulation of hyperelastic transverse isotropy based on polyconvex free energy functions[END_REF][START_REF] Schröder | Anisotropic polyconvex energies on the basis of crystallographic motivated structural tensors[END_REF].

ii. Coercivity. For all x in B (r) and [F ] in M 3 + , there exist k (r) 1 > 0, p 2 and q 3/2 such that

The strong polyconvexity implies that the stored energy function is uniformly elliptic, that is, it satisfies the (strong Legendre-Hadamard) condition

1 . Following existence theorems in finite elasticity [START_REF] Ball | Convexity conditions and existence theorems in nonlinear elasticity[END_REF][START_REF] Müller | On a new class of elastic deformations not allowing for cavitation[END_REF], the above properties ensure the existence of solutions to the nonlinear boundary value problem.