{'fraction_non_alphanumeric': 0.03287304204451773, 'fraction_numerical': 0.013345012366034624, 'mean_word_length': 4.695129107981221, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
LIMIT ANALYSIS OF MICROSTRUCTURES BASED ON HOMOGENIZATION THEORY AND THE ELEMENT-FREE GALERKIN METHOD

This paper presents a novel numerical formulation of computational homogenization analysis of materials at limit state. The fluctuating displacement field are approximated using the Element-Free Galerkin (EFG) meshless method. The estimated yield surface of materials can be determined by handling the multiscale (macro-micro) transition. Taking advantage of high-order EFG shape function and the second-order cone programming, the resulting optimization problem can be solved rapidly with the great accuracy. Several benchmark examples will be investigated to demonstrate the computational efficiency of proposed method.


INTRODUCTION
As the increasing use of composite and heterogeneous materials in practical engineering structures, the estimation of their effective properties plays a vital role in safety assessment as well as structural design. The elastic-plastic incremental method can be employed to predict the ultimate load and collapse mechanism of structures. However, direct method, e.g. limit analysis shows more effectively, i.e. the critical status of structures can be determined without any knowledge of whole loading path history [1,2].
Dealing with the microstructures, the multi-scale homogenization technique or socalled global-local analysis firstly introduced in [3] has been widely exploited in recent years. The computational homogenization methodology has been mostly applied to the periodic composite and heterogeneous materials. Based on the concept of representative volume element (RVE) and homogenization technique in [3], the first formulation of limit analysis in terms of solving the composite structure at micro-scale has been proposed in [4]. Then, the theoretical formulation was developed by [5,6] for fiber-reinforced composite using Drucker-Prager, Mohr-Coulomb or von Mises yield criterion. The first numerical implementation for this field has been reported in [7] with the use of finite element method and linear mathematical programming. By means of static direct methods, a three-dimensions finite element procedure for analysis of isotropic microstructures was developed by [8,9]. Besides, using the similar approach, a quasi-lower bound formulation for periodic composite and heterogeneous materials using the nonlinear programming was presented in [10]. In contrast, the kinematic formulations in combination with nonlinear algorithms can be found in [11][12][13][14][15]. In those works, both of isotropic and anisotropic materials obeying the von Mises or elliptic yield criterion were considered. In order to improve the computational aspect, a numerical procedure based on combination of kinematic limit analysis and homogenization theory for periodic materials proposed in [16]. However, these studies only focused on the application of finite element method.
This study aims to develop a novel computational homogenization approach for upper bound limit analysis of microstructures using Element-Free Galerkin method known as one of most successful procedures for yield design problems. The stability conforming nodal integration (SCNI) technique proposed in [17] is employed to improve the performance of the numerical formulation. In addition, the plastic dissipation will be transformed into the form of a sum of norms and the resulting optimization are then formulated as conic one. The benchmark numerical examples will be considered and the good agreement in comparison to previous procedures proves the performance of present method.

Homogenization theory
Consider a heterogeneous representative volume element (RVE) Ω ∈ R 2 at every material point x ∈ V, where V ∈ R 2 denotes the heterogeneous macroscopic-continuum. The micro-structure is subjected to the body force f, the surface load t on the static boundary Γ t and fixed by the displacement field u on the kinematic boundary Γ u .
The micro-scale problem can be treated as the boundary value one in solid mechanics, where the overall strain E are transferred to micro-structure in form of kinematic boundary constrains. At microscopic scale, the local fields is decomposed into two parts including mean and fluctuation terms. The displacement, strain and stress fields of microstructure are now given by where matrix X consists of the position components of each material point in the computational domain; Σ is the overall stress;ũ(x),˜ (x) andσ(x) denote the fluctuation parts of displacement, strain and stress rate.
For the enforcement of boundary conditions, the most efficient in terms of convergence rate so-called periodic procedure, where there are the periodicity of fluctuation displacement field and anti-periodicity of traction field on RVE boundary, is employed in this studyũ whereũ x + andũ x − are the fluctuation displacement field, t x + and t x − are the traction field of positive and negative boundaries, respectively. Note that regarding to the periodic characteristic of the fluctuation terms, the average of˜ (x) andσ(x) over the RVE should vanish, it means where the operation · stands the volume average of fields over the RVE. Denoting |Ω| for the area of RVE, the macroscopic quantities can be calculated from the microscopic ones via the average relations For any admissible velocity and stress field satisfying the periodic and anti-periodic conditions on boundary, the principle of macroscopic virtual work can be expressed as

Fundamental of limit analysis
Denoting X and Y for the appropriate spaces of statically admissible stress state and of kinematically admissible velocity state, respectively. The exact collapse multiplier will be obtained if one of following optimization problems is solved ; the plastic dissipation rate D(u) and the external work F(u) can be expressed in terms of σ and u as follow It is worth noting that most of yield criterion can be expressed in the following form where P is the coefficient matrix consisting of strength properties of materials. For isotropic materials, the von Mises criterion is frequently applied, and matrix P for plane stress problem can be expressed as where σ p is the isotropic yield stress of material.
In framework of limit analysis, the strain rates are assumed to obey the normality rule. Denoting Θ for the inversion matrix of P, the power of dissipation can be formulated in terms of strain rates as

Kinematic limit analysis based on homogenization theory
Assuming that all constitutions of ductile composite are rigid-perfectly plastic and the strain of constitutions obey the normality rule. The kinematic approach in framework of limit analysis for computation homogenization described in [11,12,14,15,18,19] will be taken into account in this paper. According to the homogenization theory, all variables related to the microscopic structures are split into two parts including mean fields averaged over RVE and fluctuation fields.
Omitting the body force f and applying the principle of microscopic virtual work (5), the normalization condition of external power can be rewritten as Now, the kinematic limit formulation of computational homogenization analysis for a periodic micro-structure can be expressed as Solving the nonlinear problem (12), the upper-bound of macroscopic limit strength λ + Σ will be determined. It should be noted that present study only considers the continuous velocity fields, and in cases when the velocity fields are assumed to be discontinuous, the dissipated power generated by discontinuities must be taken into account.

ELEMENT-FREE GALERKIN METHOD
The approximate function at point x reflected via a set of scattered nodes in arbitrary domain Ω using moving least square (MLS) technique can be expressed as where the coefficient vector a(x) and the basis function vector p(x) are defined by The complete polynomial basis is given by The unknown coefficient is minimized, where u i denotes the nodal value of approximate function at point i th , N is number of points in the neighborhood of x, w(x − x i ) is the weighting function depending on the influent domain of point x. This study employs the isotropic quadratic spline Minimizing J(x) in Eq. (16) leads to the system linear equation with Taking an inversion, a(x) can be obtained by solving Eq. (18) as Now, the approximate function can be expressed as where the EFG shape function is defined by In order to improve the computational efficiency of the EFG method, the stabilized conforming nodal integration presented in [17] is applied here to obtain stabilized shape function derivatives asΦ where A C is the area of a Voronoi smoothing domain with boundary Γ C , n α is the outward surface normal vector. Detailed calculation of the smoothing shape function derivatives Φ i,α can be found in [19].

DISCRETE FORMULATION USING EFG METHOD
In the kinematic formulation, the local microscopic fluctuation displacementũ(x) and fluctuation strain˜ (x) at point x can be approximated using EFG method as where N(x) denotes the iRBF shape function; B(x) is the strain-displacement matrix consisting the smoothed version of shape function derivatives; and d is the vector containing the nodal fluctuation displacement components.
With the use of SCNI technique for the numerical integration, the plastic dissipation well-known as the objective function of the optimization problems can be expressed as where σ p is the yield stress of material, A i is the area of the i th nodal representative domain created using Voronoi diagrams. In this study, the optimization problem will be formulated in the form of second order cone programming (SOCP) ensuring that it can be solved using the highly efficient solves. Hence, a form of sum of norm can be used to calculate the internal dissipation power as where · denotes the Euclidean norms and σ is the vector of additional variables defined by with Q denotes the Cholesky factor of Θ. Next, the periodic feature of the fluctuation displacement for nodes on the boundary of RVE need to be enforced. Denoting Γ + and Γ − for the positive and negative boundary such that Γ + ∪ Γ − = Γ and Γ + ∩ Γ − = ∅, the periodic boundary condition for each pair x + , x − of boundary material points can be expressed as Assembling to the global matrix, Eq. (28) can be rewritten as Finally, by introducing the auxiliary variables (t 1 , t 2 , . . . , t N ), the optimization problem can be formulated in the form of conic programming as follows

NUMERICAL EXAMPLES
In this section, the computational aspect of proposed method will be investigated. An unit RVE of a × a = 1 × 1 mm is used for all examples. The problems are considered in plane stress, and thus number of variables N var is equal to 6 × N + 3. The resultant optimization problems are solved using the commercial software package Mosek on a 2.8 GHz Intel Core i7 PC running Window 10. such that + [ = and + \ = ?, the periodic boundary condition for each pair {x + , x } of boundary material points can be expressed as Assembling to the global matrix, equation (28) can be rewritten as Finally, by introducing the auxiliary variables (t 1 , t 2 , . . . , t N ), the optimization problem can be formulated in the form of conic programming as follows s.t

NUMERICAL EXAMPLES
In this section, the computational aspect of proposed method will be investigated. An unit RVE of a ⇥ a = 1 ⇥ 1 mm is used for all examples. The problems are considered in plane stress, and thus number of variables N var is equal to 6 ⇥ N + 3. The resultant optimization problems are solved using the commercial software package Mosek on a 2.8 GHz Intel Core i7 PC running Window 10. The e↵ect of microscopic hole on the overall strength of perforated materials is investigated. Two perforated material models including a rectangular hole of size L 1 ⇥ L 2 and a circle hole of radius R at center are considered. The problems have been investigated using kinematic formulation in [11,12,16] and quasi-static formulation in [10]. The RVE is subjected to an orthogonal macroscopic stress (⌃ 11 , ⌃ 22 ) in plane (x 1 , x 2 ) as shown in (a) Rectangular hole RVE such that + [ = and + \ = ?, the periodic boundary condition for each pair {x + , x } of boundary material points can be expressed as Assembling to the global matrix, equation (28) can be rewritten as Finally, by introducing the auxiliary variables (t 1 , t 2 , . . . , t N ), the optimization problem can be formulated in the form of conic programming as follows s.t

NUMERICAL EXAMPLES
In this section, the computational aspect of proposed method will be investigated. An unit RVE of a ⇥ a = 1 ⇥ 1 mm is used for all examples. The problems are considered in plane stress, and thus number of variables N var is equal to 6 ⇥ N + 3. The resultant optimization problems are solved using the commercial software package Mosek on a 2.8 GHz Intel Core i7 PC running Window 10. The e↵ect of microscopic hole on the overall strength of perforated materials is investigated. Two perforated material models including a rectangular hole of size L 1 ⇥ L 2 and a circle hole of radius R at center are considered. The problems have been investigated using kinematic formulation in [11,12,16] and quasi-static formulation in [10]. The RVE is subjected to an orthogonal macroscopic stress (⌃ 11 , ⌃ 22 ) in plane (x 1 , x 2 ) as shown in The effect of microscopic hole on the overall strength of perforated materials is investigated. Two perforated material models including a rectangular hole of size L 1 × L 2 and a circle hole of radius R at center are considered. The problems have been investigated using kinematic formulation in [11,12,16] and quasi-static formulation in [10]. The RVE is subjected to an orthogonal macroscopic stress (Σ 11 , Σ 22 ) in plane (x 1 , x 2 ) as shown in Fig. 1, where θ denotes the angle between the principle stress and x 1 -axis. The matrix material of rectangular hole RVE is aluminum Al whose yield stress is σ p = 137 MPa, while those of RVE with circular hole is mild steel St3s with yield stress σ p = 273 MPa. Material behavior is assumed to obey the von Mises yield criterion. Fig. 2 illustrates the scheme of nodal discretization using Voronoi diagram.    sdof denotes the total system degrees of freedom In case of rectangular hole RVE, the problem is considered with different dimensions of hole involving (L 1 × L 2 = 0.1 × 0.5 mm) and (L 1 × L 2 = 0.1 × 0.7 mm). Tab. 1 shows present numerical solutions compared with those in [12,16]. It can be observed that proposed procedure can prove the highly accurate solutions with low computational cost in comparison with approaches presented in [12,16]. The number of variables in present procedure is less than those used in [16], while the numerical result is approximate. Moreover, with the use of cone-based algorithm the resultant optimization problem in present study can be solved rapidly, the CPU-time required for handling the problem is much less than those in [12] using iterative algorithm. Fig. 3 plots the macroscopic uniaxial strength Σ 11 for various values of angle θ corresponding to two different sizes of rectangular hole. It can be seen from both of sub-figures that the upper-bound solutions obtained using present iRBF are slightly lower (better) than available those in other studies. The figure also shows the good agreement of present solutions in comparison with previous numerical results reported in [12,16] as well as the experimental results reported by Litewka et al. (ref [12]).
For the case of perforated material with circular hole, various values of ratio R/a and loading angle θ are considered. Fig. 4 also plots the solutions reported in [11,16].     Next, the effect of microscopic hole size on the overall strength of perforated materials is investigated. Obviously, it is reasonable that the effective macroscopic strength decreases when increasing size of holes. For rectangular hole RVE, it can be observed that for all values of angle θ, the obtained macroscopic strengths in case of L 1 × L 2 = 0.1 × 0.7 mm are lower than those in case of L 1 × L 2 = 0.1 × 0.5 mm, as seen in Fig. 3. For the case of RVE with circular hole, the solutions in Fig. 4 shows the similar trend when increasing ratio R/a. Material behavior is assumed to obey the von Mises yield criterion. Figure 2 illustrates the scheme of nodal discretization using Voronoi diagram. In case of rectangular hole RVE, the problem is considered with di↵erent dimensions of hole involving (L 1 ⇥ L 2 = 0.1 ⇥ 0.5 mm) and (L 1 ⇥ L 2 = 0.1 ⇥ 0.7 mm). Table 1 shows present numerical solutions compared with those in [12,16]. It can be observed that proposed procedure can prove the highly accurate solutions with low computational cost in comparison with approaches presented in [12,16]. The number of variables in present procedure is less than those used in [16], while the numerical result is approximate. Moreover, with the use of cone-based algorithm the resultant optimization problem in present study can be solved rapidly, the CPU-time required for handling the problem is much less than those in [12] using iterative algorithm. Figure 3 plots the macroscopic uniaxial strength ⌃ 11 for various values of angle ✓ corresponding to two di↵erent sizes of rectangular hole. It can be seen from both of sub-figures that the upper-bound solutions obtained using present iRBF are slightly lower (better) than available those in other studies. The figure Material behavior is assumed to obey the von Mises yield criterion. Figure 2 illustrates the scheme of nodal discretization using Voronoi diagram. In case of rectangular hole RVE, the problem is considered with di↵erent dimensions of hole involving (L 1 ⇥ L 2 = 0.1 ⇥ 0.5 mm) and (L 1 ⇥ L 2 = 0.1 ⇥ 0.7 mm). Table 1 shows present numerical solutions compared with those in [12,16]. It can be observed that proposed procedure can prove the highly accurate solutions with low computational cost in comparison with approaches presented in [12,16]. The number of variables in present procedure is less than those used in [16], while the numerical result is approximate. Moreover, with the use of cone-based algorithm the resultant optimization problem in present study can be solved rapidly, the CPU-time required for handling the problem is much less than those in [12] using iterative algorithm. Figure 3 plots the macroscopic uniaxial strength ⌃ 11 for various values of angle ✓ corresponding to two di↵erent sizes of rectangular hole. It can be seen from both of sub-figures that the upper-bound solutions obtained using present iRBF are slightly lower (better) than available those in other studies. The figure

CONCLUSIONS
The plastic limit strength and the collapse mechanism of materials has been studied using the combination of direct analysis and homogenization theory. By means of second-order cone programming and the EFG approximation, the resulting optimization problems are kept in minimum size and solved rapidly. The good agreement of numerical solutions in comparison with other studies shows the computational efficiency of proposed method. In future work, the plane strain or three-dimensions problems, for which the volumetric locking phenomena is required to be handled, are extended. In addition, more complicate effects, e.g. material interfaces, multiple crack or even variable loading should be considered.