A Study of the 1 + 2 Partitioning Scheme of Fibrous Unitcell under Reduced-Order Homogenization Method with Analytical Influence Functions

1  Introduction

Fibrous composite materials with the advantages of high specific strength, specific stiffness, and corrosion resistance are widely used in many advanced manufacturing sectors, and that leads to the demand for high-fidelity evaluation of mechanical properties of the composite structures [1–3]. However, the overall mechanical response of composite materials is complicated because of the large gap between the mechanical properties of each constituent, especially when the nonlinearities are considered. The various mesoscale geometric structures of fibrous composite materials and loading patterns also lead to different types of material failure, such as breakage of fibers, debonding of fiber/matrix interface, or delamination of plies. Therefore, developing an efficient framework for the constitutive model of composite materials with comprehensive consideration of various nonlinear material behaviors is important, and a series of research has been conducted in this area. Firstly, various homogenized anisotropic damage models of the composite material were proposed [4–6] and applied in engineering. Meanwhile, based on the high-resolution modeling, many investigations were conducted to capture the failure mechanism on the mesoscale [7–9]. In these homogenized damage models, the loss of microscale information usually leads to a loss of accuracy and generalizability; the increasing complexity of the homogenized model also requires extensive calibration. On the other hand, the orders of magnitude gap between the characteristic length of the mesoscale model and composite structure gives rise to unacceptable computation costs. Therefore, the so-called multiscale method based on continuum mechanics, also labeled as an upscaling method, has been developed rapidly to balance fidelity and computation efficiency. A series of works on direct multiscale methods, such as the multilevel finite element (FE2) approach [10], the first-order nonlinear computational homogenization methods [11], and the higher-order homogenization [12], have made significant improvements in addressing this dilemma of scales. However, the number of unknown variables within a single unitcell is still too large for engineering applications. This required further model simplification or model reduction at the mesoscale, leading to the proposal of reduced-order models for heterogeneous continua.

Transformation field analysis (TFA) is a typical reduced model proposed by Dvorak for the analysis of elastic composite [13] and then extended for the inelastic cases [14]. The scale transitions of strain and stress in heterogeneous structures were expressed explicitly by this approach, where plastic strain and thermal expansion were treated as eigenstrain. In TFA, the transformation influence functions and concentration factors are precomputed based on the mechanical parameters of different phases and the geometric structure of the composite. The number of variables was reduced significantly due to the assumption of piecewise uniform fields. Therefore, this approach was utilized and developed for the analyses of composite structures by many scholars over the past three decades [15–17]. However, some deficiencies of classical TFA were recognized in practice. First, the prediction of macroscopic stress is usually too stiff due to the assumption of piecewise uniform plastic strain [18,19]. In addition, expectedly, the accuracy of the two-phase system in inelastic analysis is unsatisfactory when following the principle of so-called “1-phase-1-partition”. That leads to the finer division of each individual phase, in other words, a larger number of internal variables in overall constitutive relations [20]. Instead of finer division in the space domain of each phase, decomposing the plastic strain into a finite combination set of plastic modes, the so-called nonuniform transformation field analysis (NTFA) [21,22], is another option for reaching higher precision. However, the extra work for calibrating empirical laws and extensive exploration of the deformation space under plastic deformation can counteract the advantage of NTFA on computational efficiency.

Liu et al. proposed the self-consistent clustering analysis (SCA) method [23] to overcome the limitations of TFA and its derivatives. In the SCA method, the constituents of composites are divided into a specific number of subdomains by k-means clustering based on the linear response of the unitcell, while a self-consistent implicit scheme is adopted for updating the macroscopic tangent modulus. Due to the high computation efficiency on heterogeneous materials, SCA is absorbing great attention and has been developed for a series of inelastic response analyses, including elastoplasticity, strain-softening, fatigue-strength prediction, and so on [24–27]. The extensions of SCA were also proposed for local mechanical analysis of composites and simplification of implementation [28,29]. In addition, k-means clustering is also combined with other reduced-order methods, such as NTFA [30,31]. On the other hand, other cluster-discretization methods, including the virtual clustering analysis (VCA) [32] and the FEM clustering analysis (FCA) [33], were developed as optional solutions of interaction tensors. Naturally, these cluster-discretization methods have been utilized for the research of machine learning-based multiscale modeling, especially in the generation of databases. SCA was applied for data mining and uncertainty analysis in a unified data-driven computational framework to design and model materials and structures [34]. A reduced-order model involving three scales was developed for woven composites using a data-driven SCA × SCA multiscale simulation framework [35]. Li et al. built a framework in which cluster-discretization methods are applied to develop a material behavior database for training neural networks, including feed-forward neural networks (FFNN) and convolutional neural networks (CNN) [36]. They presented its application on multiscale topology optimization.

Other machine learning-based multiscale techniques have also been developed as an important direction of model reduction. Zhang et al. [37] developed a reduced-order model of the hierarchical deep-learning neural networks (HiDeNN) based on a tensor decomposition (TD). A reduced-order machine learning finite element method based on proper generalized decomposition (PGD) reduced HiDeNN presents a good performance in the multiscale analysis of composite materials, topology optimization, and additive manufacturing [38]. Fish et al. [39,40] proposed a data-physics-driven reduced-order homogenization (dpROH) approach, which significantly improves the accuracy of physics-based reduced-order homogenization (pROH) by combining data generated from high-fidelity models with Bayesian inference. The artificial neural network (ANN) and deep material network (DMN) were also applied to homogenize the nonlinear mechanical properties of heterogeneous materials [41,42]. Masrouri et al. [43] utilized generative AI to predict the stress distribution of bicontinuous composites with complex irregular distribution.

The motivation of the present work is to develop a reduced-order homogenization (ROH) multiscale method [44–47] for fibrous composite material, in which the overestimation of strength due to “1-phase-1-partition” is improved while no extra internal variables of the constitutive model of each phase are introduced; in addition, simple and explicit partitioning of the subdomains of constituents is preferred over to the complicated and implicit partitioning of cluster-discretization methods. Therefore, a framework of constitutive model based on the ROH method is proposed, in which an adaptive method for partitioning the matrix phase is involved. This model can capture the stress concentration on the mesoscale and predict the overall response of fibrous composite with better precision since the matrix phase is divided into inner and outer matrix partitions with an independent nonlinear evolution process. The volume fraction of the inner matrix partition is predicted by a numerical calculation method, whose object function targets the minimization of error in strain energy density compared to the results from a modified Eshelby method [48–51]. In addition, the influence functions of elastic strain and eigenstrain are calculated analytically so that the offline computation is simplified drastically. In this manuscript, the framework of reduced-order homogenization is briefly introduced in Section 2; the mathematical manipulation of the analytical solution of influence functions for 1 + 2 partitions is presented in Section 3, along with the numerical calculation method for partitioning of the matrix phase; the numerical results for the verification of new model under different volume fractions of fiber are exhibited in Section 4; finally, the discussion of the numerical results and the corresponding conclusions are provided in Sections 5 and 6, respectively.

2  Reduced-Order Homogenization

The multiscale computational homogenization method based on the unitcell discretized by finite elements is usually computationally prohibitive for general macroscopic structures. It requires averaging and homogenization procedures over a converged unitcell problem at every macroscopic integration point during each Newton-Raphson iteration. In addition, the unitcell problem involving material nonlinearities needs to be solved using a Newton-Raphson iterative method. Model reduction technique is necessary to reduce the degrees of freedom by applying a specific distribution of mechanical response over phases at the mesoscale. The reduced-order homogenization (ROH) approach is a typical multiscale method with low computational cost, adopting a piecewise constant approximation of mechanical quantities over each partition.

In addition, the concept of eigenstrain is introduced to consider the material nonlinear response that

μ(α)=ε(α)−M(α):σ(α),α=1,2,⋯,M(1)

where M is the number of partitions, M(α)=[Mijkl(α)] is the compliance coefficient tensor for the α-th partition. The reduced-order homogenization method assumes that the material response, including the strain, stress, and state variables, keeps the same within one partition all the time and thus introduces the following reduced equilibrium equations (see Chapter 4 of the book [52]):

ε(β)−∑α=1MP(βα):μ(α)=E(β):εc,β=1,2,⋯,M(2)

where E(α)=[Eijkl(α)],α=1,2,⋯,M are the elastic strain influence functions for the α-th partition, and P(βα)=[Pijkl(βα)],α,β=1,2,⋯,M are the eigenstrain influence functions measuring the interaction between the α-th and β-th partitions, and εc is the macroscopic strain.

The elastic strain influence functions satisfy the following two constraints:

∑α=1McαE(α)=I(3)

∑α=1McαL(α):E(α)=Lc(4)

E(1)=1c1(L(1)−L(2))−1:(Lc−L(2))(5)

E(2)=1c2(L(2)−L(1))−1:(Lc−L(1))(6)

A two-partition unitcell model can be used for a general inclusion-matrix heterogeneous material under the assumption that each phase material always keeps the same response. Thus, it significantly reduces computational costs compared to the so-called direct-homogenization.

3  Fibrous Unitcell with 1 + 2 Partitions

The “1-phase-1-partition” reduced-order homogenization (ROH) method significantly decreases computational cost. However, under general multi-axial deformation, the distribution of material status within one phase is not uniform, which can give a stiffer material response than the direct-homogenization method due to the over-constraint of the material points. This study takes a fibrous unitcell as an example, as depicted in Fig. 1, and the matrix phase is partitioned into the inner and outer partitions. It visualizes the stress distribution under different loading modes, including uniaxial tension in longitudinal (Fig. 2a,b) and transverse (Fig. 2e,f) directions, shear in the longitudinal-transverse plane (Fig. 2c,d) and transverse plane (Fig. 2g,h). These numerical results show different stress magnitudes between the inner and outer partitions at the matrix phase.

Figure 1: A sketch of a three-partition fibrous unitcell with finite element discretization: (a) fibrous unitcell; (b) fiber phase; (c) inner matrix partition; and (d) outer matrix partition

Figure 2: A sketch of von Mises stress of 1st and 2nd matrix partition under different loading modes: (a and b) uniaxial tension along the longitudinal direction; (c and d) simple shear at the longitudinal-transverse plane; (e and f) uniaxial tension along the transverse direction; and (g and h) simple shear at the transverse plane

This study uses the ROH method for the fibrous unitcell with two partitions at the matrix phase. We define this partition as a “1 + 2” partition. Due to the transverse isotropy symmetry of fibrous composite material, this study divided the matrix phase into the inner and outer partitions, as depicted in Fig. 1c,d, respectively. In the rest of the paper, the properties and status associated with the fiber phase, matrix phase, inner matrix partition, outer matrix partition, unitcell (mesoscale), and coarse (macroscopic) scale are denoted as ◻f, ◻m, ◻a, ◻b, ◻uc, and ◻c, respectively. This study aims to obtain analytical formulations for the elastic strain and eigenstrain influence functions of the “1 + 2” partitions.

3.1 Elastic Strain Influence Functions

From the point of view of macroscopic elastic strain energy density, the following can be obtained:

𝒰c=12σc:εc=12εc:Lc:εc(7)

On the other hand, the mesoscale partitions give the elastic strain energy density that

𝒰uc=12cfσf:εf+12caσa:εa+12cbσb:εb=12cfεf:Lf:εf+12caεa:Lm:εa+12cbεb:Lm:εb(8)

Under elastic response, the averaged strains at partitions are given as follows:

εf=Ef:εc,εa=Ea:εc,εb=Eb:εc(9)

Substituting (9) into (8) yields

𝒰uc=12cfεc:Ef:Lf:Ef:εc+12caεc:Ea:Lm:Ea:εc+12cbεc:Eb:Lm:Eb:εc(10)

Comparing (7) and (10) and asking for 𝒰c=𝒰uc holds for arbitrary εijc yield

cfEf:Lf:Ef+caEa:Lm:Ea+cbEb:Lm:Eb=Lc(11)

In addition, recall the equations for the elastic strain influence functions (3) and (4) for the fibrous unitcell that

cfEf+caEa+cbEb=I⇒cfEf+cm(cacmEa+cbcmEb)=I(12)

cfLf:Ef+caLm:Ea+cbLm:Eb=Lc  ⇒cfLf:Ef+cmLm:(cacmEa+cbcmEb)=Lc(13)

cfEf=(Lf−Lm)−1:(Lc−Lm)(14)

In addition, define

cmEm≡(Lm−Lf)−1:(Lc−Lf)(15)

so that

caEa:Lm:Ea+cbEb:Lm:Eb=Lc−cfEf:Lf:Ef(16)

caEa+cbEb=cmEm(17)

cacb{Ea}T{Lm}{Ea}+[ca{Ea}−cm{Em}]T{Lm}[ca{Ea}−cm{Em}]=cb{Lc}−cbcf{Ef}T{Lf}{Ef}(18)

Rearranging (18) yields

[{Ea}−{Em}]T{Lm}[{Ea}−{Em}]=cb/cacm[{Lc}−cf{Ef}T{Lf}{Ef}−cm{Em}T{Lm}{Em}](19)

Similarly, substituting (17) into (16) and eliminating {Ea} yield

[{Eb}−{Em}]T{Lm}[{Eb}−{Em}]=ca/cbcm[{Lc}−cf{Ef}T{Lf}{Ef}−cm{Em}T{Lm}{Em}](20)

The right-hand-side term in (19) and (20) can be proved as semi-positive-definite. The proof is given in Appendix A. Suppose

{ΔL}={Lc}−cf{Ef}T{Lf}{Ef}−cm{Em}T{Lm}{Em}={R}T{R}(21)

where {R} is obtained through singular value decomposition (SVD). The first row and column of {ΔL} are all zeros. Defining that

{{ΔEa}={Ea}−{Em}{ΔEb}={Eb}−{Em}(22)

For any macroscopic strain {εc}=[ε11c,0,0,0,0,0]T, the following equation can always be fulfilled:

{εc}T{ΔL}{εc}=cmcacb{εc}T{ΔEa}T{Lm}{ΔEa}{εc}=cmcbca{εc}T{ΔEb}T{Lm}{ΔEb}{εc}=0(23)

Due to the positive definiteness of {Lm}, the first columns of {ΔEa} and {ΔEb} are all zeros. In addition, the strain along the longitudinal direction within the matrix should be uniform. That yields that the first rows of {Em}, {Ea}, and {Eb} are the same so that the first rows of {ΔEa} and {ΔEb} are also zeros. Then {ΔEa} and {ΔEb} can be written as follows:

{ΔEa}=[00T0ΔEa]{ΔEb}=[00T0ΔEb](24)

where ΔEa and ΔEb are 5×5 matrices to be solved. Finally, without losing generosity, it obtains the solution for {Ea} and {Eb} that

{{Ea}={Em}+cb/caca+cb{Vm}{R}{Eb}={Em}−ca/cbca+cb{Vm}{R}(25)

where cm=ca+cb is substituted and the detailed process of calculation for {Vm} and {R} is given in Appendix B.

3.2 Eigenstrain Influence Functions

For the fibrous unitcell with 1 + 2 partitions, the reduced-order system of equations is derived from (2) that

εf−Pff:μf−Pfa:μa−Pfb:μb=Ef:εc(26)

εa−Paf:μf−Paa:μa−Pab:μb=Ea:εc(27)

εb−Pbf:μf−Pba:μa−Pbb:μb=Eb:εc(28)

The eigenstrain influence tensors satisfy the following constraints (Chapter 4 of the textbook [52]):

Pff+Pfa+Pfb=I−Ef(29)

Paf+Paa+Pab=I−Ea(30)

Pbf+Pba+Pbb=I−Eb(31)

This constructs two special cases with perfect plastic yielding at selected partition(s) to solve the eigenstrain influence functions. For Case I, it asks that the outer partition of the matrix gets perfect plastic yielding while the inner partition of the matrix and the fiber phase remain elastic. For Case II, it assumes that the whole matrix phase gets perfect plastic yield while the fiber is still elastic. A sketch of two cases is depicted in Fig. 3.

Figure 3: Sketch of two cases with perfect plastic yielding at selected partition(s)

3.2.1 Case I

In the first case, under macroscopic strain εc, the strains at the fiber phase and the inner partition of the matrix phase are given as follows:

εf=[1−νIf−νIf000]T⋅ε11c≡[1e1f,I]T⋅ε11c(32)

εa=[1−νIa−νIa000]T⋅ε11c≡[1e1a,I]T⋅ε11c(33)

νIf=νatf+caφ(1−2νm)(νm−νatf)2cf(1−νm)+ca(1+φ(1−2νm))(34)

νIa=νm−cf(1−2νm)(νm−νatf)2cf(1−νm)+ca(1+φ(1−2νm))+cfcaln⁡cf+cacf(cf+ca)(νm−νatf)2cf(1−νm)+ca(1+φ(1−2νm))(35)

The derivation of the vector components can be seen in Appendix C.

Thus, it can define the so-called perfect plastic influence function for the fiber phase and the inner partition of the matrix phase that

{E~f,I}=[10Te1f,IO5](36)

and

{E~a,I}=[10Te1a,IO5](37)

{E~b,I}=[10Te1b,IE~b,I](38)

For the fiber-dominated mode, i.e., the longitudinal direction, e1b,I is determined by the constraint condition (3) that

e1b,I=1cb[−cfe1f,I−cae1a,I]=[−νIb−νIb000]T(39)

with

νIb=−cfνatf+caνmcb−cacfcb(φ−1)(1−2νm)(νm−νatf)2cf(1−νm)+ca(1+φ(1−2νm))−cfcbln⁡cf+cacf(cf+ca)(νm−νatf)2cf(1−νm)+ca(1+φ(1−2νm))(40)

The sub-matrix E~b,I defined in (38) is calculated by the effort of eliminating the so-called “inclusion-locking” phenomenon. With arbitrary macroscopic incremental strain Δεc with Δε11c=0, it asks for the contribution of the fiber phase and inner partition of the matrix phase to the macroscopic stress is 0, i.e.,

Lc:Δεc=cbLm:Δεb=cbLm:E~b,I:Δεc(41)

with the arbitrariness of Δεc with Δε11c=0, it can obtain E~b,I through (41) satisfying

{[1(l1c)Tl1cLc]−cb[1(l1m)Tl1mLm][10Te1b,IE~b,I]}[0Δε22Δε33Δε23Δε13Δε12]=[000000],∀Δεij,ij=22,33,23,13,12(42)

Once the perfect plastic influence functions for all partitions are determined, the eigenstrain influence functions Pbb, Pba, and Pbf can be computed one after another. The reduced-order model is supposed given by an incremental macroscopic strain Δεc. By the assumption that the fiber and the inner partition of the matrix phase remain elastic while the outer partition of the matrix phase gets yield, it has Δμf=Δμa=0, and Δεb=Δμb. Accordingly, (26)–(28) are rewritten as follows:

Δεf−Pfb:Δεb=Ef:Δεc(43)

Δεa−Pab:Δεb=Ea:Δεc(44)

Δεb−Pbb:Δεb=Eb:Δεc(45)

From (45), it can find

Δεb=(I−Pbb)−1:Eb:Δεc=E~b,I:Δεc⇒Pbb=I−Eb:(E~b,I)−1(46)

Through (43) and (44), it can obtain

Δεf=(Pfb:E~b,I+Ef):Δεc=E~f,I:Δεc⇒Pfb=(E~f,I−Ef):(E~b,I)−1(47)

and

Δεa=(Pab:E~b,I+Ea):Δεc=E~a,I:Δεc⇒Pab=(E~a,I−Ea):(E~b,I)−1(48)

respectively.

3.2.2 Case II

Case II is studied, where the whole matrix phase gets perfect plastic yield while the fiber is still elastic. Similarly, under macroscopic strain εc, the strains at the fiber phase and the inner partition of the matrix phase are given as follows:

εf=[1−νIIf−νIIf000]T⋅ε11c≡[1e1f,II]T⋅ε11c(49)

εa=[1−νIIa−νIIa000]T⋅ε11c≡[1e1a,II]T⋅ε11c(50)

νIIf=νatf(51)

νIIa=νm−cf(1−2νm)(νm−νatf)2cf(1−νm)+ca+cfcaln⁡cf+cacf(cf+ca)(νm−νatf)2cf(1−νm)+ca(52)

Thus, it can define the perfect plastic influence function for the fiber phase and the inner partition of the matrix phase that

{E~f,II}=[10Te1f,IIO5](53)

and

{E~a,II}=[10Te1a,IIO5](54)

The strain influence function for the outer partition of the matrix phase is first written as follows:

{E~b,II}=[10Te1b,IIE~b,II](55)

For the fiber-dominated mode, i.e., the longitudinal direction, e1b,II is determined by the constraint condition (3) that

e1b,II=1cb(−cfe1f,II−cae1a,II)=[−νIIb−νIIb000]T(56)

with

νIIb=−cfνatf+caνmcb+cacfcb(1−2νm)(νm−νatf)2cf(1−νm)+ca−cfcbln⁡cf+cacf(cf+ca)(νm−νatf)2cf(1−νm)+ca(57)

The sub-matrix E~b,II defined in (38) is calculated by the effort of eliminating the so-called “inclusion-locking” phenomenon. With arbitrary macroscopic incremental strain Δεc with Δε11c=0, it asks for the contribution of the fiber phase and inner partition of the matrix phase to the macroscopic stress is 0, i.e.,

Lc:Δεc=cbLm:Δεb=cbLm:E~b,II:Δεc(58)

with the arbitrariness of Δεc with Δε11c=0, it can obtain E~b,II through (58).

3.3 Volume Fraction of Inner and Outer Matrix Partitions

The purpose of searching for the optimum volume fraction of inner matrix partition ca is to minimize the average error of strain energy density within the matrix phase domain between 1 + 2 partition ROH model and Eshelby tensor solution in a given unit macroscopic strain loading space, written as ℒ. The average of squared error of strain energy density within matrix domain 𝒟¯u under given ca and macroscopic strain loading εc is expressed as follows:

𝒟¯u(ca,εc)=1Θm[∫Θa(𝒰m(x,εc)−𝒰¯a(ca,εc))2dΘ+∫Θb(𝒰m(x,εc)−𝒰¯b(ca,εc))2dΘ](59)

where

𝒰m(x,εc)=12εm(x,εc):Lm:εm(x,εc),x∈Θm(60)

𝒰¯a(ca,εc)=12εc:Ea(ca):Lm:Ea(ca):εc(61)

𝒰¯b(ca,εc)=12εc:Eb(ca):Lm:Eb(ca):εc(62)

𝒰m(x,εc) is the strain energy density distribution within the matrix domain. 𝒰¯a and 𝒰¯b are the average strain energy densities of inner and outer matrix partitions, respectively, obtained through the 1 + 2 partition ROH model. The strain distribution within the matrix phase domain εm(x) is calculated by an Eshelby tensor method for finite-domain problems [48]:

εm(x)=εm+εmD(x)=εm+(Sm(x)+SBF(x)):ε∗,x∈Θm(63)

Sf and Sm(x) (x∈Θm) are classical Eshelby tensors of the fiber domain and matrix domain, respectively, for cylindrical inclusion problem based on infinite matrix assumption, and SBF(x) is the modified tensor for the finite matrix domain condition. Their non-zero components are given in Appendix D. εmD(x) is the strain disturbance within the matrix phase. ε∗ is the uniform equivalent eigenstrain within the fiber domain induced by inhomogeneity. The study assumes that component ε11∗=0. As the uniform strain field within the fiber domain can be obtained directly by the elastic strain influence functions, the strain disturbance within fiber εfD and ε∗ can be calculated by the following equation:

(Sf+⟨SBF(x)⟩Θf):ε∗=εfD=εf−εm=[Ef:(Em)−1−I]:εm(64)

Then (64) yields ε∗=A:εm.

The average 𝒟¯u(ca,εc) within loading space ℒ is written as follows:

𝒟¯u(ca)=1Ωℒ∫Ωℒ𝒟¯u(ca,εc)|εc∈ΩℒdΩ(65)

Therefore, the optimum inner matrix partition volume fraction can be calculated by solving the following equation numerically to obtain the extreme point of 𝒟¯u(ca):

∂𝒟¯u(ca)∂ca=0(66)