Shoqing Tng , , Xi Zhu

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

Keywords:Clustering Finite element method Numerical homogenization

ABSTRACT Based on strain-clustering via k -means, we decompose computational domain into clusters of possibly disjoint cells. Assuming cells in each cluster take the same strain, we reconstruct displacement field.We further propose a new way to condensate the governing equations of displacement-based finite element method to reduce the complexity while maintain the accuracy. Numerical examples are presented to illustrate the efficiency of the clustering solver. Numerical convergence studies are performed for the examples. Avoiding complexities which is common in existing clustering analysis methods, the proposed clustering solver is easy to implement, particularly for numerical homogenization using commercial softwares.

Finite element method (FEM) represents a prevailing methodology for numerical explorations in many scientific and engineering disciplines [1] . While being superior in capturing details and complexities of geometrical structures and material properties, it often results in immense algebraic system, which greatly challenges computing and memory resources, and sometimes even becomes prohibitive. In a rather sharp contrast, only a much smaller amount/aspect of data are extracted and presented for a number of practical applications. One such example is numerical homogenization.

The goal of homogenization is to predict responses under loading and unloading for material of given structure. Only accounting for the composition ratio, some analytical approaches, such as Voigt/Reuss/Hashin-Shtrikman bounds and Mori-Tanaka theory,provide quick and explicit answer in the linear elastic regime [2,3] .To get a more accurate prediction and for a larger regime of interest, numerical homogenization is performed usually by FEM at the expense of fine resolution for constituents. To alleviate the heavy computing cost thereby caused, model reduction strategies have been proposed, such as transformation field analysis (TFA),non-uniform transformation field analysis (NTFA), proper generalized decomposition (PGD), etc. One decomposes the state variables into modes, and predicts the average behaviours based on such decompositions[ 4–6 ]. More recently, Liu et al [ 7 ]. proposed a novel data-driven model reduction strategy called self-consistent clustering analysis (SCA). Afterwards, other clustering analysis methods have been developed, including virtual clustering analysis (VCA)[8] , FEM-based clustering analysis (FCA) [9] , and Hashin-Shtrikman type finite element method (HSFE) [10] . For more studies along this line, please refer to Refs. [11,12] .

As pointed out in Ref. [8] , the key assumption “once respond similarly, always similarly” enables a domain decomposition and therefore the model reduction. Behind this is the observation that strain and stress, not the displacement, are intrinsic quantities for the material. In an off-line stage, detailed numerical simulations for one or several selected loadings form a database of strain, with which clustering is performed via k-means or self-organizing map(SOM). This domain decomposition, unlike traditional ones, yields a prescribed number of clusters, each consisting disjoint finite element cells. In another word, we take a Lebesgue point of view instead of a Riemannian one.

The second key idea is a weak form of governing system, using a reference material of homogeneous linear elasticity. The true material response is substituted by a body force in terms of stress polarization. As the linear reference problem with body force is readily solved with Green’s functions, this leads to an integral Lippmann-Schwinger equation. Under the above assumption of constant strain in each cluster, numerical discretization gives a system of discrete Lippmann-Schwinger equation in much reduced degrees of freedom (only a multiple of cluster numbers). Substantially enhanced efficiency is then gained in the online prediction stage.

The reference material property determines an interaction tensor (discrete Green’s operator), representing the influence across different clusters. Accurate evaluation of the interaction tensor is expensive, as it involves operations in the order of a power of the finite element cell number. In such a computation with already substantial reduction compared with the original FEM, we found that the evaluation dominates the CPU time in off-line stage, and is much longer than the online prediction. In view of this, fast evaluation was proposed to alleviate the heavy load through coarsening,which further improves considerably the efficiency [13] . Moreover,the choice of reference material property affects numerical homogenization results through boundary terms, though rather weakly in the elastic regime [14] . One either uses a self-consistent iteration or adaptively adjusts Young’s modulus, both with additional cost.

In this letter, we use strain clustering strategy for domain decomposition, and condensate the displacement formulation of finite element analysis. In this way, we avoid the empirical or iterative choice of reference material, the derivation of Lippmann-Schwinger equation, and the expensive evaluation of interaction tensor. Last but not least, as displacement formulation is most commonly used in commercial softwares, this allows a fast and easy implementation for our clustering solver.

The whole idea is to express the displacement through constant strain over each cluster. Then the stiffness matrix is transformed to that for clusterwise strain. In the following, we shall first describe such a construction in one and two space dimension(s). Then we form the governing algebraic system for the clusterwise strain. Finally, we describe several numerical examples and make some discussions.

By definition, strain is the symmetrized gradient of displacement.

Under this assumption, we reconstruct the displacement field. It is emphasized thatεI(I= 1,2,···,k)are taken as the main variables in the algebraic system formulated below.

In one space dimension, strain is a scalar and we denote it asε. In thei-th interval, by integration we may readily get

In a rectangular domain [0,Lx] ×[0,Ly] , displacement may be obtained through

Here the boundary termsu(0,y),v(x,0)are chosen according to the loading/boundary conditions. If appropriately chosen,εxyevaluated from the displacement field automatically satisfies the consistency condition Eq. (6) , and should comply with the given strain field.

Numerically, we construct piecewise linear displacement field in(i1,i2)-th cell, centered at(xi1-1/2,yi2-1/2)=((i1 -1 / 2)Δx,(i2 -1 / 2)Δy), as follows.

First along the liney=yi2, we take the mean strain on its two sides, then integrate inxdirection to get forx∈ [xi1-1,xi] ,

The partial derivativevx(xi1-1/2,yi2-1/2)may be computed accordingly, which together with Eq. (11) forms the shear strainεxy.We take the cluster average ofεxyto define the clusterwise shear strain, and reassign it to all cells of this cluster. This step does not make any difference to the effective modulus in the linear elastic regime. But when nonlinearity appears, this averaging and reassignment step is suggested.

Notice thatu(0,yi2),v(xi1,0)may be adjusted to comply with the loading conditions. However, we choose them to be 0 in the following numerical examples for simplicity.

The expressions are similar for three space dimensions.

Fig. 1. Numerical example in one space dimension. a Strain and 8 clusters; b Displacement; c Stress; d Relative error with k = 4 , 8 , 16 , 32 , 64 .

In displacement-based FEM, the governing algebraic system is composed of equations for boundary conditions, and those for energy minimization. The latter is the major part, which may lead to a large system in general, hence becomes the target for model reduction.

For the sake of clarity, here we formulate explicitly the reduced governing system only in one space dimension. Extensions to multiple dimensions are straightforward. In this case, a typical displacement formulation with minimization reads

As usual, we denote U =(u1,u2,···,uN)Tcollecting all nodal displacements (except the assigned 0 displacement atx0 = 0 ). The discrete form for minimization yields

Here D corresponds to the gradient operator, leading to a rank deficiency in K . The linear system therefore is supplemented by boundary conditions to be uniquely solvable. For instance, it can be

We formally denote such condition(s) as BU = ～ b, and add it to Eq. (14) as a penalty to get a full-ranked system

whereαis a large number.

The reconstruction relates displacement with the clusterwise strain X =(ε1,ε2,···,εk)T. The linear reconstruction may be expressed in terms of a matrix P such that U = PX . Substituting this into the minimization problem we get a reduced system

This comes from the minimization of a discrete energy

We illustrate the proposed method by two numerical examples,in one and two space dimension(s), respectively.

Fig. 2. Material structure in two space dimensions: Each unitary square grid takes a dimension 15 ×15 , no-mark grids stand for matrix, marked ones stand for two types of inclusions.

Example 1

Applying a strain ˉε= 0.05 , i.e.,u(-1)= 0,u(3)= 0.2 , the analytical solution gives ˉσ= 0.015 . The numerically solved strain with 8 clusters is displayed in Fig. 1 a. We may reconstruct the displacement out of this clusterwise strain, which is piecewise linear and plotted in Fig. 1 b. With only 8 clusters (or 8 slopes), it very well reproduces the reference solution. The stress distribution corresponding to the piecewise constant strain is displayed in Fig. 1 c.We see fluctuations around the exact solutionσ= 0.015 . The average stress is 0.0151, amounting to a relative error of 0.64% . This is also the relative error for the effective elastic modulus. When we increase the number of clusters, prediction for the effective modulus becomes more accurate. In Fig. 1 d, we see a convergence at a rate 2.2955. Note that in VCA, the rate is 1.96.

Example 2

We consider a material in the domain [0,75] ×[0,75] . See Fig. 2 . The matrix takes Young’s modulusE1= 100 MPa and Poisson’s ratioν1= 0.3 . The two inclusions takeE2 = 500 MPa,ν2=0.19 andE3= 300 MPa,ν3= 0.25 in grids marked by 2 and 3, respectively.

The boundary condition corresponding to a uniaxial extension is prescribed as

TakingΔx=Δy= 1 , we perform FEM computation to get the reference solution. Then we cluster the 75 ×75 cells into three categories. We putk1= 40,k2= 16,k3= 4 clusters in the matrix and inclusion sub-domains, respectively. See Fig. 3 for the reference solution and clustering.

Fig. 3. Example 2 computed by FEM: a Strain component ε 11 ; b Clustering with k 1 = 40 , k 2 = 16 , k 3 = 4 ; c Displacement u ; d Displacement v .

Fig. 4. Example 2 computed by clustering solver ( k 1 = 40 , k 2 = 16 , k 3 = 4 ): a Strain component ε 11 ; b Stress component σ11 ; c Displacement u ; d Displacement v .

The corresponding numerical results by the clustering solver is displayed in Fig. 4 . We observe a rather good agreement with the reference solution on the strain componentε11, which is the dominant one under such loading. Notice that the degrees of freedom is now less than 2% of the FEM algebraic system. The reconstructed displacementualso well reproduces the reference solution. The vertical one, which is in a much smaller magnitude,is not very well captured. This induces a certain amount of error in the stress field. Nevertheless, after averaging, it does not cause much problem in the numerical homogenization result. As a matter of fact, with the reference average stressσ11= 8.421546 MPa,σ22= 3.290498 MPa,σ12= 0.030550 MPa , the clustering solver givesσ11= 8.551933 MPa,σ22= 3.236085 MPa,σ12= 0.0 6 6151 MPa .Except for the shear stress, which is again rather small, the other two components take error below 1.7%.

The numerical error may be further reduced if we increase the number of clusters. As a limiting case, if we assign a cluster for each cell, i.e., 75 ×75 clusters altogether,the numerical results amount toσ11= 8.375391 MPa,σ22=3.29030 6 MPa,σ12= 0.030 620 MPa , corresponding to relative error 0.55%,0.006%,0.23% , respectively. Such error stems from our numerical reconstruction and reformulation of the governing equations. It manifests the convergence of the proposed clustering solver.

We remark that our preliminary numerical study suggests that enough of clusters should be assigned to the matrix to maintain the accuracy.

In this letter, we propose an easy-to-implement approach for clustering analysis. Same as SCA and VCA, we use one set of well resolved strain data to cluster the FEM cells. By a simple reconstruction, the displacements may be expressed in terms of clusterwise strain. The governing system on such clusterwise strain reduces substantially the degrees of freedom, which is formulated through minimization. Preliminary numerical examples in one and two space dimension(s) illustrate the efficiency of the proposed clustering solver. In particular, convergence studies in both examples manifest the accuracy. The most distinct feature of this approach is its readiness to be implemented within commercial softwares. Without the complexities of forming interaction tensors and choosing suitable reference materials, the numerical cost is also much reduced.

The potential application of this method are two-folded. First,it may be adopted for numerical homogenization, where the effective/average stress and strain are under consideration. Secondly, it may also serve as an upscaling strategy for solving the displacement (and therefore stress and strain) fields. The raw solution by reconstruction provides a good guess for iterative schemes to solve the original problem.

This letter only illustrates the potential of the method. In fact,there are many aspects for further extensions and improvements.For one thing, the boundary condition for displacement reconstruction leaves a big room for better complying with true mechanical loading. That will enhance the stress field resolution as well.Meanwhile, with the provided numerical examples and previous experiences on clustering analysis, it is promising to predict effectively homogenization results for elasticity, hyper-elasticity, plasticity and large deformation problems[ 15 ]. These are our on-going research, and we shall report the results in future publications.

Declaration of Competing Interest

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

Acknowledgments

This work was supported in part by the National Natural Science Foundation of China (Grant Nos. 11832001, 11890681 and 11988102).