## Abstract

This study presents a family of estimates for effective potentials in nonlinear polycrystals. Noting that these potentials are given as averages, several quadrature formulae are investigated to express these integrals of nonlinear functions of local fields in terms of the moments of these fields. Two of these quadrature formulae reduce to known schemes, including a recent proposition (Ponte Castañeda 2015 *Proc. R. Soc. A* **471**, 20150665 (doi:10.1098/rspa.2015.0665)) obtained by completely different means. Other formulae are also reviewed that make use of statistical information on the fields beyond their first and second moments. These quadrature formulae are applied to the estimation of effective potentials in polycrystals governed by two potentials, by means of a reduced-order model proposed by the authors (non-uniform transformation field analysis). It is shown how the quadrature formulae improve on the tangent second-order approximation in porous crystals at high stress triaxiality. It is found that, in order to retrieve a satisfactory accuracy for highly nonlinear porous crystals under high stress triaxiality, a quadrature formula of higher order is required.

## 1. Introduction

A common engineering practice in the analysis of heterogeneous materials, either composite materials or polycrystals, is to use effective or homogenized material properties instead of taking into account all details of the geometrical arrangement of the individual phases. This has motivated the development of theoretical/analytical/computational homogenization techniques to predict the effective properties of such heterogeneous materials directly from those of the individual constituents and from their microstructure.

Variational principles play a central role in these approaches for several reasons. First, it is always easier to handle and to approximate scalar quantities, such as potentials, rather than tensorial quantities. Second, when a variational principle is available, rigorous bounds on the effective properties can be deduced by appropriately choosing the trial fields, delimiting a zone of potentially attainable properties in the space of effective parameters [1,2]. These variational principles are of current use for *linear composites* where, most of the time, the elastic moduli derive from a quadratic, convex, strain-energy function.

The situation in *nonlinear composites or polycrystals* is more complex. The behaviour of most common materials is partly reversible (part of the energy provided to the system is stored) and partly dissipative (part of the energy is converted into heat). However, the two extreme cases, purely (hyper)elastic materials and purely creeping materials corresponding to the behaviour of engineering materials under specific loading conditions, are of great interest. Both cases are usually governed by a single potential, strain-energy function in the former case and a viscous dissipation potential in the latter one. Then, the effective behaviour of the composite is again governed by a single effective potential which is the average of the minimum potential of a representative volume element (RVE) under prescribed overall strain or stress (this is the *average variational ‘principle’* of Whitham [3], which can be found in a different form in Ponte Castañeda & Suquet [4]). The effective response of the composite is again governed by a single variational principle. This property has led to improved bounds and estimates for nonlinear composites in the two simplified situations where either the energy dissipated or the energy stored in the deformation of the composite is neglected. Bounds for viscoplastic polycrystals improving on the Taylor bound were obtained by Dendievel *et al.* [5], who generalized the Hashin–Shtrikman variational principles to nonlinear media, and by Ponte Castañeda and co-workers ([6–8] and references therein), who introduced and used systematically the notion of a linear comparison polycrystal. In conclusion, the effective potential of a polycrystal, or an approximation of it, is obtained by computing, or approximating, the average of a nonlinear (and non-quadratic) function of the trial field which is then minimized. This average is an integral and computing or approximating this integral is the link between effective potentials and quadrature formulae (see §3).

When the composite consists of individual constituents governed by two potentials, the free energy and the dissipation potential, accounting for reversible and irreversible processes, respectively, the same need for determining two effective potentials arises. Note that, unlike the case of linear composites, the existence of potentials at the level of individual constituents is not guaranteed for all materials. It defines a subclass of materials called *generalized standard materials*, which covers most, but not all, elasto-viscoplastic or elasto-plastic materials. Restricting attention to this subclass, it can be shown [9–11] that the effective constitutive relations of the composite still derive from two effective potentials, the effective free energy, obtained as the minimum of the average of the local free energy when all state variables of the system are fixed, and a dissipation potential, which is the average of the local dissipation potential in an evolution of the state variables. In the situations that we have in mind, the potential that is hard to determine is not the one given by a minimization problem (in contrast with the single potential case), but rather the dissipation potential. Again, for this purpose one has to compute the average of a nonlinear (and non-quadratic in general) function of a field which links the determination of effective potentials to quadrature formulae.

This averaging is performed by an approach based on numerical integration. However, although the vocabulary used in the present article (numerical integration, Gauss quadrature, etc.) sounds very similar to that used in other works on reduced-order models ([12–14] among others), the present approach is quite different. Here, we work in the space of probability distributions of fields, and the quadrature formulae are approximations in this space. In particular, the integration points (Gauss points) in our approach are *not* spatial locations within the RVE, but are combinations of moments of the fields. The above-mentioned references deal with integration in *space* over a well-defined RVE, and their Gauss points (or integration points) are actual locations within the RVE.

The use of our integration rules at different orders is illustrated in the context of the reduced-order model non-uniform transformation field analysis (NTFA), initially proposed in Michel & Suquet [15,16] and further developed by several authors to approximate the exact homogenized model, and more specifically the two above-mentioned effective potentials. In full generality, the nonlinearity of the exact effective potential dissipation impedes a full decoupling between the microscopic and macroscopic scales. In earlier work [9,17] the present authors established an approximation of these evolution equations, achieving this decoupling by means of a linearization technique directly imported from nonlinear homogenization. It was based on the tangent-second-order (TSO) method, which performs a Taylor expansion of the potential in each phase in the neighbourhood of the average of the fields [9,17]. The resulting reduced model (NTFA-TSO) has been shown to be efficient in composites and polycrystals where all phases undergo incompressible plastic deformations [9,17]. However, the TSO linearization is known to be inaccurate for porous materials. The present study proposes a physically sound approximation for porous materials, whose accuracy depends on the material nonlinearity.

The paper is organized as follows. The role of effective potentials is recalled in §2 for crystal plasticity laws in the simplest form sufficient for our purpose here (without isotropic or kinematic hardening). The link with quadrature formulae (or numerical integration rules) is made in §3, where several quadrature formulae of different degrees of exactness are given. These quadrature rules are applied in §4 to estimate the effective potentials in nonlinear polycrystals. It is worth noting that one of these approximations has been previously derived by Ponte Castañeda [8] by completely different means, in the context of his ‘fully optimized’ variational method for constituents governed by a single potential. The last section is devoted to porous single crystals and shows how the quadrature formulae of different orders improve over the tangent approximation, especially at high stress triaxiality.

## 2. Crystal plasticity and effective potentials

### (a) Crystal (elasto-visco)plasticity

Polycrystals are aggregates of single crystals which differ by their orientation and shape but which obey the same constitutive relations, up to a rotation. The (visco)plastic strain in a single crystal results from slip on different slip systems
*M*_{e} is the fourth-order tensor of elastic compliance of the single crystal, *S* is the total number of slip systems, *m*_{s} is the Schmid tensor of system *s*, with slip plane normal *n*_{s} and slip direction *b*_{s} (orthogonal to *n*_{s}); the symbol : denotes contraction over two indices, and *γ*_{s} evolve. By *τ*_{s} on the same slip system and of other variables, which we will not take into account in the interest of simplicity. Without loss of generality, this relation can be written as
*ψ*_{s} may depend on the slip system. A typical example of potential corresponding to power-law creep is (the subscript (*s*) is omitted for simplicity)

### (b) One or two potentials

In creep problems, elasticity can be neglected and the total strain rate reduces to the viscoplastic strain rate. The constitutive relations at the single-crystal level are completely specified by only one potential *ψ*_{s} for each slip system,

However, when elastic and viscous deformations are of the same order (for cyclic loadings for instance), two potentials are required. The first potential is the free energy *w* (in the present case, the elastic energy) and the second potential is the viscous potential *ψ* (called force potential). We shall only consider the simplest case where

### (c) Polycrystals

A *polycrystal* is an aggregate of a large number *N* of single crystals with *different orientations* (grains). The size of the grains is small compared with that of the aggregate. Let *V* denote the volume occupied by an RVE of the polycrystal and *V* ^{(r)}, *r*=1,…,*N*, denote the domains occupied by the *N* grains. The characteristic functions of the domains *V* ^{(r)} are denoted by *χ*^{(r)}, and spatial averages over *V* and *V* ^{(r)} are denoted by 〈.〉 and 〈.〉^{(r)}, respectively; *c*^{(r)}=〈*χ*^{(r)}〉 is the volume fraction of phase *r* (grain *r*). The compact notation

All grains are assumed to be homogeneous, perfectly bonded at grain boundaries. Their slip systems *w*^{(r)} and *w* and *ψ* depend on the position ** x** inside

*V*,

#### (i) One potential

When there is only one potential, the overall strain rate

It is apparent from the expression of the local potential *ψ* and from the variational principle (2.7) that, in order to estimate accurately the effective potential

#### (ii) Two potentials

When the individual slips in the grains are governed by two potentials, it is still necessary to compute the average of the force potential 〈*ψ*(** σ**)〉, where

**is now the actual stress field with no immediate variational property.**

*σ*This can be further illustrated in the context of the NTFA, a reduced-order model initially proposed by Michel & Suquet [15]. It consists of two steps. In a first step, the viscous strain field *ε*_{v}(** x**) is decomposed on a finite basis of modes

*μ*

^{(k)}(

**),**

*x**k*=1,…,

*M*,

In a second step, evolution equations for the new unknowns of the model, the kinematic variables ** ξ** or the associated forces

**is the elastic strain localization tensor and**

*A*

*ρ*^{(k)}are elastic fields generated by the modes

*μ*^{(k)}. Unfortunately, the cost of evaluating each of the individual integrals

**(**

*σ***), followed by the evaluation of**

*x**ψ*(

**(**

*σ***)), which is a nonlinear local function, which is then averaged.**

*x*This has motivated the authors [9,17] to use techniques from nonlinear homogenization to approximate these integrals in terms of quantities which can be pre-computed. In these previous studies, the authors used the TSO approximation [19].

The present study is devoted to a family of other approximations, including approximations of higher order, with an unexpected connection with another recently proposed variational estimate [8]. The following section explores different ‘integration rules’ to evaluate integrals in the form *V* ^{(r)}.

It follows from (2.9) and (2.14) that, in crystalline and polycrystalline materials, the potentials *ψ* which are to be averaged are functions of the resolved shear stress *τ*_{s} on each slip system and are therefore functions of a *scalar field*. As will be seen with the help of probability distributions, the integration rules which are needed are approximations of integrals over a unidimensional variable. The extension of these integration rules to tensorial fields (with six components for instance) is not straightforward and would require numerical integration rules in a space of higher dimension (

## 3. Approximation of an integral by quadrature

For simplicity, the subscript *s* and the superscript *r* will be dropped in what follows.

Let *p*_{τ}(*t*) be the (unknown) probability density function (PDF) of the field *τ*(** x**) in

*V*, such that

The moments *E*_{n}(*τ*) of the field *τ*(** x**) and its central moments

*C*

_{n}(

*τ*) can be expressed as

Two of these moments are commonly used in nonlinear homogenization, under slightly different notations

The averaged potential can be rewritten as

Approximations of the last integral (3.4) using moments of *τ* of different orders are now derived, based on a procedure inspired from numerical integration (or quadrature), and more specifically Gauss quadrature, as recalled in appendix A. Note, however, that problem (3.4) has similarities to, but also differences from, that of numerical integration (A 4) described in appendix A. In numerical integration, the weight function *w* is known and the function *f* is not prescribed. Here, the ‘weight’ function *p*_{τ}(*t*) is unknown, whereas the potential *ψ* is given.

### (a) Quadrature formulae

It is considered here that the field *τ*(** x**) is given, as well as its PDF

*p*

_{τ}(

*t*). All moments,

*τ*

_{i}are the integration points and

*α*

_{i}are the weights. The approximation rule (3.5) is of

*interpolary type*because it only makes use of the values of

*ψ*at integration points. The integration points can be chosen arbitrarily or optimized. The

*quadrature formula*(3.5) has

*degree of exactness*if it is exact when

*q**ψ*is any polynomial of degree less than or equal to

*q*. Exactness to degree

*n*(at least) can be achieved for any choice of the

*τ*

_{i}’s by choosing

*l*

_{i}is the

*i*th characteristic Lagrange polynomial such that

*l*

_{i}(

*τ*

_{j})=

*δ*

_{ij}for

*i*,

*j*=0,…,

*n*. The interpolatory quadrature formulae (3.5) and (3.6) have degree of exactness

*n*for any choice of the nodes

*τ*

_{i}.

It is, however, possible to gain more accuracy (a higher degree of exactness *n*+*m*) by specific choices of the integration points. The optimization of the integration points is precisely the aim of Gauss quadrature in numerical integration and the maximum degree of exactness is 2*n*+1 corresponding to *m*=*n*+1. A similar procedure is followed here, although the optimality of the degree of exactness is not explored.

*One integration point, exactness to degree 1:* *n*=0, *m*=1. The quadrature formula (3.5) with a single integration point reads
*α*_{0} and *τ*_{0} have to be determined so that (3.7) is exact for *ψ* constant and linear. For this, (3.7) is written with *ψ*(*t*)=1 and

The optimal quadrature formula with one point is therefore

*Two integration points, exactness to degree 2:* *n*=1, *m*=1. The quadrature formulae with two points read as
*ψ* is a polynomial of degree less than or equal to 2. Defining new unknowns *h*_{i} as

Multiplying the second equation of (3.11) by *h*_{0}+*h*_{1} and using the first and third relations, we obtain that
*h*_{0} and *h*_{1} have opposite signs. Unfortunately, only three equations are available to determine the four unknowns, *α*_{i} and *h*_{i}; *α*_{0} can be chosen arbitrarily. The other three unknowns are determined as

Note that, for the *h*_{i} to be properly defined, *α*_{0} has to be chosen strictly between 0 and 1. In conclusion, with any choice of *α*_{0} between 0 and 1, the following two quadrature formulae are exact to order 2 (i.e. with equality whenever *ψ* is a polynomial of degree less than or equal to 2):

The two integration rules (3.14) and (3.15), which look different, yield in fact the same approximation for 〈*ψ*(*τ*)〉. This can be seen by interchanging *α*_{0} and 1−*α*_{0} in (3.13). In the example of §5, *α*_{0} was chosen equal to 1/2.

*Two integration points, exactness to third order:* *n*=1,*m*=2. The number of integration points is still *n*+1=2 but now the integration formula (3.5) is required to be exact for polynomials of degree 3, or less. Taking

This additional equation requires information on the third moment of *τ*. Then, combining (3.11) and (3.16), after multiplying the third equation in (3.11) by *h*_{0}+*h*_{1} and using (3.16), it is found that

Once again, there is an indeterminacy in the *h*_{i}’s, but it can be checked that the resulting integration rules coincide, so that there is no indeterminacy in the final approximation of 〈*ψ*(*τ*)〉. The resulting quadrature formula (3.18) is exact at order 3 (for all polynomials of degree less than or equal to 3). One degree of accuracy is gained by the knowledge of the third moment of *τ*. Note, in addition, that the *α*_{i}’s and the *h*_{i}’s which were left undetermined in the preceding paragraph are now completely determined.

### Remark

When *C*_{3}=0 (which is in particular the case when the PDF of the field *τ*(** x**) is centred around its mean), then

### (b) Higher-order approximations

Approximations of higher order can also be derived using the same procedure. Obviously, increasing the number of integration points will increase the degree of accuracy of the numerical integration formula (3.5), but will also require additional information about higher-order moments of the field *τ*(** x**). With

*n*+1 integration points

*n*+1), the following system of equations is obtained by enforcing equality in (3.5) with

*C*

_{0}=1). In principle, these 2

*n*+2 equations should permit determination of 2

*n*+2 unknowns

*α*

_{i}and

*h*

_{i}. However, the algebraic equations for the

*h*

_{i}’s are of high degree and it is not clear, in general, how the systems (3.20) should be solved. This is left for future work.

### (c) ‘Square convex’ potentials

It is often the case that the potential *ψ* is a function of the square of its argument. In other words, *ψ* can be written as

If, in addition, *f* is a convex function, the potential *ψ* is said to be ‘square convex’. This is, in particular, the case for power-law materials, with *n*≥1 to ensure convexity of *f*.

When the ‘square’ condition (3.21) holds, the different approximation rules (3.9), (3.14) and (3.18) can be applied to the potential *ψ* itself, as well as to *f*. The corresponding approximations for the averaged potential are different; they make use of different moments of *τ* and have different degrees of exactness. For instance, the approximation rule (3.18) for 〈*ψ*(*τ*)〉 makes use of the first three moments of *τ* and is exact whenever *ψ* is a polynomial of degree less than or equal to 3. The corresponding rule for 〈 *f*(*τ*^{2})〉 makes use of the first three moments of *τ*^{2} and is exact whenever *f* is a polynomial of degree less than or equal to 3, i.e. roughly speaking, when *ψ* is an even polynomial of degree less than or equal to 6.

The quadrature formulae for 〈 *f*〉 are easily deduced from (3.5).

*One integration point, exactness to degree 1 in* *f*,

*Two integration points, exactness to degree 2 in* *f*: *α* being arbitrary between 0 and 1,

Recall that there is also another solution with a change of sign in front of the square root in the expressions of the *h*_{i}’s. It should also be noted that nothing prevents *θ*_{0} from being negative. This occurs when the fluctuations of *τ*^{2} are larger than its average. Then, *f*(*x*) should be extended to negative *x*. This extension will be given for power-law materials in §5.

*Two integration points, exactness to degree 3 in* *f*: the integration rule is again (3.23) and (3.24) with

In (3.24) and (3.25), *C*_{2}(*τ*^{2}) and *C*_{3}(*τ*^{2}) are the centred moments of *τ*^{2} of degree 2 and 3, respectively. They can be expressed in terms of the moments of *τ* as

The quadrature formulae (3.23) and (3.24) (with *α*_{0} being arbitrary between 0 and 1) involves the moments of *τ*(** x**) up to order 4 and will be denoted as 〈

*ψ*

_{QF4}(

*τ*)〉. Similarly, the approximation (3.23) and (3.24) with (3.25) involves the moments of

*τ*(

**) up to order 6 and will be denoted as 〈**

*x**ψ*

_{QF6}(

*τ*)〉.

## 4. Approximate effective potentials in polycrystals

We come back to the estimation of effective potentials in polycrystals. The approximation rules of §3 with different degree of exactness are now reviewed.

### (a) Integration rule for 〈*ψ*(*τ*)〉 with one integration point

The simplest integration rule with one integration point to evaluate 〈*ψ*(*τ*)〉 is

This crude approximation will not be pursued here.

### (b) Integration rule for 〈 *f*(*τ*^{2})〉 with one integration point

Under assumption (3.21), approximation (3.9) written for *f* and not *ψ* with a single integration point reads as

The average potential *for an arbitrary stress field σ* is therefore approximated by

When the constitutive relations for each slip system at the single-crystal level derive from a single potential, the approximation (4.3) yields an approximation for the effective potential

The Euler–Lagrange equations associated with the variational problem at the right-hand side of (4.4) can be derived following a procedure identical to the one used in [4] and lead to a linear elasticity problem for a linear comparison composite (LCC). When the *ψ*_{s} are square convex, the right-hand side of (4.4) is actually an upper bound for the effective potential and the corresponding LCC is the optimal LCC, optimal in the sense of the variational theory of [6]. It should however be noted that the approximation rule (4.3) (which is an upper bound when the *ψ*_{s} are square convex) can be applied to *any* stress field ** σ**.

### (c) Integration rule for 〈*ψ*(*τ*)〉 with two integration points: link with the ‘fully optimized’ method

We consider the integration rule (3.14) with two integration points and exactness to degree 2. To simplify notations and to make the link with the ‘fully optimized’ variational estimate of Ponte Castañeda [8], we adopt the following notations: *α*_{0}, *τ*_{0} and *τ*_{1} in grain *r* and on system *s* are, respectively, denoted as ** σ** (not necessarily belonging to

**at this degree of exactness), and**

*σ*This approximation was first derived by Ponte Castañeda [8] in the context of materials with a single potential by completely different means. In addition, it was given for the solution of a specific variational problem, whereas it is given (and used) here for an arbitrary stress field ** σ**.

It is possible to make the link between the approximation rule (4.5) and the LCC of [8]. Indeed, when the approximate potential 〈*ψ*_{FO}(** σ**)〉 of (4.5) is substituted into the variational problem (2.7), the Euler–Lagrange equation deriving from this approximate potential reads as

**(**

*σ***) denote Frechet derivatives with respect to the field**

*x***. More specifically,**

*σ*Define

Then, the constitutive relations contained in the Euler–Lagrange equations associated with the approximate effective potential read as (in addition to equilibrium of ** σ** and compatibility of

**)**

*ε*These equations coincide with those derived by Ponte Castañeda [8] for his optimal linear comparison polycrystal.

### (d) Integration rules for 〈 *f*(*τ*^{2})〉 with two integration points

The two-point integration rules (3.23) and (3.24) and (3.23)–(3.25) for *f*(*τ*^{2}) (instead of *ψ*(*τ*)) can also be used to generate approximations of the effective potentials for arbitrary stress fields. The corresponding integration rules make use of the second and fourth moment of *τ* when the degree of exactness for *f* is 2, and of the second, fourth and sixth moment of *τ* when the degree of exactness for *f* is 3. These integration rules applied to crystals with a single potential will therefore lead to comparison composites which are no longer linear. However, they are useful in the context of the NTFA for crystals with two potentials.

The first integration rule corresponds to (3.23) and (3.24) applied to each system and reads as
** σ** at this degree of exactness), and

*QF*4 in the approximate potential refers to the fact that it makes use of the moments of the field up to order 4. A similar approximation 〈

*ψ*

_{QF6}(

**)〉 making use of moments up to order 6 is obtained by applying the integration rule (3.23)–(3.25) to**

*σ**f*for each slip system and each individual grain.

## 5. Sample example: non-uniform transformation field analysis for porous single crystals

This section illustrates the above approximation rules for crystals governed by two potentials in the context of the NTFA. In previous works [9,17], the authors proposed an approximation of the effective potential

The approximations of §4 should improve significantly on the TSO for porous materials. This is why porous single crystals are considered. But in contrast with [8], this is done in the context of the NTFA, allowing us to illustrate the hierarchy between the different integration rules of §4, and making use of the different moments of the stress field.

### (a) Porous single crystals

A volume *V* of a single crystal is considered. This volume contains 50 spherical, equi-sized, non-overlapping pores which are randomly distributed throughout *V* (with periodicity conditions when the pores intersect the boundary of *V*). The total porosity is *f*=9%.

The matrix material is a single crystal with hexagonal symmetry with 12 slip systems: three basal, three prismatic and six pyramidal (see [9] for details). The *c*-axis of the hexagonal crystalline structure is aligned with *e*_{3}. The elastic stiffness matrix reads as (all moduli in MPa)

The viscous slips on the different systems are governed by the power-law relation (2.3), with the same exponent *n*, the same reference slip rate *τ*_{0} on all systems. Three different sets of material parameters are considered,

### (b) Training paths: elasto-viscoplastic regime

In accordance with the general procedure described in detail in Michel and Suquet [9,16,17], the modes *μ*^{(k)} required by the NTFA are generated by the so-called snapshot proper orthogonal decomposition (POD) exploiting full-field simulations along specific loading paths called *training paths*. First, a parametric study of the influence of the mesh refinement has been performed. In all simulations, quadratic 10-node tetrahedral elements (with four integration points/element) were used. The meshes were generated using the open source platform SALOME (www.salome-platform.org). Figure 1*a* shows the dependence of the asymptotic overall stress on the mesh refinement in the most discriminating case of hydrostatic loading with the largest nonlinearity exponent *n*=10. Figure 1*b* shows the mesh selected in the subsequent full-field simulations (second-last most refined mesh with a total of 31 319 nodes). This discretization has been found to be sufficient to obtain accurate results in the most nonlinear case, *n*=10. Further refinement of the mesh did not lead to significant differences in the overall response (figure 1*a*).

The macroscopic stress is imposed along a fixed direction in stress space,

To avoid over-shooting of the asymptotic stress, the ‘control parameter’ *Σ*^{0} are outputs of the computation. The directions of overall stress *Σ*^{0} correspond to axisymmetric stress states
^{−3} s^{−1}.

Full-field simulations are performed for three different values of the power-law exponent, *n*=1,3 and 10. Twenty-five snapshots are stored along each training path and used to generate modes by means of the POD (see [16] for more details). The threshold used in the POD to select the modes is 5×10^{−5}. Two training paths were first used with stress triaxiality ratios *X*=0 (axisymmetric shear) and *n*) and two or three modes along the path of hydrostatic tension (two modes when *n*=1, three modes when *n*=3 and *n*=10). For reasons which will be given below, an additional training path was added for nonlinear materials (*n*=3,10) along a direction corresponding to a stress triaxiality ratio *X*=cot(*π*/10)≃3. The POD algorithm selected two modes along this path for both values of *n*. The first modes obtained for the three training paths (*n*=10) with stress triaxiality ratios *X*=cot(*π*/10), *X*=0, respectively, are shown in figure 2 for the purpose of illustration.

Once the modes for each training path are selected, it is first checked that the reduced-order model reproduces accurately the full-field simulations along these paths. This is shown in figure 3 for hydrostatic loading, which is the most discriminant loading case. The reduced-order models are referenced as follows: NTFA-TSO refers to the prediction of the previous model of [17] based on the TSO linearization. NTFA-FO refers to the approximation (4.5) with *QF*4 refers to the approximation (4.11) with all the *QF*6 refers to the approximation (4.12) where the *f* corresponding to the power-law potential (2.3) was adopted. Let *E*(*n*) denote the integer part of *n* and Frac(*n*)=*n*−*E*(*n*); then *f* is defined as

The modes are identical in all NTFA models. The following comments can be made.

— When

*n*=1, all four NTFA models (TSO, FO,*QF*4,*QF*6) coincide and give the exact average potential for a given stress field. As can be seen in figure 3*a*, their common prediction is in very good agreement with the full-field simulations, which shows that the decomposition (2.10) on a basis of modes is pertinent.— Under hydrostatic loading, and for rate-sensitivity exponents different from

*n*=1, the TSO linearization predicts a linear elastic response of the porous crystal. As is well known, this unrealistic prediction is due to the fact that the tangent operator of the potential*ψ*, computed at the average stress in the matrix which is hydrostatic, is infinitely stiff when*n*≠1.— The three reduced-order models (FO,

*QF*4,*QF*6) predict a non-elastic response for the porous crystal and therefore improve on the prediction by the TSO. The FO rule (4.5) is however not very accurate even for a moderate nonlinearity exponent*n*=3. The predictions*QF*4 and*QF*6 are quite similar for moderate nonlinearity but*QF*6 is clearly more accurate when*n*=10. The higher the information on the moments, the more accurate are the predictions.

### (c) Validation paths: gauge surfaces

Once the modes are identified and once it has been checked that the reduced-order models reproduce correctly the overall stress–strain response of the porous crystal along the training paths, it remains to assess the accuracy of the predictions along paths which were not used to identify the modes. These other paths are *validation paths*.

The accuracy of the reduced-order models is assessed by exploring loading conditions different from the training paths. Attention is focused on the purely viscous response of the porous crystal characterized in stress space by the *gauge surface* of the porous material. This surface is defined as [20]

This surface can be found by first exploring radial paths in stress space at different triaxiality ratios. When the strain along the loading direction *Σ*^{0} is large enough, the overall stress reaches a plateau *Σ*^{0} and lying on the gauge surface is determined by solving

Owing to the homogeneity of *ψ*), it is found that

The average potential

The gauge surfaces predicted by the different schemes are compared with those obtained by full-field simulations. The gauge surfaces are plotted here in the two-dimensional space of axisymmetric loadings (5.3) and the classical invariants *Σ*_{m} and *Σ*_{eq} are used (even though the full gauge surface is not isotropic in the six-dimensional space of overall stresses). First, the reduced-order models were implemented with modes determined with only two training paths corresponding to hydrostatic (*X*=0) states. Depending on the value of *n*, four or five modes were selected by the POD (four modes when *n*=1, five modes when *n*=3 and *n*=10). The results are shown in figure 4. When *n*=1, all reduced-order models agree and are in excellent agreement with the full-field simulations. When *n*=3,10, the predictions of the reduced-order models differ and the gap between the predictions increases with the stress triaxiality ratio, being maximal for hydrostatic loading (*QF*6.

However, when *n*=10, a gap can still be observed at intermediate stress triaxialities. This is why a third training path was used corresponding to a stress triaxiality ratio *X*=cot(*π*/10). Then, the POD added two modes for this path which brought the number of modes for *n*=3 and *n*=10 to 7. The predictions of the reduced-order models with seven modes are shown in figure 5. The gap with the full-field simulations at all stress triaxialities is reduced, at least with the integration rule based on the sixth moment of the resolved shear stress.

### (d) Speed-up

The CPU times for the different methods (full-field simulations, NTFA-TSO, NTFA-FO and NTFA-*QF* of different degrees of exactness) are compared in table 1 for purely deviatoric loading (the CPU times for the other loading paths are similar). The different speed-ups, measured by the CPU ratios given in table 1, can be commented on as follows.

— The fastest NTFA models are based on the TSO and the FO approximations with a speed-up of more than 10

^{6}. Regarding the TSO, as already commented upon, its predictions for hydrostatic loading are physically wrong and are inaccurate for all loadings with a significant amount of hydrostatic stress.— Among the three NTFA models based on integration rules, the fastest is clearly the FO model, which requires only the computation of moments up to order 2 of the resolved shear stress on all systems. The speed-up is between 10

^{6}and 2×10^{6}depending on the number of modes. The speed-up for*QF*4, which requires moments up to order 4, varies between 1.5×10^{4}and 3×10^{4}. The speed-up for*QF*6, which requires moments up to order 6, varies between 30 and 100. The price paid to increase the accuracy by introducing higher-order moments is quite significant. The different moments required by the different schemes are recalled in appendix B. The number of operations required by FO scales as (6+*M*)^{2}, whereas it scales as (6+*M*)^{4}for*QF*4 and (6+*M*)^{6}for*QF*6, where*M*denotes the number of modes in the decomposition (2.10).*QF*4 should therefore be (6+*M*)^{2}times slower than FO, which is consistent with the CPU times reported in table 1 (about 100 times slower for five or seven modes). Regarding*QF*6, it involves the computation of the*α*_{i}’s, in addition to the*h*_{i}’s. Therefore,*QF*6 is more than (6+*M*)^{2}slower than*QF*4 (300 times slower rather than 100 times).— The number of modes also has a significant impact on the number of operations. Referring again to the formulae of appendix B, the number of operations between

*M*_{1}modes and*M*_{2}modes scales as*r*^{2},*r*^{4},*r*^{6}, with*r*=(6+*M*_{2})/(6+*M*_{1}). For instance, when*M*_{1}=5 and*M*_{2}=7, then*r*≃1.18. So the FO,*QF*4 and*QF*6 with seven modes are about 1.4, 2 and 3 times slower, respectively, than the same schemes with five modes.

## 6. Conclusion

The conclusion of this study is twofold.

— First, quadrature formulae offer a different perspective on the recently proposed ‘fully optimized’ method of Ponte Castañeda [8] to estimate effective potentials in nonlinear composites. These quadrature rules are also helpful to generate a family of estimates for effective potentials generalizing this previous estimate. In particular, statistical information of high order on the local fields, when it is available, can be incorporated in these new approximations.

— Second, it has been shown how formulae at different orders remove two limitations of the reduced-order model based on the TSO approximation [9]. The TSO approximation was found to be accurate for incompressible polycrystals with moderate nonlinearities, but its accuracy deteriorates when the nonlinearity becomes stronger. Another well-known limitation of the TSO approximation is that it fails for porous crystals subject to highly triaxial stress states (this is because the linearization is made about the average stress in each crystal). These two limitations are partially removed in the present study. Not surprisingly, it is found that, in order to retrieve a satisfactory accuracy for highly nonlinear porous crystals under high stress triaxiality, a quadrature formula of higher order is required. Considering higher-order moments of the fields comes with a computational cost and therefore a compromise must be found between cost and accuracy.

## Data accessibility

This article has no additional data.

## Authors' contributions

P.S. was more specifically in charge of the theoretical developments on quadrature formulae. The numerical simulations were performed by J.-C.M. Both authors participated in the reduced-order model and in the preparation of this manuscript. They gave their final approval for publication.

## Competing interests

We declare we have no competing interests.

## Funding

The authors acknowledge the support of the Labex MEC and of A*Midex through grant nos. ANR-11-LABX-0092 and ANR-11-IDEX-0001-02

## Acknowledgements

The authors acknowledge fruitful discussions with M. Idiart.

## Appendix A. Reminder on numerical integration

Consider the weighted integral
*w* is the weight. We see that *I*(*f*) can be approximated by a quadrature formula

The above equation is said to have *degree of exactness r if I_{n}(f)=I(f) for any * where

*r*. A natural approximation of

*I*(

*f*) is

*Π*

_{n}(

*f*) is the interpolating Lagrange polynomial of

*f*over a set of

*n*+1 distinct nodes {

*x*

_{i}},

*i*=0,1,…,

*n*, and

*i*th characteristic Lagrange polynomial such that

*l*

_{i}(

*x*

_{j})=

*δ*

_{ij}for

*i*,

*j*=0,…,

*n*. As

*Π*

_{n}(

*f*)=

*f*for any

*n*for any choice of the nodes

*x*

_{i}. It could however be possible that more accuracy (a higher degree of exactness

*n*+

*m*) is gained by specific choices of nodes. The aim of Gauss quadrature is precisely to optimize the choice of these nodes [21], theorem 10.1.

*For a given* *m*>0, *the quadrature formula*
*has degree of exactness* *n*+*m* *iff it is of interpolatory type and the nodal polynomial* *ω*_{n+1} *associated with the nodes* *x*_{i} *is such that*
*where* *denotes the set of polynomials of degree* ≤*m*−1 *and*

The maximum degree of exactness that one can get is 2*n*+1 (corresponding to *m*=*n*+1).

## Appendix B. Moments

According to (2.15), the resolved shear stress *τ*_{s} on system (*s*) reads as

Then

The following elementary moments are therefore computed once and for all and stored:

## Footnotes

1 Other boundary conditions can be considered [11].

- Received March 24, 2017.
- Accepted August 1, 2017.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.