Jean Feydy
email: jean.feydy@ens.fr

Alain Trouvé
email: trouve@cmla.ens-cachan.fr

Global divergences between measures: from Hausdorff distance to Optimal Transport

Keywords: shape registration, kernel, energy distance, hausdorff distance, optimal transport, GPU

The data fidelity term is a key component of shape registration pipelines: computed at every step, its gradient is the vector field that drives a deformed model towards its target. Unfortunately, most classical formulas are at most semi-local: their gradients saturate and stop being informative above some given distance, with appalling consequences on the robustness of shape analysis pipelines. In this paper, we build on recent theoretical advances on Sinkhorn entropies and divergences [6] to present a unified view of three fidelities between measures that alleviate this problem: the Energy Distance from statistics; the (weighted) Hausdorff distance from computer graphics; the Wasserstein distance from Optimal Transport theory. The ε-Hausdorff and ε-Sinkhorn divergences are positive fidelities that interpolate between these three quantities, and we implement them through efficient, freely available GPU routines. They should allow the shape analyst to handle large deformations without hassle.

Introduction

Shape registration as a variational problem. Given a source shape A and a target B, a key problem in medical image analysis is to register the former onto the latter. That is, to estimate a mapping ϕ (a change of coordinates) that maps the source A into a model ϕ(A) which is "close enough" to the target.

Most classical registration algorithms strive to minimize an energy E(ϕ) = Reg(ϕ) regularizer + d (ϕ(A), B) fidelity which is the sum of a regularization term -encoding a prior on acceptable mappings -and a data attachment term -or fidelity -that measures how far the model ϕ(A) is from the target B.

The need for robust fidelities and gradients. Unfortunately, as of today, most fidelities can at best be described as semi-local. Relying on small convolution filters or kernel functions that saturate at long range [START_REF] Kaltenmark | A general framework for curve and surface comparison and registration with oriented varifolds[END_REF], they stop being informative when parts of the shapes are far away from each other. In recent years, finely crafted formulas have been proposed to alleviate this problem [START_REF] Lombaert | Spectral log-demons: diffeomorphic image registration with very large deformations[END_REF]; but they were probably too hard to implement and did not meet widespread adoption. As a result, most users today still rely on finely tuned coarse-to-fine schemes to register shape populations.

Contribution. At MICCAI 2017, we introduced the theory of Optimal Transport to the medical imaging community [START_REF] Feydy | Optimal transport for diffeomorphic registration[END_REF]. Leveraging the ideas and algorithms presented in [11], we showed that using globally optimal spring systems to drive a registration routine is tractable, and improves the robustness of pipelines to large deformations. The present paper is about taking advantage of new advances in the field [START_REF] Feydy | Sinkhorn entropies and divergences[END_REF] that let us bridge the gap between Optimal Transport and the standard shape analysis toolkit.

In section 1, we review the standard theory of measures and kernel distances (also known as blurred Sums of Squared Distances). We stress the relevance of the scale-invariant kernel k(x) = -x , which induces the global Energy Distance between shapes. We also notice that kernel distances rely on linear potentials (influence fields) generated by the shapes.

In section 2, we show how to use and compute non-linear potentials. We introduce a family of cheap fidelities between measures, the ε-SoftMin costs, that interpolate between the Energy Distance (ε = +∞) and the weighted Hausdorff distance (ε = 0) borrowed from computer graphics.

Finally, in section 3, we come back to the optimal transport cost and show that it is nothing but a "Hausdorff" distance under a mass repartition constraint. We interpret the celebrated Sinkhorn algorithm as a balancing scheme on distance fields and put an emphasis on two fidelities: the cheap ε-Hausdorff and the high-quality ε-Sinkhorn divergence, with a guarantee of positivity for both. Fig. 1: We focus this paper on the registration of thin segmented volumes (a, from the OsteoArthritis Initiative) encoded as measures on the ambient space (b). We provide efficient GPU routines to compute long-range gradients, from cheap kernel distances (c) to high-quality optimal transport plans (d).

In practice. Most importantly, we provide efficient CUDA routines -with Matlab, numpy and pytorch bindings -that can be used to implement these new data attachment terms. As shown in section 4, our KeOps library [START_REF] Feydy | Kernel operations on the gpu, with autodiff, without memory overflows. www.kernel-operations[END_REF] allows users to process curves, surfaces and segmentation maps with up to 100,000 actives vertices on a cheap laptop's GPU. Our code is freely available:

Please visit github.com/jeanfeydy/global-divergences.

These new tools fit seamlessly into the standard shape analyst's toolkit; they should help the reader to improve with little to no overhead the robustness to large deformations of its shape analysis pipeline.

Representing shapes as measures on a space of features

In this paper. We choose to focus this paper on a setting that is understood well by all researchers in medical image analysis: the registration of normalized density maps. Our source bitmap A (in red) and target B (in blue) will be encoded as measures

α = N i=1 α i δ xi and β = M j=1 β j δ yj , with N i=1 α i = 1 = M j=1 β j ,
where the x i 's (respectively y j 's) are the coordinates of the N (resp. M) nonzero pixels of A (resp. B), with positive weights α i (resp. β j ) summing up to one.

In most figures, we will display the gradient ∇ xi d(α, β) of a fidelity "d" as a green vector field supported by the x i 's. This descent direction is meant to be used by registration algorithms and is thus the primary information to look at in our pictures. In the background, depending on the section, we also display the level lines of the linear potential "k (α -β)" (in blue) or of the influence fields "a" (in red) and "b" (in blue) -more about that later.

Extensions. The results presented in this paper can be extended to other use cases fairly easily. First, we may wish to use an image-based registration of segmentation maps instead of the mass preserving "Jacobian-free" action. To do so, we should simply compute the gradient ∇ αi d(α, β) of fidelities with respect to the weights of the atomic dirac masses; the presence of long-range interactions is equally important to the robustness of the registration algorithm, with mass contraction (i.e. deletion) replacing the spreading out phenomenon observed in Figure 3.(a-b).

Most of our results still hold when the source and the target don't have the same mass -the only noticeable changes would be located in section 3, and we recommend [11] as an introduction to the theory of unbalanced optimal transport. Going further, these new tools and GPU routines can also be used to handle fiber tracks, curves and surfaces through the varifold framework presented in [START_REF] Kaltenmark | A general framework for curve and surface comparison and registration with oriented varifolds[END_REF].

Notations. In order to let our results be useful to researchers working with curves and surfaces -which are best represented as measures on a product space (position,orientation,curvature) -we will refer to the ambient space R 2 or R 3 as to an abstract feature space X . The letters x, y and z will denote points in the feature space, while α, β or µ stand for finitely supported positive measures; finally, a, b and m denote real-valued functions on X understood as influence fields generated by their respective measures.

If (z i ) i∈[[1,N]
] is a collection of N points in X and if m : X → R is a function on the feature space, we will also write "m zi " to denote the length-N vector (m

(z i )) i∈[[1,N]] of values of m sampled on the point cloud z i .
Finally, if µ = N i=1 µ i δ zi is a finitely supported measure and if m : X → R is a function on the feature space, we will write

µ , m = ( µ i | m zi ) = N i=1 µ i m(z i ).
Here,

(µ i ) i∈[[1,N]] and (m zi ) i∈[[1,N]] are two vectors of R N : the measure-function duality bracket • , • is thus understood as a simple scalar product ( • | • ) in R N .

Kernel distances

If α and β represent two shapes in the feature space X , using standard informationtheoretic fidelities such as the symmetrised Kullback-Leibler divergence

KL sym (α, β) = 1 2 KL(α, β) + 1 2 KL(β, α) = 1 2 α -β , log dα dβ 0 (1)
is not recommended: shape analysis routines should take into account the geometry of the feature space.

Kernel norms. A common way of doing so is to endow the feature space X with a symmetric kernel function k : X × X → R and to use

d k (α, β) = 1 2 α -β , k (α -β) = 1 2 α -β , b k -a k , (2)
where

a k (z) = -(k α)(z) = - N i=1 α i k(x i , z) (3)
and b k (z) = -(k β)(z) = - M j=1 β j k(y j , z). (4)
In practice, these summations can be implemented as matrix-vector products or, as advocated in Figure 11, by using the online map-reduce routines of the KeOps library [START_REF] Feydy | Kernel operations on the gpu, with autodiff, without memory overflows. www.kernel-operations[END_REF]. Popular choices include the Gaussian and Laplacian kernels:

Gaussian σ (x -y) = exp(-x -y 2 /σ 2 )
and Laplacian σ (x -y) = exp(-x -y /σ).

However, as ∇ xi d k (α, β) is given by the gradient of the linear potential (b k -a k ) sampled on the x i 's and weighted by the α i 's, we argue in Figures 23that a more robust baseline could be given by the Energy Distance kernel from statistics [START_REF] Székely | Energy statistics: A class of statistics based on distances[END_REF]:

Energy(x -y) = -x -y .

Computing non-linear potentials

The log-sum-exp trick. In order to build tractable algorithms, restricting ourselves to potentials a and b that depend linearly on the measures α and β seems to be a necessary evil... But we can go further. Indeed, on top of the "summation" operation of Eqs. [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF][START_REF] Fedorov | 3d slicer as an image computing platform for the quantitative imaging network[END_REF], we can implement on the GPU another differentiable reduction operator: the log-sum-exp or SoftMax, defined through log

N i=1 exp(v i ) = V + log N i=1 exp(v i -V
), with V = max i v i taken out of the expression for numerical stability. The KeOps library implements an online variant of this "max-factorization" trick, and lets us scale this operation to large values of N -see Figure 11.

Definition. Then, we propose to endow the ambient space X with a symmetric cost function C : X × X → C(x, y) -say, x -y -a regularization strength ε > 0 and a kernel function

k ε = exp(-1 ε C( • , • )) to define min ε x∼α C(x, z) = -ε log(k ε α)(z) = -ε log N i=1 exp log(α i ) -1 ε C(x i , z
) . Mimicking Eq. ( 2), we propose to see the SoftMin functions as non-linear influence fields, analogous to the linear potentials a k and b k . Hence, we introduce the ε-SoftMin cost through

d ε-SoftMin (α, β) = 1 2 α -β , b ε -a ε = 1 2 α i | b ε xi -a ε xi + 1 2 β j | a ε yj -b ε yj , where a ε (z) = min ε x∼α C(x, z) and b ε (z) = min ε y∼β C(y, z).
Interpretation. Simple calculations show that if C(x, y) = x -y , the ε-SoftMin cost converges towards the Energy distance as ε goes to infinity. At the other end of the spectrum, if C(x, y) 0 with equality if x = y,

d ε-SoftMin (α, β) ε→0 ----→ 1 2 N i=1 α i min j C(x i , y j ) + 1 2 M j=1 β j min i C(x i , y j )
As shown in Figures 56, the SoftMin operators is thus allowing us to interpolate between statistics and computer graphics.

Positivity. Unfortunately, one cannot guarantee the positivity of the ε-SoftMin fidelity: linearizing the cost, we find pairs of measures such that d ε-SoftMin (α + δα, α) < 0. However, if λ is a reference measure on the feature space X (say, the Lebesgue measure on R D ), then

ε KL sym ((k ε α) • λ, (k ε β) • λ) = 1 2 λ • k ε (α -β) , ε log k ε α k ε β = 1 2 λ • k ε (α -β) , b ε -a ε 0. In practice, if (α -β) is close enough to its ε-blurred image λ • k ε (α -β), d ε-SoftMin (α, β) is thus positive too.

Balancing distance fields: the Sinkhorn algorithm

As seen in Figure 7, adding a mass distribution constraint to SoftMin distances can improve the quality of our descent directions. Thankfully, this is now possible thanks to the theory of Optimal Transport, which generalizes the Wasserstein distance -we recommend the recent handbook [11] for reference.

Primal "Monge" problem. As illustrated in Figure 8, Optimal Transport is about solving a convex registration problem: for ε > 0, we strive to minimize a primal cost OT ε defined through

OT ε (α, β) = min πx i ,y j ∈R N×M 0 i,j π xi,yj C(x i , y j ) transport cost π, C + ε i,j π xi,yj log π xi,yj α i β j -π xi,yj + α i β j entropic regularization, ε KL( π, α ⊗ β)
, under a linear constraint -α should be fully transported onto β:

∀ i ∈ [[1, N]] , α i = M j=1 π xi,yj and ∀ j ∈ [[1, M]] , β j = N i=1 π xi,yj .
The Sinkhorn algorithm. Strong duality holds on OT ε . The major contribution from [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF] was to show that the dual problem can be solved efficiently on the GPU. In a nutshell, we run the following algorithm:

Algorithm 1 Sinkhorn Iterative Algorithm: Sink(α i , x i , β j , y j )

Parameters : symmetric cost C : (x, y) → C(x, y), regularization ε > 0 Input : source α = N i=1 αiδx i , target β = M j=1 βjδy j Output : influence fields a α→β and b β→α , sampled on the yj's and xi's respectively 1: ay j ← zeros(M) ; bx i ← zeros(N)

Vectors of size M and N, respectively 2: for it = 1 to nits do

In practice, nits = 10 to 30 is enough 3:

ay j ← minε,x∼α[C(x, yj) -b(x)] = -ε log N i=1 exp log(αi) -1 ε (C(xi, yj) -bx i ) 4: bx i ← min ε,y∼β [C(xi, y) -a(y)] = -ε log M j=1 exp log(βj) -1 ε (C(xi,
) = Sink(α i , x i , β j , y j ), we get OT ε (α, β) = α , b β→α + β , a α→β = α i | b β→α xi + β j | a α→β yj .
Interpretation. The Sinkhorn algorithm is a block-coordinate ascent on the dual variables. Mathematically speaking, these are Lipschitz functions defined on the ambient space X and sampled on the measures' supports. As illustrated in Fig. 9 and detailed in [START_REF] Feydy | Sinkhorn entropies and divergences[END_REF], we propose to understand them as influence fields a α→β and b β→α that encode an implicit transport plan

π = exp 1 ε (b ⊕ a -C) • α ⊗ β.
Fig. 8: Looking for a low-cost mapping, from α to β. As we solve the Optimal Transport problem, we find a transport plan π whose marginals are equal to α and β respectively. In practice, we solve a regularized problem whose dual solution is easy to compute on the GPU. Fig. 9: The (standard) Sinkhorn algorithm brings balance to the force. On top of α, β, a α→β (in red) and b β→α (in blue), we display the mean "springs" linking the x i 's to β (in red) and the y j 's to α (in blue). Algorithm 1 is all about normalizing the blue (line 3) and red (line 4) springs until reaching equilibrium. Towards a positive Optimal Transport cost. In [START_REF] Feydy | Optimal transport for diffeomorphic registration[END_REF], we advocated the use of OT ε (α, β) as a data attachment term for diffeomorphic registration. Unfortunately though, just as if we used an ε-SoftMin fidelity, the minimum of the loss functional α → OT ε (α, β) is not reached when α is exactly equal to β.

With collaborators [START_REF] Feydy | Sinkhorn entropies and divergences[END_REF], we thus decided to shift our attention towards a new geometric entropy:

F ε (α) = -1 2 OT ε (α, α) = ε min µi∈R N 0 α , log dα dµ + 1 2 µ , k ε µ -1 2 , (5)
with µ = N i=1 µ i δ xi -this identity stands thanks to a change of variable "µ i = exp(a xi /ε) α i ". We let SymSink denotes the symmetrized Sinkhorn algorithm:

ax i ← zeros(N)
Vector of size N 2: for it = 1 to nits -1 do

In practice, nits = 3 is enough 3:

ax i ← 1 2 (ax i + minε,x∼α [C(xi, x) -a(x)]) = 1 2 ax i -ε log N k=1 exp log(α k ) -1 ε (C(xi, x k ) -ax k ) 4: a x i ← minε,x∼α[C(xi, x) -a(x)] = -ε log N k=1 exp log(α k ) -1 ε (C(xi, x k ) -ax k ) 5: a y j ← minε,x∼α[C(yj, x) -a(x)] = -ε log N k=1 exp log(α k ) -1 ε (C(yj, x k ) -ax k )
6: return a y j , a x i Vectors of size M and N, respectively

The ε-Sinkhorn divergence. Then, we can define

d ε-Hausdorff (α, β) = 1 2 α -β , ∇F ε (α) -∇F ε (β) = 1 2 α i | b β↔β xi -a α↔α xi + 1 2 β j | a α↔α yj -b β↔β yj , d ε-Sinkhorn (α, β) = OT ε (α, β) -1 2 OT ε (α, α) -1 2 OT ε (β, β) = α i | b β→α xi -a α↔α xi + β j | a α→β
) = Sink(α i , x i , β j , y j ).
The ε-Hausdorff divergence is the symmetrized Bregman divergence associated to F ε on the space of probability measures on X , and can be shown to behave like the ε-SoftMin cost at the small and large ε limits. Meanwhile, the ε-Sinkhorn divergence is an "unbiased" Optimal Transport cost that has been recently introduced in the Machine Learning community [START_REF] Genevay | Learning generative models with sinkhorn divergences[END_REF].

The intuition here is that since OT ε (α, β) converges towards a kernel scalar product α, C β when ε goes to infinity, adding the self-correlation corrective terms lets us converge towards a genuine kernel squared norm 1 2 α -β 2 -Csay, the Energy Distance if C(x, y) = x -y . Most importantly, we are able to prove that both formulas define positive divergences for ε > 0 : Proof. The proof of this result is given in [START_REF] Feydy | Sinkhorn entropies and divergences[END_REF].

In a nutshell: the first inequality relies on the positivity of the kernel k ε , as it ensures the convexity of the potential F ε -Eq. ( 5) -and the positivity of the associated Bregman divergence. The second inequality derives from the convexity of OT ε (α, β) with respect to α and β varying independently.

Conclusion

Overview. All things considered, we introduced three positive divergences to the shape analysis community: the cheap and global Energy Distance; the highquality ε-Sinkhorn cost; and, sitting in-between, a brand new ε-Hausdorff divergence inspired by computer graphics. All of them define well-posed, differentiable loss functions for registration problems.

As we linked these theories with each other in sections 2 and 3, we were able to provide important theoretical guarantees and efficient GPU routines. In practice, we advocate the use of the PyTorch + KeOps combination [START_REF] Paszke | Automatic differentiation in pytorch[END_REF][START_REF] Feydy | Kernel operations on the gpu, with autodiff, without memory overflows. www.kernel-operations[END_REF] that provides automatic differentiation and scalability to shapes with up to 100,000 active vertices. Going further. Now, which one of these formulas should we use in practice? As seen in Figure 10, using an ε-Sinkhorn divergence is equivalent to performing a full convex registration -with no guarantee of topology preservation -every time we need a descent direction... Do we really need to go that far?

The answer to this question is highly dependent on the remainder of the registration pipeline. In months to come, we thus plan to test our new fidelities in a wide range of settings -from standard LDDMM to Deep Learning based methods -as we strive to provide our colleagues with reliable tools.

Fig. 2 :

 2 Fig. 2: The linear potential k (α -β), for standard kernel functions. Here, α and β are sampled from the standard Lebesgue measures on the segments [.2, .35] and [.65, .8], respectively. Out of these three curves, the third is the only one whose (minus) gradient always points from α towards β.

Fig. 3 :

 3 Fig. 3: The Energy Distance is scale-invariant and robust to large deformations. This is the 2D equivalent of Figure 2, with level lines of k (α -β) displayed in the background. Notice the spreading out effect in (a-b).

Fig. 4 :

 4 Fig. 4: Linear potentials can only take you so far. (a) As it faces a mass imbalance, the global gradient of the Energy Distance tries to split up the largest red mass into pieces. (b-c) The SoftMin fidelity, introduced in section 2, allows us to induce a more "focused" behavior into our algorithms.

Fig. 5 :

 5 Fig. 5: The SoftMin operator. The log-sum-exp trick allows us to interpolate between two kinds of distance fields to a measure: the "Energy potential" | • | α -for ε = +∞ -and the distance field to the support {.25, .75} of α -for ε = 0.

Fig. 6 :

 6 Fig. 6: The ε-SoftMin fidelity interpolates between the Energy Distance and a weighted Hausdorff distance between the supports [1]. Here, we use the simple Euclidean cost C(x, y) = x -y .

1 Fig. 7 :

 17 Fig. 7: Naive projection isn't the panacea. Using an ε-SoftMin cost is equivalent to encoding our shapes through their distance images -the so-called influence fields, displayed in the background. (a) Unfortunately, such a cost is prone to giving a disproportionate importance to the extremities of both shapes, as points are only influenced by their nearest neighbors. (b) Presented in section 3, the Sinkhorn loop lets us introduce a mass distribution constraint to alleviate this problem: we shift the influence fields a (in red) and b (in blue) to retrieve balanced gradient fields. Figures computed with C(x, y) = x -y 2 .

  yj) -ay j ) 5: return ay j , bx i Vectors of size M and N, respectively And at convergence, with (a α→β yj , b β→α xi

1 Fig. 10 :

 110 Fig.10: Computing ε-Hausdorff and ε-Sinkhorn divergences. On this page, we use a quadratic cost C(x, y) = x -y 2 so that ε is homogeneous to a squared distance. As evidenced in (b), the gradient of the ε-Sinkhorn divergence is good enough for one-shot registration in simple cases.

Algorithm 2

 2 Symmetric Sinkhorn Algorithm: SymSink(α i , x i , y j )Parameters : symmetric cost C : (x, y) → C(x, y), regularization ε > 0 Input : source α = N i=1 αiδx i , target point cloud (yj) j∈[[1,M]]Output : influence field a α↔α sampled on the xi's and the yj's 1:

yj

  (α i , x i , y j ), (b β↔β xi , b β↔β yj ) = SymSink(β j , y j , x i ) and (a α→β yj , b β→α xi

1

 1 Fig. 11: The KeOps library lets shape analysis routines scale up to real data. Performances on a cheap laptop's GPU (GTX 960M). (a) As it provides CUDA routines for online map-reduce operations, our "KErnelOPerationS" library -developed with Benjamin Charlier and Joan A. Glaunès -allows Matlab, numpy and pytorch users to compute huge N-by-N convolutions without having to store large kernel matrices in the GPU memory. (b) Experiments performed on point clouds in R 3 , endowed with a Euclidean cost C(x, y) = x -y .

Fig. 12 :

 12 Fig.12: On real data. Our routines could be used to registrate thin structures such as these knee caps from the OsteoArthritis Initiative -special thanks to Zhenlin Xu and Marc Niethammer for letting us know about this dataset. Here, the source and target volumes are respectively made up of 52,319 and 34,966 voxels -out of a 192-192-160 volume. As advertised in Figure11, this Energy Distance's gradient was computed in half a second on the author's laptop.


Martin Hirsekorn
email: martin.hirsekorn@onera.fr

Lionel Marcin

Thierry Godon

Multi-scale modeling of the viscoelastic behavior of 3D woven composites

Keywords: A. 3-dimensional reinforcement, B. Thermomechanical, B. Stress relaxation, C. Computational modeling

A method is presented to predict numerically the homogenized viscoelastic behavior of 3D woven composites using only its constituents' behavior. It is based on elastic homogenizations applied to the Laplace-Carson transform of the time-dependent viscoelastic behavior of the constituents. Two scale chages are necessary: from micro-to meso-scale to obtain the homogenized behavior of the consolidated yarns, and from mesoto macro-scale. The temperature and cure dependent viscoelastic behavior of the matrix is identified from experimental data, using the time-temperature superposition principle with the cure dependent glass transition temperature as reference temperature. The meso-scale representative unit cell of the composite is extracted from X-ray microtomography images. The homogenized viscoelastic behavior is used to calculate the evolution of the apparent moduli of the composite with temperature between -50 • C and 200 • C. The results are in good agreement with experimental data over the temperature range where the matrix behavior was properly identified.

Introduction

During the manufacturing process of composites with thermosetting polymer matrices, residual stresses in the material cannot be avoided, because the thermosetting polymers are usually cured at much higher temperatures than the typical in-service temperature of the composite parts. These residual stresses may cause damage and influence its evolution [START_REF] Fiedler | The influence of thermal residual stresses on the transverse strength of CFRP using FEM[END_REF][START_REF] Maier | Performance enhancements of polymer-matrix composites by changing of residual stresses[END_REF][START_REF] Correa | Effect of thermal residual stresses on the matrix failure under transverse compression at micromechanical level -A numerical and experimental study[END_REF][START_REF] Agius | Rapidly cured epoxy/anhydride composites: Effect of residual stress on laminate shear strength[END_REF] and change the final shape of composite parts [START_REF] Fernlund | Experimental and numerical study of the effect of cure cycle, tool surface, geometry, and lay-up on the dimensional fidelity of autoclave-processed composite parts[END_REF][START_REF] Wisnom | Mechanisms generating residual stresses and distortion during manufacture of polymer-matrix composite structures[END_REF][START_REF] Ersoy | Development of spring-in angle during cure of a thermosetting composite[END_REF][START_REF] Causse | Spring-in behavior of curved composites manufactured by Flexible Injection[END_REF][START_REF] Kappel | Manufacturing distortions of a CFRP box-structure -A semi-numerical prediction approach[END_REF].

The major part of the residual stresses is due to a mismatch in the thermal expansion coefficients of the constituent materials (fibers and matrix), and thus generated during the cooling of the composite after the end of the curing cycle [START_REF] Fernlund | Experimental and numerical study of the effect of cure cycle, tool surface, geometry, and lay-up on the dimensional fidelity of autoclave-processed composite parts[END_REF][START_REF] Darrow | Isolating components of processing induced warpage in laminated composites[END_REF][START_REF] Wisnom | Mechanisms generating residual stresses and distortion during manufacture of polymer-matrix composite structures[END_REF]. Furthermore, a thermosetting polymer shrinks during curing [START_REF] Li | In-situ measurement of chemical shrinkage of MY750 epoxy resin by a novel gravimetric method[END_REF][START_REF] Billotte | Chemical shrinkage and thermomechanical characterization of an epoxy resin during cure by a novel in situ measurement method[END_REF], which may also contribute to the formation of residual stresses [START_REF] Russell | A new method to reduce cure-induced stresses in thermoset polymer composites, Part III: Correlating stress history to viscosity, degree of cure, and cure shrinkage[END_REF][START_REF] Darrow | Isolating components of processing induced warpage in laminated composites[END_REF][START_REF] Wisnom | Mechanisms generating residual stresses and distortion during manufacture of polymer-matrix composite structures[END_REF].

While the mechanical behavior of carbon or glass fibers is linear elastic in a good approximation and almost independent of temperature, at least below the glass transition temperature (T g ) of typical thermosetting polymers (usually not exceeding about 200 • C), the mechanical behavior of the polymer matrix depends strongly on temperature and on the degree of cure [START_REF] Kim | Stress relaxation behavior of 3501-6 epoxy resin during cure[END_REF][START_REF] Simon | Modeling the evolution of the dynamic mechanical properties of a commercial epoxy during cure after gelation[END_REF][START_REF] O'brien | Viscoelastic properties of an epoxy resin during cure[END_REF][START_REF] Ruiz | Numerical analysis of cure temperature and internal stresses in thin and thick RTM parts[END_REF][START_REF] Ersoy | Development of the properties of a carbon fibre reinforced thermosetting composite through cure[END_REF][START_REF] Billotte | Chemical shrinkage and thermomechanical characterization of an epoxy resin during cure by a novel in situ measurement method[END_REF]. Moreover, it is time dependent, i.e., creep and relaxation phenomena are observed, which may have a significant impact on the residual stresses and the final shape of composite parts [START_REF] Russell | A new method to reduce cure-induced stresses in thermoset polymer composites, Part III: Correlating stress history to viscosity, degree of cure, and cure shrinkage[END_REF][START_REF] Hobbiebrunken | Microscopic yielding of CF/epoxy composites and the effect on the formation of thermal residual stresses[END_REF][START_REF] Ding | A thermo-viscoelastic model of process-induced residual stresses in composite structures with considering thermal dependence[END_REF][START_REF] Benavente | Viscoelastic distortion in asymmetric plates during post curing[END_REF]. As a consequence, residual stresses can be reduced by adapting the cure cycle such that relaxation phenomena are intensified [START_REF] Russell | A new method to reduce cure-induced stresses in thermoset polymer composites, Part III: Correlating stress history to viscosity, degree of cure, and cure shrinkage[END_REF][START_REF] Kim | Reduction of residual stresses in thick-walled composite cylinders by smart cure cycle with cooling and reheating[END_REF].

Since in the aeronautical industry highly precise part shapes are required, it is currently common practice to adapt the mold shape several times until the composite part fits the shape tolerances. This procedure of designing a composite part is very expensive, but can be greatly optimized if precise modeling tools are available that predict the generation of residual stresses and the final shape of composites with complex reinforcements, such as three-dimensional (3D) interlock weavings. To improve the accuracy of currently available modeling tools, it is necessary to take into account creep and relaxation and their dependence on temperature and degree of cure during the modeling of the curing and cooling process. The experimental characterization of a composite material required to determine the model parameters is complex and expensive, as 3D strain rate, temperature, and cure dependent mechanical behaviors have to be measured. Furthermore, these properties have to be re-identified, each time the reinforcement architecture is changed. Significant cost reduction is possible if only the constituent materials are fully characterized and the composite properties are predicted using multi-scale models.

At small strains, creep, relaxation, and strain rate dependence of polymer resins are well described by viscoelastic models of the generalized Maxwell type [START_REF] Biot | Theory of stress-strain relations in anisotropic viscoelasticity and relaxation phenomena[END_REF][START_REF] Kim | Stress relaxation behavior of 3501-6 epoxy resin during cure[END_REF][START_REF] Simon | Modeling the evolution of the dynamic mechanical properties of a commercial epoxy during cure after gelation[END_REF][START_REF] O'brien | Viscoelastic properties of an epoxy resin during cure[END_REF]. This kind of model was used by several authors to simulate the generation of residual stresses in composites with both thermosetting and thermoplastic polymer matrices [START_REF] Sunderland | A thermoviscoelastic analysis of process-induced internal stresses in thermoplastic matrix composites[END_REF][START_REF] Svanberg | Prediction of shape distortions Part I. FE-implementation of a path dependent constitutive model[END_REF][START_REF] Clifford | Thermoviscoelastic anisotropic analysis of process induced residual stresses and dimensional stability in real polymer matrix composite components[END_REF][START_REF] Ding | A thermo-viscoelastic model of process-induced residual stresses in composite structures with considering thermal dependence[END_REF][START_REF] Zhang | Residual stresses created during curing of a polymer matrix composite using a viscoelastic model[END_REF]. In all these works, the time-temperature superposition principle used, as it was shown to describe well the temperature dependence of the viscoelastic behavior of polymers [START_REF] Kim | Stress relaxation behavior of 3501-6 epoxy resin during cure[END_REF][START_REF] Simon | Modeling the evolution of the dynamic mechanical properties of a commercial epoxy during cure after gelation[END_REF][START_REF] O'brien | Viscoelastic properties of an epoxy resin during cure[END_REF][START_REF] Zarrelli | Toward a constitutive model for cure-dependent modulus of a high temperature epoxy during the cure[END_REF]. The influence of the degree of cure was introduced in some works as an explicit dependence of the relaxation times on degree of cure [START_REF] Kim | Stress relaxation behavior of 3501-6 epoxy resin during cure[END_REF][START_REF] O'brien | Viscoelastic properties of an epoxy resin during cure[END_REF][START_REF] Zhang | Residual stresses created during curing of a polymer matrix composite using a viscoelastic model[END_REF], but it can also be modeled through a cure dependent T g used as reference temperature in the time-temperature superposition principle [START_REF] Simon | Modeling the evolution of the dynamic mechanical properties of a commercial epoxy during cure after gelation[END_REF][START_REF] Svanberg | Prediction of shape distortions Part I. FE-implementation of a path dependent constitutive model[END_REF][START_REF] Zarrelli | Toward a constitutive model for cure-dependent modulus of a high temperature epoxy during the cure[END_REF][START_REF] Courtois | Viscoelastic behavior of an epoxy resin during cure below the glass transition temperature: Characterization and modeling[END_REF].

In the present study we focus on composites with 3D woven reinforcements [START_REF] Mouritz | Review of applications for advanced three-dimensional fibre textile composites[END_REF][START_REF] Kamiya | Some recent advances in the fabrication and design of threedimensional textile preforms: a review[END_REF][START_REF] Ansar | Modeling strategies of 3D woven composites: A review[END_REF]. This type of composites became an interesting option for industrial applications since the development of modern looms capable of weaving complex reinforcement structures [START_REF] Mouritz | Review of applications for advanced three-dimensional fibre textile composites[END_REF][START_REF] Kamiya | Some recent advances in the fabrication and design of threedimensional textile preforms: a review[END_REF]. It is therefore possible to manufacture complex composite parts with very few assembling steps and thus short production times [START_REF] Mountasir | Competitive manufacturing of 3D thermoplastic composite panels based on multi-layered woven structures for lightweight engineering[END_REF], reducing the number of joints, which are often critical for damage onset. Moreover, composites with 3D woven reinforcements are less prone to delamination and thus more resistant to damage caused by impact [START_REF] Mouritz | Review of applications for advanced three-dimensional fibre textile composites[END_REF][START_REF] Gerlach | In-plane and through-thickness properties, failure modes, damage and delamination in 3D woven carbon fibre composites subjected to impact loading[END_REF][START_REF] Turner | Collapse of 3D orthogonal woven carbon fibre composites under in-plane tension/compression and out-of-plane bending[END_REF] compared to multi-layer composites of 2D weavings or unidirectional (UD) plies. In composites with woven reinforcements, usually three characteristic scales are distinguished [START_REF] Lomov | Meso-FE modelling of textile composites: Road map, data flow and algorithms[END_REF]: the micro-scale, at which each fiber and the matrix is modeled individually, the meso-scale, at which the reinforcement architecture is modeled and yarns consisting of several thousand fibers are treated as homogeneous material, and the macro-scale, at which the structure is modeled and the composite material is treated as homogeneous.

Modeling the generation of residual stresses in these materials is therefore a multi-scale problem, and suitable techniques are needed to determine the homogenized behavior of a heterogeneous material at the next higher scale. A common technique to calculate the homogenized behavior of heterogeneous viscoelastic materials is to apply the Laplace-Carson (LC) transform to the integral [START_REF] Brenner | A "quasi-elastic" affine formulation for the homogenised behaviour of nonlinear viscoelastic polycrystals and composites[END_REF][START_REF] Lévesque | Numerical inversion of the Laplace-Carson transform applied to homogenization of randomly reinforced linear viscoelastic media[END_REF][START_REF] Ricaud | Effective properties of linear viscoelastic heterogeneous media: Internal variables formulation and extension to ageing behaviours[END_REF][START_REF] Pichler | Upscaling of viscoelastic properties of highly-filled composites: Investigation of matrix-inclusiontype morphologies with power-law viscoelastic material response[END_REF][START_REF] Mourid | Comparison between analytical and numerical predictions for the linearly viscoelastic behavior of textile composites[END_REF] or the differential [START_REF] Beurthey | Structural morphology and relaxation spectra of viscoelastic heterogeneous materials[END_REF][START_REF] Thionnet | Prévision du comportement en fluage viscoélastique endommageable des structures à matrice organique[END_REF] formulation of the viscoelastic behavior to obtain a linear relationship between the transformed stress and strain history (more details will be given in section 3). In this way, classical elastic homogenization methods can be used to determine the homogenized viscoelastic behavior.

In the present article, a homogenization strategy for composites with 3D woven reinforcements is presented based on the LC transform. It allows obtaining the viscoelastic behavior of the composite material from the behavior of its constituents, including the dependence on temperature and degree of cure. In section 2, the formulation of the viscoelastic models used in this work is presented. A 3D generalization of the viscoelastic model developed by Courtois et al. [START_REF] Courtois | Viscoelastic behavior of an epoxy resin during cure below the glass transition temperature: Characterization and modeling[END_REF] for the polymer matrix of the studied composite is given in section 2.2. The proposed homogenization strategy is described in section 3 and validated numerically at the microscopic scale by comparison to full-field Finite Element (FE) calculations. In section 4, the strategy is applied to the case of a 3D woven composite using a representative unit cell (RUC) of the reinforcement architecture obtained from X-ray micro-computed tomography (µ-CT). The model is again validated numerically by comparison to full-field FE calculations. The obtained homogenized viscoelastic behavior of the composite is used to calculate the temperature dependence of the apparent elastic moduli of the composite. The results are compared to experimental observations.

Viscoelastic constitutive model

General anisotropic formulation

To model the viscoelastic behavior of the composite and its constituents at the different scales (yarns, matrix) the generalized Maxwell model is used. This model can be related to the orthotropic generalization [START_REF] Sawant | A thermo-mechanical viscoelastic analysis of orthotropic materials[END_REF] of Schapery's [START_REF] Schapery | On the characterization of nonlinear viscoelastic materials[END_REF] integral form of viscoelastic behaviors: For a material that is in a completely relaxed state at the time t 0 , the stress response to a strain history starting at t 0 can be written in the integral form

σ(t) = t t0 E(t -t ) : ε(t )dt ( 1
)
where E is the relaxation stiffness tensor (a 3D generalization of the relaxation modulus) and ε(t ) is the strain rate at a time between t 0 and t. In Eq. ( 1) and in the folllowing, single underlining indicates second order tensors, double underlining indicates fourth order tensors, and the colon indicates contraction over two indices.

The relaxation modulus is often written in terms of a Prony series [START_REF] Kim | Stress relaxation behavior of 3501-6 epoxy resin during cure[END_REF][START_REF] Simon | Modeling the evolution of the dynamic mechanical properties of a commercial epoxy during cure after gelation[END_REF][START_REF] O'brien | Viscoelastic properties of an epoxy resin during cure[END_REF], which can be done in the same way for E:

E(t) = E ∞ + N k=1 E k e - t-t 0 τ k (2)
where τ k is called the relaxation time of the k th relaxation mechanism. Using this expression, the integral form resolves to

σ(t) = E 0 : ε(t) - N k=1 E k : ε ve k (t) (3)
with the instantaneous elastic stiffness tensor

E 0 = E ∞ + N k=1 E k (4)
The tensorial internal variables ε ve i associated to the relaxation times τ i are called history variables [START_REF] Sawant | A thermo-mechanical viscoelastic analysis of orthotropic materials[END_REF]. Their use avoids the integration over the strain history (Eq. 1).

The evolution in time of the ε ve i is given by differential equations of the form

εve k (t) = 1 τ k (ε(t) -ε ve k (t)) (5)
They describe the evolution in time of the strains of the viscoelastic elements in a generalized Maxwell model. This means that under the assumption of a Prony series description for the relaxation stiffness, the integral form is equivalent to a generalized Maxwell model, with each term in the Prony series corresponding to one viscoelastic Maxwell element (with a 3D elastic behavior described by the tensor E k ). This form of an anisotropic generalized Maxwell model is used in the following to model the viscoelastic behavior of the composite and its constituents. The composite behavior is orthotropic, the yarns obey transverse isotropy around the fiber axis, and the matrix is isotropic. In the respective viscoelastic models, the tensors E k and E ∞ take the appropriate form with 9 (orthotropy), 5 (transverse isotropy), or 2 (isotropy) independent coefficients.

Specification for the polymer matrix

The polymer used as a matrix in the composite material described in section 4 is an isotropic material.

The fourth order tensors in Eq. ( 3) can thus be written as

E k = K k 1 3 1 ⊗ 1 + G k 1 - 1 3 1 ⊗ 1 E ∞ = K ∞ 1 3 1 ⊗ 1 + G ∞ 1 - 1 3 1 ⊗ 1 ( 6
)
where 1 is the second order identity tensor and 1 the fourth order identity tensor with minor symmetry.

Courtois et al. [START_REF] Courtois | Viscoelastic behavior of an epoxy resin during cure below the glass transition temperature: Characterization and modeling[END_REF] carried out relaxation tests on fully and partially cured specimens of the epoxy resin used in the studied composite at different temperatures below its glass transition temperature (T g ). Based upon their data, they proposed a 1D model for the tensile viscoelastic behavior of the resin as a function of temperature and degree of cure. In a composite material, the matrix is exposed to complex 3D strain states. To predict the homogenized behavior of the composite, we therefore need a full 3D model of the viscoelastic behavior of the matrix. If we use 1D model of Courtois et al. [START_REF] Courtois | Viscoelastic behavior of an epoxy resin during cure below the glass transition temperature: Characterization and modeling[END_REF] directly to model the evolution of the tensile modulus in a 3D model, this would imply the assumption of a constant Poisson ratio during relaxation. However, it was observed by Tschoegl et al. [START_REF] Tschoegl | Poisson's ratio in linear viscoelasticity -A critical review[END_REF], O'Brien et al. [START_REF] O'brien | Cure-dependent viscoelastic poisson's ratio of epoxy[END_REF] that the Poisson ratio increases during relaxation approaching values close to 0.5 for a fully relaxed material at high temperatures.

We therefore assume that the viscoelastic strain is isovolumetric, which corresponds to setting all K k in Eq.

(6) to zero (and, as a consequence from Eq. ( 4),

K 0 = K ∞ ).
Under this assumption, the stress relaxation at a constant tensile strain (with all other strain components unconstrained) is accompanied by an increasing transverse strain, thus reproducing the phenomena observed in [START_REF] Tschoegl | Poisson's ratio in linear viscoelasticity -A critical review[END_REF][START_REF] O'brien | Cure-dependent viscoelastic poisson's ratio of epoxy[END_REF]. This implies that, while the 1D model was identified under the assumption of a constant strain, this is no longer the case for the 3D strain tensor. We therefore re-identified the model parameters by simulating the experiments of Courtois et al. [START_REF] Courtois | Viscoelastic behavior of an epoxy resin during cure below the glass transition temperature: Characterization and modeling[END_REF] using the 3D model with the assumption of isovolumetric viscoelastic strain. This means that instead of identifying the tensile moduli E k of the Maxwell elements, we identify their shear moduli G k .

Courtois et al. [START_REF] Courtois | Viscoelastic behavior of an epoxy resin during cure below the glass transition temperature: Characterization and modeling[END_REF] proposed a continuous function to approximate the tensile viscous moduli as a function of the relaxation times τ k . In the present work, we use the same continuous function to describe the viscous shear moduli G k as a function of the relaxation times τ k :

G(τ k ) = A • exp - log 10 τ k -log 10 τ peak l peak + B 2 • 1 -erf log 10 τ k -log 10 τ peak l peak (7)
where τ peak , l peak , and A are the position, width, and height of a Gauss-curve like peak on the logarithmic relaxation time scale, and B is the height of the plateau of the function towards short relaxation times. erf is the Gaussian error function. The relaxation times τ k are given in minutes.

The temperature dependence of the viscoelastic behavior is modeled by means of the principle of timetemperature superposition. An Arrhenius model is used for the shift factors [START_REF] Simon | Modeling the evolution of the dynamic mechanical properties of a commercial epoxy during cure after gelation[END_REF][START_REF] Sadeghinia | Characterization and modeling the thermo-mechanical cure-dependent properties of epoxy molding compound[END_REF][START_REF] Courtois | Viscoelastic behavior of an epoxy resin during cure below the glass transition temperature: Characterization and modeling[END_REF]:

a T (T, c) = H(c) R • ln10 1 T - 1 T g (c) (8)
with the cure dependent T g as reference temperature. H is an activation energy, which was found to decrease linearly with the degree of cure [START_REF] Courtois | Viscoelastic behavior of an epoxy resin during cure below the glass transition temperature: Characterization and modeling[END_REF], and R the idealized gas constant.

The cure dependence is introduced via T g by means of the DiBenedetto equation [START_REF] Pascault | Relationships between glass transition temperature and conversion[END_REF] T

g (c) -T g0 T g1 -T g0 = λc 1 -c (1 -λ) (9)
where T g1 is the glass transition temperature of the fully cured resin, and T g0 is the glass transition temperature of the uncured resin. c designates the degree of cure.

The cure and temperature dependent relaxation times are thus given by:

log 10 τ (T, c) = log 10 τ (T g (c)) + a T (T, T g (c)) (10)
In this model, the instantaneous elastic modulus is independent of the degree of cure, as confirmed experimentally by [START_REF] Kim | Stress relaxation behavior of 3501-6 epoxy resin during cure[END_REF][START_REF] O'brien | Viscoelastic properties of an epoxy resin during cure[END_REF]. However, since the relaxation curves are shifted towards longer relaxation times for increasing degrees of cure, the apparent elastic modulus (slope of the stress-strain curve over a given strain range in the linear domain) increases with cure, and the Poisson ratio decreases from values close to 0.5 at the gel point to the values observed for the fully cured resin at room temperature, as observed by [START_REF] O'brien | Cure-dependent viscoelastic poisson's ratio of epoxy[END_REF].

The parameters of the model are summarized in Table 1. Note that they are different from those obtained by Courtois et al. [START_REF] Courtois | Viscoelastic behavior of an epoxy resin during cure below the glass transition temperature: Characterization and modeling[END_REF], because in the present work, we introduce viscoelasticity into the shear behavior (instead of the tensile behavior) to model the same experimental data, which implies that the 3D strain tensor is not constant during the relaxation phases. The model is numerically implemented with one Maxwell element per decade of relaxation times [START_REF] Courtois | Viscoelastic behavior of an epoxy resin during cure below the glass transition temperature: Characterization and modeling[END_REF], which is sufficient to well reproduce the observed relaxation curves. The relaxation times are chosen such that they reach from the highest tested strain rate at the coldest temperature (shortest shifted relaxation time) to the longest tested relaxation time at the highest temperature (longest shifted relaxation time).

Viscoelastic homogenization

Material behavior in the Laplace-Carson space

The Laplace-Carson (LC) transform of a function f in time t is defined as

f (p) = p ∞ 0 f (t)e -pt dt (11)
Applied to the integral form of a viscoelastic behavior (Eq. 1), a linear relation between the transformed stress and strain is obtained:

σ(p) = Ê(p) : ε(p) (12)
This means that in the transformed space, the viscoelastic constitutive equation corresponds to Hooke's law with stiffness tensors depending on the transformation parameter p. Therefore, common homogenization techniques for elastic materials can be used in the transformed space to calculate the LC transform of the homogenized viscoelastic behavior. To this result, the inverse LC transform has to be applied in order to determine the homogenized viscoelastic behavior in the time space.

The LC transform of a relaxation stiffness written in terms of a Prony series (Eq. 2) is given by

Ê(p) = E ∞ + N k=1 p p + 1 τ k E k (13)
For the special values of p = 0 and p = ∞, the relaxed ( Ê(0) = E ∞ ) and the instantaneous stiffness ( Ê(∞) = E 0 , see Eq. 4) are obtained. This means that the relaxed and the instantaneous stiffness tensor of the homogenized material can be calculated directly by elastic homogenization using the relaxed and instantaneous stiffness tensors of the constituents, respectively. The LC transform is only needed to calculate the E k of the viscoelastic Maxwell elements.

Homogenization strategy

The relaxation spectrum of the matrix can well be described by the continuous function given in Eq. ( 7). The function is sufficiently smooth that one relaxation time per decade is enough to well reproduce the experimentally observed relaxation curves [START_REF] Courtois | Viscoelastic behavior of an epoxy resin during cure below the glass transition temperature: Characterization and modeling[END_REF]. We suppose in the following that the yarns and the composite also have a continuous relaxation spectrum that can be described by one relaxation time per decade. In this case, the model described in section 2 can be used in its anisotropic form to model the yarns and the composite, and we choose the same relaxation times for the yarns and the composite as for the matrix. Consequently, the homogenized model will cover the same range of temperatures and strain rates as the matrix model. This approach is similar to those presented in [START_REF] Brenner | A "quasi-elastic" affine formulation for the homogenised behaviour of nonlinear viscoelastic polycrystals and composites[END_REF][START_REF] Ricaud | Effective properties of linear viscoelastic heterogeneous media: Internal variables formulation and extension to ageing behaviours[END_REF].

Lévesque et al. [START_REF] Lévesque | Numerical inversion of the Laplace-Carson transform applied to homogenization of randomly reinforced linear viscoelastic media[END_REF] discussed the accuracy of the numerical inversion of the LC transform and proposed an algorithm that determines the relaxation times of the homogenized material such that the accuracy of the inverse LC transform meets an acceptance criterion. According to Lévesque et al. [START_REF] Lévesque | Numerical inversion of the Laplace-Carson transform applied to homogenization of randomly reinforced linear viscoelastic media[END_REF], a priori definition of the relaxation times affects the quality of the inversion of the LC transform. However, if the relaxation spectrum for both the constituents and the composite is continuous and sufficiently smooth, the choice of the same, regularly spaced, relaxation times will provide an accurate reproduction of the homogenized relaxation curves, if the stiffness tensors of the Maxwell elements are appropriately determined. This is not obvious, since the underlying morphology of the multi-phase material may strongly modify the resultant spectrum shape [START_REF] Beurthey | Structural morphology and relaxation spectra of viscoelastic heterogeneous materials[END_REF].

We therefore control for the conditioning of the inverse LC transform and propose an optimum choice of the parameters p, at which the transformed stiffnesses should be calculated if the same regularly spaced relaxation times are used for the composite and all its constituents (section 3.3). To validate the hypothesis of continuous homogenized spectra and the accuracy of the proposed homogenization method, the homogenized viscoelastic behavior is compared with FE simulations of the viscoelastic behavior of the full heterogeneous RUC used for homogenization (section 3.4 for the micro-scale and section 4.2 for the meso-scale).

The proposed homogenization strategy is thus summarized as follows:

• Calculate the LC-transformed stiffness tensors for each constituent for M different values of p using Eq. ( 13)

• Calculate the homogenized LC-transformed stiffness tensors for each p by elastic homogenization

• Calculate the homogenized relaxed stiffness tensor E ∞ by elastic homogenization

• Calculate the homogenized E k by solving the equation system given by Eq. ( 13) using the homogenized Ê(p) and E ∞

The third point of the strategy is actually equivalent to the second point for p = 0, and the elastic homogenization of the instantaneous stiffness is equivalent to the second point for p = ∞. The last step corresponds to the numerical inversion of the LC transform. Its accuracy thus depends on the conditioning of the equation system Eq. ( 13). In order to calculate the N tensors E k , this equation system needs full rank. Therefore, the Ê(p) of each constituent have to be calculated for at least N different values of p, i.e., M ≥ N . If more p are used, the Moore-Penrose pseudoinverse of the transformation matrix is calculated to determine the homogenized E k .

Conditioning of the inverse Laplace-Carson transform

The conditioning of the inverse LC transform is analyzed for a homogenization starting from a micro-scale representative unit cell (RUC) of a yarn. At this scale, the RUC is relatively simple and thus computational costs are low. A hexagonal periodic distribution of parallel fibers was chosen as RUC for the yarns (Figure 1), because it yields a transverse isotropic elastic behavior. The fiber volume fraction in the RUC was fixed such that the overall fiber volume fraction and the warp to weft yarn weight ratio of the composite specimens studied in section 4.3 is obtained with the mesoscopic composite model used in section 4.1. In the following, the method is demonstrated for the case of a warp yarn (74.6% fiber volume fraction).

A FE mesh of the hexagonal unit cell is generated using second order prismatic elements (15 nodes) with reduced integration. The unit cell for a warp yarn shown in Figure 1 contains 3068 elements, 656 in the fiber and 2412 in the matrix. The matrix discretization was chosen such as to ensure at least 6 matrix elements between two neighboring fibers. Since the fiber is much stiffer than the matrix, the strain field in the fiber is relatively homogeneous. The element size in the fiber can thus be chosen much larger than in the matrix. Periodic boundary conditions (PBC) are applied in the three directions of periodicity in the plane perpendicular to the fiber, and in fiber direction. The FE problems are solved using an implicit scheme solved using the Newton-Raphson method. Convergence is reached if the norm of the internal forces is smaller than 10 -5 times the norm of the external forces.

For each elastic homogenization, the RUC is loaded such that its average strain corresponds successively to the 6 elementary strains (only one component ε kl = ε lk of the strain tensor is non-zero). The components of the homogenized elastic stiffness tensor E are then obtained from average stress response σ ij to the average strain ε kl as E ijkl = σ ij /ε kl . This method is used to calculate the homogenized transformed stiffness tensors Ê(p i ) for different choices of p i , i ∈ {1, . . . , M } with M ≥ N . For each p i , the transformed matrix stiffness is obtained from Eq. 13. Since the fiber behavior is elastic, its transformed stiffness is independent of p and equal to its stiffness in time space. The elastic properties of the fiber (transverse isotropic) are given in Table 2.

From the homogenized Ê(p i ), the parameters of the homogenized generalized Maxwell model are obtained by solving the equation system [START_REF] Russell | A new method to reduce cure-induced stresses in thermoset polymer composites, Part III: Correlating stress history to viscosity, degree of cure, and cure shrinkage[END_REF]. A measure for the precision of the solution is the condition number of the matrix

L ik = pi pi+ 1 τ k
. It only depends on the relaxation times τ k and the choice of the p i . Therefore, the results presented in the following for the case of a micro-scale homogenization also apply to the mesoscopic scale (section 4).

A convienient choice for one p i is p = ∞, because Ê(∞) corresponds to the instantaneous elastic stiffness (see Eq. 13). The homogenized Ê(∞) can therefore be obtained by simple elastic homogenization using the instantaneous stiffnesses of the constituents (see sections 3.1 and 3.2). This p = ∞ will be used in all of the following homogenizations. The other p i are chosen proportional to the inverse of the relaxation times, which depend on temperature T and degree of cure c through the shift factors a T (see Eq. 10). In this way, the p i depend on T and c through the negative shift factors -a T . As a consequence, the matrix L becomes independent of T and c. Since the stiffness tensors E k and E ∞ of the polymer resin are also independent of T and c (and the fiber behavior is independent of T and c as well), the homogenized E k and E ∞ are also independent of T and c. The proposed homogenization technique (in which the same relaxation times are used for the homogenized material and for its constituents) thus transfers the T and c dependence of the constituents directly to the homogenized behavior, with same shift factors for the relaxation times.

It is therefore sufficient to carry out the homogenization at the reference temperature T g (c). One relaxation time per decade is chosen such that the whole strain rate and temperature range of the experiments presented in section 4.3 is covered (τ k = 10 k min, k ∈ {-22, . . . , 2} at T g (c)). The simplest option ensuring that L is independent of T and c, is to choose p i = 1 τi in addition to p = ∞. However, the resulting homogenized spectrum exhibits some oscillations in the shear moduli (Figure 2), which for the case of the transverse shear modulus even includes some negative values. Adding one larger and one smaller p i (10 23 and 10 -3 ) has no notable influence on the results. If some information at intermediate times is added to the system by using p i = 10 i min -1 , i ∈ {∞, 22.5, 22.0, 21.5, . . . , -2.5}, the spectra become smooth and the condition number of the matrix L is reduced by more than a factor of 4. The conditioning becomes even better, if the p i are limited to times that lie in between the relaxation times (p i = 10 i min -1 , i ∈ {∞, 22.5, 21.5, 20.5, . . . , -2.5}). This choice of the p i ensures smooth spectra, and only a total of M = N + 3 elastic homogenizations are necessary to obtain the homogenized viscoelastic behavior.

Comparison between the homogenized behavior and FE simulations

The homogenized behavior obtained with this choice of p i is validated by comparison with FE simulations using the heterogeneous RUC shown in Figure 1 with the polymer behavior presented in section 2.2 for the matrix and the properties given in Table 2 for the fiber. The temperature is set to T g (c) (which yields a T = 0 from Eq. 8). Accordingly, the relaxation times are set to τ k = 10 k min, k ∈ {-22, . . . , 2} at T g (c).

An average strain of 0.1% is applied to the RUC within 10 -17 min and held constant for 10 4 min (in order to cover 21 orders of magnitude in time, a default range in the used FE software). The other average strain components are free. The average stress response of the RUC corresponds to the relaxation master curve over the simulated time scale. Tracing the free strain components in time allows to validate the evolution of the Poisson ratio during relaxation. The same FE mesh, PBC, and convergence criteria are used as in the homogenization procedure (section 3.3).

Figure 3 shows that the homogenized behavior is in very good agreement with the full viscoelastic FE simulation for all load cases. We also observe that due to the isovolumetric viscous strain of the matrix, the Poisson effect increases during relaxation. The agreement over the whole relaxation time range confirms that the assumption of a continuous viscoelastic spectrum of the homogenized yarn material (section 3.2) was correct.

Viscoelastic behavior of the composite

Representative Unit Cell

In this section, the homogenization method presented in section 3 is used to calculate the viscoelastic behavior of a composite material with a 3D woven interlock reinforcement. Building a FE mesh whose element faces are consistent with the internal material geometry (interfaces between yarns in contact and between yarns and matrix pockets) is a difficult task [START_REF] Grail | Consistent Finite Element mesh generation for meso-scale modeling of textile composites with preformed and compacted reinforcements[END_REF], in particular for highly complex 3D reinforcements.

On the other hand, voxelized meshes are appropriate to calculate homogenized elastic properties [START_REF] Doitrand | Comparison between voxel and consistent meso-scale models of woven composites[END_REF][START_REF] Isart | Internal geometric modelling of 3D woven composites: A comparison between different approaches[END_REF] if the voxel size is chosen sufficiently small such that the volume fractions of warp and weft yarns and the matrix pockets are accurately represented. Since the proposed homogenization method is based on elastic FE calculations only, we use voxel meshes in the following to avoid the highly complex meshing task.

Voxel meshes can be built directly from µ-CT images of the fibrous reinforcement in the composite [START_REF] Schneider | A meso-FE voxel model of an interlock woven composite[END_REF][START_REF] Hello | Numerical simulations of woven composite materials with voxel-FE models[END_REF][START_REF] Naouar | 3D composite reinforcement meso F.E. analyses based on X-ray computed tomography[END_REF]. By analyzing the gray scale levels of the µ-CT images and the orientations of the grayscale gradients, each voxel can be assigned to either a matrix pocket or a specific warp or weft yarn. The resulting voxel mesh of a meso-scale RUC of the studied composite is shown in Figure 4. It contains a total of 373932 voxels (141 in warp direction, 102 in weft direction, and 26 in out-of-plane direction). The effective volume of the voxels assigned to the warp yarns, the weft yarns, and the matrix is calculated. The fiber volume fractions in the warp and weft yarns that are used in the following calculations are adjusted such that the overall fiber volume fraction is 59% and the weight ratio of the fibers in the warp and weft yarns is 79/21.

With the voxel mesh shown in Figure 4, the fiber volume fraction in the warp yarns is 74.6% and in the weft yarns 79.6%. For each slice of the voxel mesh, the center of the cross-section of each yarn is determined.

These points are used to build the yarn center line by spline interpolation. The fiber orientation in each yarn slice is then set to the direction of the tangent to the yarn center line in the respective slice. The elastic properties assigned to the yarn voxel elements are obtained by rotating the stiffness tensor assigned to the yarn such that the axial behavior (L) is obtained in the local fiber direction.

The same homogenization method as presented in section 3.1 is used with the optimum choice of p i given in section 3.3. Each voxel is simulated by a linear hexahedral element (8 nodes) with full integration. The same convergence criterion for the implicit scheme is used as at the microscale. 3D periodic boundary conditions (PBC) are used at the meso-scale as well, even if the reinforcement geometry obtained from micro-CT is not perfectly periodic. However, mixed uniform boundary conditions (MUBC) as used by Hello et al. [START_REF] Hello | Numerical simulations of woven composite materials with voxel-FE models[END_REF] on similar materials lead to an overestimation of the homogenized stiffness, because constant normal displacements are imposed on the RUC boundaries for the tensile load cases, which prohibit relative displacements of the yarns and matrix pockets in these directions. With PBC, we expect a slight underestimation of the composite stiffness, because the opposite ends of yarns with not perfectly periodic paths on the RUC boundaries are not attached entirely to each other. Since the voxel mesh is periodic due to its regular structure, the application of 3D periodic boundaries is straightforward.

Numerical validation of the homogenized behavior

The obtained homogenized viscoelastic behavior is compared to full FE simulations of the viscoelastic behavior of the RUC in the same way as at the microscopic scale (section 3.4), using the same mesh, boundary conditions, and convergence criteria as for the homogenization procedure (section 4.1). The results are shown in Figure 5. Very good agreement is observed for the tensile and shear relaxation curves in all directions, and for the evolution of the Poisson ratios in time, which shows that the assumption of a continuous viscoelastic spectrum for the composite material (section 3.2) was correct. There are some small differences towards the end of the relaxation, which illustrate the limits of the presented method. In fact, if viscoelastic elements with very low weights are included, oscillations can be observed in the homogenized spectra around the corresponding relaxation times. Therefore, the continuous function (Eq. 7) used to calculate the shear moduli of the Maxwell elements in section 3 was truncated for very long relaxation times, when the weights become small (less than 0.1% of the instantaneous elastic modulus). This leads to a relatively abrupt end of some components of the viscoelastic spectrum of the yarn, as can be seen from the transverse shear spectrum in Figure 2 (this phenomenon is less important for the matrix and for the axial shear spectrum of the yarn). The longest relaxation times are therefore not sufficiently well represented in the homogenized behavior. Nevertheless, these differences are small in the studied case, and the longest relaxation times become only relevant above T g , which is not the most important temperature range for our applications (as above T g , the stresses are almost immediately relaxed, this temperature range does not contribute much to the generation of residual stresses). The elastic homogenizations in the LC space needed to calculate the homogenized viscoelastic behavior were also carried out using a mesh with voxels of half the size as in the mesh shown in Figure 4 (i.e., 8 times as many voxels). The results are practically the same as with the mesh shown in Figure 5. We therefore conclude that this voxel mesh is fine enough for our purposes.

Comparison with experimental observations

Composites specimens were manufactured by resin injection molding (RTM) using the interlock woven reinforcement shown in Figure 4 and the epoxy resin tested by Courtois et al. [START_REF] Courtois | Viscoelastic behavior of an epoxy resin during cure below the glass transition temperature: Characterization and modeling[END_REF], on which the model presented in section 2.2 was identified. Tensile tests were carried out according to the ASTM standard D3039 using a universal testing machine equipped with hydraulic grips and an oven. The specimen strain was measured using strain gauges, and the temperature controlled by thermocouples glued to the specimens.

Acoustic emission sensors were placed on the specimens to detect the damage thresholds. The apparent elastic properties of the specimens were measured at different temperatures reaching from -50 • C to 200 • C, by determining the secant modulus of the stress-strain curve between strains of 1 • 10 -3 and 4 • 10 -3 . It was ensured that the material is not damaged in the strain range, over which the elastic moduli are determined.

Only fully cured specimens were tested. The experimental validation of the homogenization procedure is thus limited to the temperature dependence.

The tensile moduli were measured in direction parallel to the warp (index 1) and the weft (index 2) yarns, respectively. The in-plane shear moduli (indices 12) were measured by Iosipescu shear tests based on the ASTM standard D5379, while the out-of-plane shear moduli (indices 23 and 31) where measured using V-notched rail shear tests based on the ASTM standard D7078. The shear strain in the specimens was measured using strain gauge rosettes on both surfaces of the specimens. The apparent shear moduli were taken as the secant moduli between shear strains of 1 • 10 -3 and 4 • 10 -3 . In each of the testing setups, the loading speed was adjusted such that a strain rate of 10 -3 s -1 was measured with the strain gauges.

The homogenized viscoelastic behavior shown in Figure 5 was used to calculate the evolution of the apparent elastic properties of the composite for strain rates of 10 -3 s and 10 -4 s. To this aim, the stress response of the homogenized model under a monotonically increasing uniaxial strain was calculated using the equations given in section 2.1. The apparent tensile moduli are determined from the resulting stress-strain curve in the same way as for the experimental data. Likewise, the apparent shear moduli are obtained by applying a monotonically increasing pure shear strain. The relaxation times of the homogenized viscoelastic model are calculated using Eq. 10 with the shift factors of the matrix (as described in Section 3.3), taking T g of the fully cured resin (see Table 1) as the reference temperature. In this way, the evolution of the apparent moduli with temperature can be calculated using the homogenized viscoelastic model of the composite.

The results are compared to the experimental observations in Figure 6. Due to confidentiality obligations, the numerical and experimental moduli are normalized by the tensile modulus in warp direction measured at 25 • C. The numerical results agree well with the experiments up to temperatures of about 120 • C, where the experimentally measured moduli start to decrease significantly due to the vicinity to the T g of the polymer matrix. On average the moduli are slightly underestimated up to 120 • C, which may be due to the use of periodic boundary conditions. However, other reasons are possible, like for example that the polymerization of the matrix may be influenced by the presence of the fibers (even if the same temperature cycle was applied to the composite and the resin specimens) or that the fiber volume fractions vary within the yarns (in our model, we assumed a constant fiber volume fraction everywhere in the warp yarns and a different constant fiber volume fraction everywhere in the weft yarns). Significant differences around the T g of the matrix are not surprising, since only very few experimental data was available to identify the model presented in section 2.2 at temperatures above 130 • C (see [START_REF] Courtois | Viscoelastic behavior of an epoxy resin during cure below the glass transition temperature: Characterization and modeling[END_REF]), and therefore, the matrix model is probably not very reliable in this temperature range. Above T g , the matrix becomes very soft, such that interactions between the yarns become important. These interactions are not taken into account in the meso-scale composite model, which is based on a small strain formulation. Therefore, the calculated relaxed moduli are lower than observed experimentally. In fact, above T g , the behavior of the composite is expected to become closer to that of the reinforcement without the matrix. On the whole, considering that we only used the properties of the constituents as input data (the viscoelastic model for the epoxy matrix identified independently, a voxelized FE mesh of the composite identified from µ-CT images, and elastic properties for the carbon fibers that were proven to give accurate results at room temperature), the proposed method provides an excellent estimation of the evolution of the elastic properties of a composite with a complex 3D woven reinforcement with temperature.

Conclusions

A homogenization technique was proposed exploiting the linear relationship between stress and strain in a linear viscoelastic model in Schapery's integral form [START_REF] Schapery | On the characterization of nonlinear viscoelastic materials[END_REF] after application of the LC transform. Under the assumption of continuous viscoelastic spectra of all constituents and the composite material, the same relaxation times were chosen for the constituents and the composite (one per decade covering the range of interest). The LC parameters p, at which the LC transform of the viscoelastic behavior of the constituents is calculated and the elastic homogenizations are carried out in the transformed space, were chosen at intermediate values with respect to the inverse of the relaxation times. This choice reduces the condition number of the equation system to be solved to calculate the homogenized viscoelastic behavior via the inverse of the LC transform and avoids oscillations in the homogenized spectra observed with other choices of p. In this case, the proposed homogenization technique consists of N + 3 elastic homogenizations in the transformed space, where N is the number of relaxation times of the constituents. Since in the matrix model, temperature T and degree of cure c only influence the relaxation times and the parameters p are based on the inverse of the relaxation times, the LC transform matrix is independent of T and c, and hence the same T and c dependence of the relaxation times is obtained for the homogenized behavior.

The proposed method was applied to the case of a composite made of a 3D interlock woven reinforcement and an epoxy matrix, whose viscoelastic behavior was determined experimentally elsewhere [START_REF] Courtois | Viscoelastic behavior of an epoxy resin during cure below the glass transition temperature: Characterization and modeling[END_REF]. First, a homogenization at the microscopic scale was carried out to obtain the viscoelastic behavior of the warp and weft yarns. In a second step, these behaviors were used together with the matrix behavior and a mesoscopic RUC obtained from micro-CT images of the composite to calculate the homogenized viscoelastic behavior of the composite. Numerical comparison of the obtained homogenized behavior and full scale FE simulations show excellent agreement over the whole relaxation time scale, both at the microscopic and the mesoscopic scale.

The resulting homogenized viscoelastic model was used to predict the temperature dependence of the mechanical behavior of the composite between -50 • C and 200 • C. Good agreement with experimentally tested composite moduli was observed, in particular in the temperature range where the matrix behavior was properly identified. This result is very promising, as some major complexities, such as a varying cure cycle of the resin within the composite compared to pure resin specimens (due to the presence of the reinforcement) or fiber volume fractions that vary within different yarns were neglected.

As the model of the resin also includes the dependence of the viscoelastic behavior on degree of cure, it should be possible to estimate the viscoelastic properties of partially cured composites with the proposed homogenization method. To validate this aspect, composite specimens manufactured via interrupted cure cycles will be tested mechanically in a future study. The fully validated procedure can then be used to model of the generation of residual stresses during the curing process in composites with 3D textile reinforcements.

Very promising results were already obtained by Benavente et al. [START_REF] Benavente | Numerical analysis of viscoelastic process-induced residual distortions during manufacturing and post-curing[END_REF] in using the presented homogenized model to predict the evolution of the distortion of composite plates and L-angles due to creep and relaxation phenomena linked to the viscoelastic nature of the matrix material. (continuous lines) and 0.0001s -1 (dashed lines) compared with experimentally measured moduli of the composite at a strain rate of 0.001s -1 (dots). The moduli are normalized by the experimentally measured tensile modulus in warp direction at 25 • C.

Figure 1 :i

 1 Figure 1: FE mesh of the micro-scale RUC (fiber in red, matrix in gray).

Figure 2 :

 2 Figure 2: (a) homogenized axial shear moduli (G LT (τ )) and (b) homogenized transverse shear moduli (G T T (τ )) of the Maxwell elements as a function of the relaxation time τ for the homogenized viscoelastic behavior of the warp yarns for different choices of p i = 10 i min -1 . The corresponding condition numbers of the LC transform matrix L are (in order of appearance): 1240.9, 1286.1, 299.3, and 221.0.

Figure 3 :

 3 Figure 3: Comparison between the stress response of the homogenized viscoelastic model for the warp yarns and the FE simulation using the RUC shown in Figure 1: stress in loading direction under constant tensile strain in axial (a) and transverse (b) direction with respect to the fiber, (c) shear stress evolution under constant shear strain in the axial (LT) and the transverse (TT) plane, and (d) evolution of the ratio between the negative transverse strain and the constant applied tensile strain in time axial (LT) and transverse (TT) direction.

Figure 6 :

 6 Figure 6: Secant moduli calculated using the homogenized viscoelastic behavior of the composite at a strain rate of 0.001s -1

Table 1 :

 1 Parameters of the viscoelastic model for the epoxy resin.E L [GPa] E T [GPa] G LT [GPa] ν LT ν T T

	A [MPa]	B [MPa]	log 10 τ peak	l peak	G ∞ [MPa]
	302.6	50.85		-0.3380	1.3638	2.089
	ν 0	H(c) [kJ/mol]	λ	T g0 [ • C]	T g1 [ • C]
	0.35	288.9 -111.4c	0.2603	155.0	46.78
		276.0	10.3	27.9	0.26 0.3

Table 2 :

 2 Transverse isotropic elastic properties of the fibers: E are Young's moduli, G the shear modulus, and ν the Poisson ratios. The index L indicates the direction parallel to the fiber, T indicates the transverse direction.

Acknowledgements

The research presented in this article was funded by the Safran Group, France.

 Figure 4: Voxel mesh of the meso-scale RUC of the composite (warp yarns in red, weft yarns in blue, matrix in gray).

Time [min] Shear stress [MPa]

FE calculation [START_REF] Billotte | Chemical shrinkage and thermomechanical characterization of an epoxy resin during cure by a novel in situ measurement method[END_REF] Homogenized behavior (