Lei Zhang, Shaoqiang Tang

HEDPS and LTCS, College of Engineering, Peking University, Beijing 100871, China

Keywords: Homogenization Self-consistent clustering analysis Lippmann-Schwinger equation Green’s function

ABSTRACT Recently proposed clustering-based methods provide an efficient way for homogenizing heterogeneous materials, yet without concerning the detailed distribution of the mechanical responses.With coarse fields of the clustering-based methods as an initial guess, we develop an iteration strategy to fastly and accurately resolve the displacement, strain and stress based on the Lippmann-Schwinger equation, thereby benefiting the local mechanical analysis such as the detection of the stress concentration.From a simple elastic case, we explore the convergence of the method and give an instruction for the selection of the reference material.Numerical tests show the efficiency and fast convergence of the reconstruction method in both elastic and hyper-elastic materials.?2021 The Author(s).Published by Elsevier Ltd on behalf of The Chinese Society of Theoretical and Applied Mechanics.

Macroscopic materials are composed of microscopic fine struc- tures, which are often heterogeneous.To predict quickly and faith- fully the macroscopic behavior based on microscopic structures, re- duced order models, such as the transformation field analysis (TFA) [1] , the nonuniform transformation field analysis (NTFA) [ 2 , 3 ] and proper orthogonal decomposition (POD) [4] have been proposed.

More recently, data-driven clustering-based homogenization methods, such as self-consistent clustering analysis (SCA) [5-10] , virtual clustering analysis (VCA) [11-13] and finite element method (FEM) clustering analysis (FCA) [14] , were developed to further im- prove the efficiency.These methods decompose the microscopic computational domain into several sub-regions called as clusters by machine learning techniques according to the strain similar- ity in pre-computed high-fidelity results.Assuming uniform strain and stress in each cluster during a loading/unloading process, clustering-based homogenization is able to provide the macro- scopic/average properties efficiently, by ignoring the detailed dis- tribution of the responses.While demonstrated effective in many applications including fractue, it may fail.

On the other hand, when detailed distributions are demanded, one may adopt finite element or meshfree methods.Unfortunately for heterogeneous materials, they usually induce immense compu- tational cost.To reduce the cost, for periodic cells, some numeri- cal methods [15-17] make use of the fast Fourier transform (FFT) to achieve higher efficiency both in terms of speed and memory.In general, these methods take zero or uniform fields for initial- ization.As a matter of fact, clustering-based methods provide a roughly correct distribution of the mechanical responses in form of piece-wise constant fields, hence serve as a good initial guess for faster convergence in iterations.When dealing with history- dependent cases, classical numerical methods might have to sim- ulate the whole loading process in an incremental way.Mean- while, clustering-based methods might also provide rough infor- mation such as state variables in the elasto-plastic materials.For the sake of efficiency, we only need homogenization rather than local detailed features in most case.The clustering-based results might help to determine the location and loading for refinement.

In this paper, we propose an iteration strategy based on the Lippmann-Schwinger equation with coarse fields of the clustering- based homogenization as the initial guess.This scheme accu- rately resolves the displacement, strain and stress, and recovers the lost distribution details of reduced-order model.We take periodic boundary conditions for illustration and use SCA results as the ba- sis.Based on the periodic Lippmann-Schwinger equation, we re- solve the strain and stress fields with iteration scheme first.Then we reconstruct the displacement field using the Green’s method.Furthermore, we explore the relationship between the convergence and the selection of the reference material.Numerical examples show that the iteration scheme converges fast when the appropri- ate reference material is chosen.

In the following, we first present a general framework for the iteration scheme, and then introduce the problem setup and the Lippmann-Schwinger equation.We briefly review the SCA method.The reconstruction of strain and displacement fields is proposed later.We then study the convergence and give an instruction for the selection of the reference material.Numerical examples verify the efficiency of the method.Some remarks are drawn finally.

Mathematical framework

Considerageneralproblemdescribedby

whereuis the unknown function of independent variablesx.Lis a linear operator andNis a nonlinear operator.Usually by some transformation such as the Gauss/Stokes theorem, we might find an adjoint operatorL*satisfying

where*denotes convolutionv*w()d,Bis a linear operator and?Ωis a boundary integral on?Ω.We then obtain the Green functionφby

Substitutev=u,w=φinto Eq.(2) , we obtain the following in- tegral equation,

This actually gives the Green’s method.Considering the unique- ness of the solution to Eq.(1) , Eq.(4) provides a way to solve the problem by iterations, i.e.,

If it converges, the limit gives the fixed point, i.e., the solution to Eq.(1) .

Methodology and the Lippmann-Schwinger equation

Consider the material in a representative volume element (RVE)Ω.The local problem inΩis governed by equilibrium equation, strain-displacement equation and constitutive equation

whereσ=σ(ε;x)is a function of the strainεand the positionx.For the sake of clarity, denote the operator.In most applications, RVE consists of a limited number of phases and there are different local constitutive relations in each phase.

In order to close the problem, we choose periodic boundary conditions as follows.The whole strain fieldε(x)is split into its averageEand a fluctuation field(x)

Correspondingly, the displacement fieldu(x)consists of two terms

where(x)is the fluctuation term corresponding to(x).Periodic boundary conditions refer to periodic(x)and hence, and anti- periodic tractionn·σ.

This problem with periodic boundary conditions can be solved by FEM [14] .In this paper, we use clustering-based methods, illus- trated with SCA, to solve it efficiently, while the result is merely a rough piece-wise constant strain field.Based on this, we propose a method to refine the strain field and reconstruct the displacement field.

We rewrite Eq.(6) as

Here,C0denotes the stiffness of a reference homogeneous linear elastic material andτ=σ?C0:εis called as the polarization ten- sor.Here, corresponding to Eq.(1) , let

By virtue of the Green’s function, the solution of Eq.(11) can be expressed by applying periodic Green functionφ(i)(x?), namely, the solution to the problem with a concentrated force alongxidi- rection at the point, associated withC0as

with

It is noticed that this equation satisfies the compatibility condi- tions directly.The explicit form in Fourier space is

forΦ, and

forΨ, whereξis the coordinate in Fourier space corresponding toxin real space,λ0andμ0are the Lamécoefficients of the refer- ence stiffnessC0.The formulation works for both 2D and 3D prob- lems.The convolution in both Eqs.(13) and (15) can be translated into a direct multiplication in Fourier space:

which can be realized by applying FFT algorithm.

Review of self-consistent clustering analysis

SCA [5] is a two-stage reduced order modeling approach con- sisting of an offline data compression stage and an online predic- tion stage.In the offline stage, the original high-fidelity RVE rep- resented by voxels or elements is compressed into clusters.In the online stage, macroscopic loading is applied to the reduced order (clustered) RVE, i.e., average strainEis given.SCA works for peri- odic problems described previously.

The key idea is to decompose the computational domain intoknon-overlapped sub-regions, so-called clusters, denoted byΩI,I= 1,2,···,kand assume uniform strain and stress in each cluster, i.e.,

distinguished mathematically by the characteristic functionχ

These clusters are obtained by machine learning technique ac- cording to material responses under limited loading cases dur- ing the offline stage.This allows one to discretize the Lippmann- Schwinger equation into

Finally, we solve the above algebraic system Eq.(23) and derive the macroscopic constitutive relation between a macroscopic stress and strain.

The flowchart is summarized as below:

Offline:

(a) Under the selected loadings, compute the mechanical re- sponses with high-fidelity numerical simulations and restore data information in tensorA.

(b) DecomposeΩintokclusters (Ω1,Ω2,···,Ωk) by a cluster- ing algorithm according toA.

(c) ComputeDIJ.

Online:

(a) At the current loading step, solve the Lippmann-Schwinger equation (23) and obtain the strainεI(I= 1,2,···,k)in each clus- ter.

(b) Compute macroscopic/average strain and stress.

Reconstruction method

Clustering-based methods such as SCA solve the homogeniza- tion problems efficiently, but they only give a rough piece-wise constant strain field, and do not offer local information correctly.We propose a reconstruction method to refine the results based on the continuous Lippmann-Schwinger equation.Based on the re- fined strain field, we can also construct the displacement field ef- ficiently according to Eq.(13) .

Refinement of the strain field

We consider the elastic and inelastic cases, with a microscopic constitutive relation

Take the rough results of SCA as the initial guess and give a recursive scheme according to Eq.(15) :

wherendenotes the iteration step andΦis the kernel function associated with the reference modulusC0used in the refinement.

The convolution is calculated in the Fourier space readily, so the algorithm is implemented with FFT algorithm.In general, the it- erations converge very fast when selecting appropriate reference material.In the following numerical examples, the refined results reach a good accuracy in about 3 iteration steps.Remark that the kernelΦhas been computed in the offline stage of SCA, so we do not need to compute it again, which improves the efficiency.

We shall discuss the relationship between the convergence and the reference modulusC0in a while.Note that if the scheme is convergent, the final result converges to the analytical solution of the Lippmann-Schwinger equation, i.e., the real solution of the system.The refined resultsεrefine(x)are continuous, give specific strain response and satisfy the compatibility conditions.

Reconstruction of the displacement field

The above results are all about the strain field, while the tra- ditional methods (e.g., FEM) give the displacement field directly.The displacement field plays a key role in the mechanical analy- sis, especially when we consider the current configuration.The dis- placement field is not easy to be derived from the strain field ac- cording to the strain-displacement Eq.(7) .According to the Green’s method Eq.(13) , we are able to reconstruct the displacement field from the strain results, i.e.,

where the average of fluctuation partis determined by the cen- troid drift or the displacements of given points.To ensure the ac- curacy, we usually reconstruct the displacement field from the re- fined strainεrefine and stressσrefine.The convolution can also be computed in the Fourier space, so Eq.(28) becomes

The reconstructed displacement field is obviously continuous.Note that if the strain and stress are exact, the displacement is also exact.

Convergence analysis

Now, we study the convergence of iteration scheme for strain refinement.SCA provides a selection of the reference material moduli, i.e., the homogeneous material moduli.Due to different solution schemes, this selection can not always ensure the con- vergence of the iteration scheme, so we might use another refer- ence material in the refinement.If the material is homogeneous, the algorithm is a linear problem.We study this simple case to give some instructions for the selection of the reference modulusC0,refineintherefinement.

Let the stiffness of the material be constantC, and that of the reference material beC0.It is remarked that we do not need the inverse Fourier transform in the iteration process, so the iteration scheme reduces to

The initial value is considered as a non-zero fluctuation field.If we denote, we can rewrite Eq.(30) as

with

Due to ?Φ(0)= 0 , the last term on the right side vanishes.Note that the convergence of this iteration scheme depends on eigenval- ues ofA(ξ)at each wave numberξ.

We focus on a simple case in 2D, i.e.,ξ= [ξ1,0]T.We have the iteration scheme in the matrix form:

The eigenvalues ofA(ξ)isWhen

whereλ(x)andμ(x)are the material properties at each material point.

In summary, based on this, we suggest to take a harder stiffnessC0, such as the hardest phase as the reference material.

Numerical tests

Twodimensionalelasticmaterials

We apply our method to 2D plane strain problem to illustrate the efficiency of the reconstruction.

In Fig.1 , the color blue represents the matrix and the color yellow represents the inclusion.The volume fraction of the inclu- sion is 20%.The Young’s moduli and Poisson’s ratios of the matrix (phase 1) and inclusion (phase 2) are:

The mesh size for the high-fidelity RVE model is 600 ×600 .Here, the reference result (DNS) is obtained with the FFT-based method with the basic scheme [ 10 , 15 ].The results under the uni- axial tension loadingε11 = 0.001 are used as data-base of cluster- ing.In an offline stage, the mesh is compressed intok1= 8 clusters in the matrix andk2clusters in the inclusion by SOM.As the vol- ume fraction of the inclusion phase is 20%, we takek2 =k1/4 in the clustering process.

Under the uniaxial tension loadingε11 = 0.001 , the problem is a linear elastic one.The strain field obtained by SCA is plotted in Fig.2 b.The reference results are shown in Fig.2 a for comparison.It is noticed that the SCA result describes the strain field roughly and misses some features.Let us focus on the part in the solid line-box in the figures.In Fig.2 a, there exists a high strain region (in yellow) in the matrix.Figure 2 b undermines such significant structure in the SCA strain field.

In the refinement, the reference material takes Young’s modulus 4 ×105MPa and the Poisson’s ratioν0.2 (λ= 1.11 ×105MPa,μ= 1.67 ×105MPa), i.e., the inclusion phase.Figure 3 presents recon- struction results at different iteration steps.After merely one it- eration step, the high strain region appears.After three iteration steps, the strain field is captured very well, as compared with the FFT result.

We take the sectiony= 300 and observe the distribution curves ofε11 of different methods alongx-direction as shown in Fig.4 .The curve of the SCA result is piece-wise constant due to the clustering-based reduced order model.It correctly captures the FFT result in the average, but can not reflect the extreme values cor- rectly.After one iteration step, the distribution in the matrix is al- ready consistent with the FFT result, yet not so in the inclusion.The refinement after three iteration steps fits well the FFT result, including the extreme values.

In Table 1 , we list the strain energy

Table 1 Strain energy and deviation from FFT under the loading ε 11 = 0 . 001 .

and the deviations of the results from DNS

of different methods.The results converge to the reference solu- tion very fast.After 3 iteration steps, the error is less than 3%.We remark that strain at some material points might be negative after only one iteration step, but the result is good enough after three iteration steps.

In SCA, the homogeneous moduli are selected as the reference one, i.e.,λ0= 5.78 ×104MPa,μ0= 3.4 ×104MPa .This set of ref- erence material properties do not satisfy the condition Eq.(35) and the iteration scheme can not converge.

Hyper-elastic materials

The iteration scheme can be extended to the finite strain prob- lem.We consider a simple hyper-elastic model, i.e., a linear re- lation between the second Piola-KirchhoffstressSand the Green strainG

with

Cis the standard fourth-order isotropic stiffness tensor.Fis the deformation gradient.

We use the same geometry as shown in Fig.1 , and the same material properties as shown in Eq.(36) .Here, the reference result (DNS) is obtained with the FFT-based method [16] .We stretch the material alongx1 -direction to 20% deformation, well beyond the limit of the small strain problem.

We still select the inclusion phase as the reference material.The strain distributions of different methods along the sectiony= 300 are plotted in Fig.5 .The deviations from FFT are listed in Table 2 .After three iteration steps, the refinement is able to capture the strain field with less than a 3% difference from the FFT result.This verifies the efficiency of the reconstruction method.

Table 2 The deviation from FFT under the loading F 11 = 1 . 2

To summarize, in this paper, we propose a reconstruction method to accurately resolve the displacement, strain and stress, and recover the lost detailed features of reduced-order model.We first propose an iteration scheme based on the Lippmann- Schwinger equation with the coarse strain field of clustering- based homogenization and achieve high resolution, and further construct an accurate displacement field using Green’s method.Then we start from a simple elastic case to explore the con- vergence and the selection of the reference material.In the numerical tests, we refine the SCA results by the reconstruc- tion method, which shows the fast convergence and accuracy of the method.The set of construction iterations can be extended to other clustering-based homogenization with different kernel functions.

We are now working on the elasto-plastic cases to expand its applications.We consider to start simulations from the clustering- based strain fields based on the rough state variables provided by clustering-based methods.Compared to the traditional way to sim- ulate the whole loading process in an incremental way from zero fields, the reconstruction based on these rough information has the potential to achieve high efficiency.

Declaration of interests

The authors declare that they have no known competing finan- cial interests or personal relationships that could have appeared to influence the work reported in this paper.

Declaration of Competing Interest

The authors declare no conflict of interest.

Acknowledgements

Lei Zhang and Shaoqiang Tang were supported partially by the National Natural Science Foundation of China under grant numbers 11832001, 11521202, 11890681 and 11988102.

Appendix: Finite strain formulation

The finite strain problem is governed by the following equations in the undeformed configuration:

wherePis the first Piola-Kirchhoffstress,Fis the deformation gra- dient at each material pointX.We consider the periodic boundary condition.

For deriving the Lippmann-Schwinger equation, the reference homogeneous reference material is described by a constitutive re- lation:

C0isthestandardfourth-orderisotropicelasticstiffnesstensor.

Similar to Eq.(15) , the finite strain system can also be reformu- lated by the Lippmann-Schwinger equation

Φis the fourth-order tensor defined by

in the Fourier space.

After clustering the material points based on DNS results, we can also obtain the discrete Lippmann-Schwinger equation