A high-order three-scale reduced asymptotic approach for thermo-mechanical problems of nonlinear heterogeneous materials with multiple spatial scales

Journal Pre-proof A high-order three-scale reduced asymptotic approach for thermo-mechanical problems of nonlinear heterogeneous materials with multiple spatial scales Zhiqiang Yang, Chuanzhou Long, Yi Sun PII:

S0997-7538(19)30774-0

DOI:

https://doi.org/10.1016/j.euromechsol.2019.103905

Reference:

EJMSOL 103905

To appear in:

European Journal of Mechanics / A Solids

Received Date: 25 September 2019 Revised Date:

30 October 2019

Accepted Date: 12 November 2019

Please cite this article as: Yang, Z., Long, C., Sun, Y., A high-order three-scale reduced asymptotic approach for thermo-mechanical problems of nonlinear heterogeneous materials with multiple spatial scales, European Journal of Mechanics / A Solids (2019), doi: https://doi.org/10.1016/ j.euromechsol.2019.103905. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Masson SAS.

A

high-order

three-scale

reduced

asymptotic

approach

for

thermo-mechanical problems of nonlinear heterogeneous materials with multiple spatial scales Zhiqiang Yang∗, Chuanzhou Long, Yi Sun Department of Astronautic Science and Mechanics, Harbin Institute of Technology, Harbin, 150001 China

Abstract An effective high-order three-scale reduced asymptotic (HTRA) approach is developed to simulate the coupled thermo-mechanical problems of nonlinear heterogeneous materials with multiple spatial scales. In this paper, the heterogeneous structures are constructed by periodical layout of unit cells on microscopic and mesoscopic scale. At first, the linear and nonlinear higher-order local cell functions defined at microscale and mesoscale are derived by multiscale asymptotic expansion approach. Further, two types of homogenized parameters are introduced, and related nonlinear homogenized solutions are obtained on the macroscopic domain. Then, the high-order three-scale temperature and displacement solutions are established by assembling the various multiscale local solutions and homogenization solutions. The key characteristic of the presented methods is an efficient reduced form for evaluating higher-order local cell problems, and hence greatly reducing the calculating amount in comparison to direct numerical simulations. Also, the new high-order nonlinear homogenization solutions that do not require higher-order continuity of the macroscale solutions. Finally, the effectiveness and accuracy of the multiscale algorithms are verified by some typical numerical examples. Keywords HTRA algorithms; Nonlinear thermo-mechanical problems; Three-scale formulae; Reduced-order homogenization; Multiple spatial scales 1. Introduction Heterogeneous materials and structures with excellent elastic/inelastic properties, such as good thermal stability, high strength, high heat resistance, etc., have been extensively utilized for some practical engineering problems. In general, these kinds of composites are usually heterogeneous at a certain scale (Kuznetsov and Fish, 2012), and the corresponding effective information can be determined by microscale or mesoscale structure. In particular, owing to some materials usually having multiscale ∗

Corresponding author. E-mail address: [email protected] (Z. Q. Yang). 1

character, and therefore it is extremely essential to develop an efficient numerical approach for investigating the macroscopic properties of the heterogeneous materials under multiscale and multi-physics fields. As far as we know, classical homogenization approaches (or zeroth-order approaches) predict macroscopic performance of heterogeneous materials by changing the original differential equations with rapid oscillating parameters to the homogenization equations with efficient parameters (Cioranescu and Donato, 1999; Oleinik et al., 1992; Babuska, 1976; Babuska et al., 1974; Li et al., 2016; Patil et al., 2017; Fish et al., 1999; E and Engquist, 2003; Abdulle and Nonnenmacher, 2011; Allaire, 2003; Hou and Wu, 1997; Zhang et al., 2010; Farhat et al., 2001; Hughes, 1995; Zabaras and Ganapathysubramanian, 2009; Geers et al., 2001; Miehe, 2002; Geers et al., 2010; Kouznetsova et al., 2004; Fish and Yuan, 2007). This idea was firstly proposed by Babuska (Babuska, 1976; Babuska et al., 1974) and later heightened by some remarkable contributions, including the heterogeneous multiscale methods (E and Engquist, 2003; Abdulle and Nonnenmacher, 2011), multiscale finite element methods (Hou and Wu, 1997; Zhang et al., 2010), the discontinuous enrichment methods (Farhat et al., 2001) and variational multiscale methods (Hughes, 1995). Further, these theories of the homogenization approach were also expanded to several nonlinear problems (Fish et al., 1999; Geers et al., 2001; Kouznetsova et al., 2004; Fish and Yuan, 2007; Zhang and Oskay, 2016; Zhang and Oskay, 2017), such as the multiscale enrichment based on partition of unity approach (Fish and Yuan, 2007), variational multiscale enrichment methods (Zhang and Oskay, 2016; Zhang and Oskay, 2017) and computational homogenization methods (Geers et al., 2010). However, these multiscale approaches and related algorithms just considered the first-order asymptotic homogenization solutions for the problems. Under some conditions (Yu and Cui, 2007; Yang et al., 2012; Yang et al., 2015; Guan et al., 2015), the first-order homogenization solutions can not effectively get the local ocillaton information for the heterogeneous materials. To resolve the shortcomings, Fish and Kuznetsov (2010) developed a called as computational continua formula by changing the positions of quadrature points. Bourgat (1977), Bacigalupo (2014), Smyshlyaev and Cherednichenko (2000), Cui and co-workers (Yu and Cui, 2007; Yang et al., 2012; Yang et al., 2015; Guan et al., 2015), Boutin (1996), Allaire and Habibi (2013) studied the distinct higher-order multiscale asymptotic homogenization approaches. A main feature of these multiscale approaches is that higher-order continuity of macroscale 2

(or coarse-scale) solutions is usually avoided due to the sequential homogenized solutions for distinct orders and reestablishment of higher-order macroscale derivatives by postprocessing. Furthermore, the higher-order multiscale approaches are also developed to investigate the thermo-mechanical coupling problems of composites with periodic oscillation coefficients (Oleinik et al., 1992; Li et al., 2016). Temizer and Wriggers (2011) analyzed a thermoelasticity problem with small periodic distribution and supplied a detailed review for the theoretical analysis of the homogenized approaches. Terada et al. (2010) applied the two-scale asymptotic approach to study a thermo-mechanical problem considering microscale heat convection effect. Li et al. (2016) and Yang and Cui (2013) developed a second-order two-scale (SOTS) approach to investigate the different dynamic thermo-mechanical coupling problems, and provided an efficient multiscale algorithm. Yu and Fish (2002) developed the classical homogenization approach to simulate the thermo-viscoelastic coupling problems by introducing multiple spatial and temporal scales, and obtained enough reasonable results. For a detailed introduction of various thermo-mechanical problems, we can refer to (Li et al. (2016); Yang and Cui (2013)). Also, since some composite structures often possess different spatial scales of heterogeneities for practical engineering applications, such as hierarchical layered structures or concrete materials, the zeroth-order homogenization approaches and second-order homogenization approaches are not sufficiently accurate to analyze the effective thermal/mechanical properties of the heterogeneous materials. To solve the problems, Bensoussan et al. (2011) extended a reiterated homogenized approach to investigate the materials with multiple scales of heterogeneities. Allaire and Briane (1996), Holmbom et al. (2005) and Abdulle and Bai (2015) studied the mathematical theory for various three-scale problems. Fabricius et al. (2008) developed the reiterated homogenization methods to analyze hydrodynamic lubrication considering three different length scales. Guan et al., (2015), Zhang et al., (2015) and Chen et al. (2016) presented a random three-scale model to study effective performance of concrete materials by a three-scale multiscale algorithm. Later, Rodríguez et al. (2016) and Nascimento et al., (2017) developed the reiterated homogenized approach to predict effective performance of heat conduction problems for the composites with multiple configurations. Cao (2005) employed an iterated two-scale approach to analyze a linear elasticity problem with various spatial scales, and provided a mathematical proof for the problems. From Cao’s work, it is obviously found that the 3

iterative homogenized approach is an effective approach to compute the macroscopic coefficients from microscopic scale to macroscopic scale. Recently, Ramírez-Torres et al. (2018, 2019) used the iterated three-scale homogenization methods to study macroscopic performance of hierarchical composites with different spatial scales. Mahnken and Dammann (2016a, 2016b) applied a three-scale homogenization model for analyzing effective parameters of heterogeneous multiscale structure. Yang et al. (2017), Yang et al. (2018, 2019) and Dong et al. (2018) established a high-order multiscale homogenization method to simulate the heat and mechanical properties of the materials with multiple configurations. Further, Dong et al. (2019) developed the high-order three-scale homogenization methods to study the thermo-mechanical coupling problems of heterogeneous materials with multiple length scales. Obviously, by the linear higher-order correctors, the local microscale and mesoscale oscillating behavior of the heterogeneous structures can be captured more correctly. In addition, owing to some composite materials often manifesting nonlinear effect, such as damage, elastoplasticity, fatigue and so on, the zeroth-order homogenized approach will become quite difficult since the nonlinear local solutions must be calculated over and over again, and hence a large number of computation time is required to accurately get the local oscillation information of displacement fields and their derivatives (Fish et al., 1999; Geers et al., 2001; Kouznetsova et al., 2004). Therefore, the developments of reduced-order models for computing a nonlinear heterogeneous problem have become an important research area. Massively parallel computing (Mosby and Matous, 2016), order reduction simulation to multiscale problems or the coupling of both are introuduced for actual engineering computation. Also, on the basis of classical homogenization approaches, some order reduction methods, such as transformation field analysis (TFA) (Dvorak, 1992), fast Fourier transforms (FFT) (Moulinec and Suquet, 1998), the reduced-order homogenization methods for two scales (Oskay and Fish, 2007; Fish et al., 2019) and more than two-scales (Yuan and Fish, 2009a; Yuan and Fish, 2009b), and proper orthogonal decomposition (Oliver et al., 2017) among others were efficient in reducing the computation time for the non-linearity. In particular, Oskay and Fish (2007) constructed a reduced order model to analyze the failure problems and efficiently eliminate the computing time. Further, Zhang and Oskay (2016) and Zhang and Oskay (2017) extended a variational multiscale

method

to

investigate

the

elasto-viscoplastic

and

inelastic

thermo-mechanical behavior, respectively. Bhattacharjee and Matouš (2016) provided 4

a manifold-based reduced order method to analyze the hyperelastic composites. Fish et al. (2019) and Yang et al. (2019) extended a second-order reduction homogenization approach to the inelastic problems and nonlinear thermo-mechanical problems, respectively. According to the order reduction models and asymptotic homogenization techniques, the information of eigenstrains, including inelastic strains or moisture-induced strains, is studied by evaluating a local elasticity problem. However, the existing lower-order nonlinear homogenization models and related numerical techniques cannot be effectively applied to simulate the nonlinear thermo-mechanical problems for the heterogeneous composites with multiple spatial configuration owing to the complicated microstructure of the materials and strong coupling of the equations. In a word, there is short of adequate studies on the nonlinear coupled problems for the heterogeneous structures with multiple configuration. Besides, we have not discovered any efficient higher-order three-scale approaches to analyze the nonlinear problems for the composites. Thus, in this study, we principally discuss the thermo-mechanical coupling problems of nonlinear heterogeneous materials with multiple periodic configurations established by periodic distributions of unit cells on microscale and mesoscale (Fig.1). The main purpose of this work is to supply a unified micro-meso-macro multiscale model with efficient reduction format for evaluating the nonlinear materials with multiple configurations, which combines a higher-order homogenization approach that does not require higher-order continuity of macroscale solutions. The remainder of this work is listed as follows. Section 2 gives nonlinear governing equations of the heterogeneous materials with multiple spatial configurations. Section 3 describes the HTRA formulae for the coupled problems. In Section 4, an efficient model reduction form is presented. Section 5 is devoted to the detailed multiscale numerical algorithms. In Section 6, some representative examples are illustrated to confirm the effectiveness of the nonlinear multiscale algorithms. Finally, some conclusions are introduced in Section 7. 2. Nonlinear thermo-mechanical models of heterogeneous materials Here, the heterogeneous materials are described by three scales or two distinct local cells, as displayed in Fig.1. The macroscopic properties of the problems at distinct scales are studied by two unit cells: the microscopic cell ε 2 Z with size ε 2 and the mesoscopic cell ε1Y with size ε1 , and ε 2 < ε1 . The microscale area ε 2 Z includes 5

two or more material phases and builds the matrix on the mesoscale ε1Y . Also, the macroscale area Ω can be constructed by mesoscale area ε1Y and microscale area

ε 2Z . Then, the coupled thermo-mechanical equations for the heterogeneous materials with multiple spatial scales are introduced as follows:

(

)

(

∂ β ε1ε 2 (x)(T (x) − T ) − ∂ a ε1ε 2 (x)eε1ε 2 (u (x)) kl ε1ε 2 ε1ε 2 ∂x j ij ∂x j ijkl ε1ε 2 + ∂ ( aijkl (x) µklε1ε 2 (x) ) = f i (x), ∂x j

)

x ∈ Ω, i = 1, 2,3,

 ∂uε ε k (x) ∂uε1ε 2l (x)  eklε1ε 2 (x) ≡ uε1ε 2 ( k , xl ) (x) = 1  1 2 + , x ∈ Ω, 2  ∂xl ∂xk  µklε1ε 2 (x) = f (eklε1ε 2 (x), σ klε1ε 2 (x)), x ∈ Ω, ∂Tε ε (x) − ∂ (kijε1ε 2 (x) 1 2 ) = g (x), x ∈ Ω, ∂x j ∂xi

(

where Tε1ε 2 (x) and uε1ε 2 (x) = uε1ε 2 1 (x), uε1ε 2 2 (x), uε1ε 2 3 (x)

(1)

)

T

displacement fields; g (x) and f ( x) = ( f1 ( x), f 2 ( x), f 3 ( x) )

denote temperature and

T

describe thermal source

and body force; eklε1ε 2 (uε1ε 2 (x)) (k , l = 1, 2,3) denote the total strains including elastic and inelastic components; µklε1ε 2 (x) is called eigenstrains (i.e., inelastic strains, moisture effects or phase transformation (Yuan and Fish, 2009; Fish et al.(2019); Yang et al. (2019)). In this work, only the inelastic strains are investigated); kijε1ε 2 (x) , εε βijε ε (x) and aijkl (x) (i, j , k , l = 1, 2, 3) represent the heat conductivity, thermal 1 2

1 2

modulus and elastic tensor, and satisfy symmetry

kijε1ε 2 (x) = k εji1ε 2 (x) ,

βijε ε (x) = β εji ε (x) , 1 2

1 2

ε1ε 2 ε1ε 2 ε1ε 2 1ε 2 aijkl (x) = aεjikl (x) = aijlk (x) = aklij (x) ,

and ellipticity

λ1γ iγ i ≤ kijε ε (x)γ iγ j ≤ λ2γ j γ j , 1 2

λ1γ iγ i ≤ βijε ε (x)γ iγ j ≤ λ2γ jγ j , 1 2

εε λη 1 ijηij ≤ aijkl ( x)ηijη kl ≤ λ2η klη kl , 1 2

where λ1 and λ2 are independent of x , ε 1 and ε 2 , ηij (i, j = 1, 2, 3) are arbitrary symmetric matrix and γ i (i = 1, 2, 3) are arbitrary vector. 6

Fig. 1. Three-scale heterogeneous structure

Further, the boundary conditions are given in the following:

v jσ ijε1ε 2 (x) = pi (x),

x ∈ Γσ ,

uε1ε 2 (x) = u (x),

x ∈ Γu ,

Tε1ε 2 (x) = T% ,

x ∈ ΓT ,

vi kijε1ε 2 (x)

∂Tε1ε 2 (x) = 0, x ∈ Γ q , ∂x j

(2)

where v = (v j ), j = 1, 2,3 denote the exterior normal on Γσ and Γ q , p(x) are the tractions on Γσ , and u(x) denote the displacement values on Γu . T% describes the temperature results on ΓT , and satisfies ΓT I Γ q = ∅ , Γσ I Γu = ∅ , as illustrated in

εy Fig. 1. y = x and z = x = 1 , x denotes the macroscale coordinate, y is the

ε1

ε2

ε2

mesoscale coordinate in ε1Y and z denotes microscale coordinate in ε 2 Z . Then, the parameters

rapidly

kijε1ε 2 (x) = kij ( x , x ) = kij (y, z )

ε1 ε 2

oscillation ,

βijε ε (x) = βij ( x , x ) = βij (y, z ) and ε1 ε 2 1 2

ε1ε 2 aijkl (x) = aijkl ( x , x ) = aijkl (y, z ) are Y × Z -periodic piecewise differential functions.

ε1 ε 2

3. High-order three-scale formulae In this section, the high-order three-scale asymptotic formulae are derived for analyzing the nonlinear thermo-mechanical problems of the heterogeneous materials with multiple spatial scales. Motivated by the ideas given in (Yang et al., 2017; Yang et al., 2018; Dong et al. 7

2018; Dong et al., 2019), uε1ε 2 (x) and Tε1ε 2 (x) can be obtained in the following forms: uε1ε 2 (x) = u 0 (x) + ε1u1 (x, y ) + ε12u 2 (x, y ) + ε 2u3 (x, y , z )

(3)

+ε1ε 2u 4 (x, y , z ) + ε 22u 5 (x, y , z ) + ⋅⋅⋅ and Tε1ε 2 (x) = T0 (x) + ε1T1 (x, y ) + ε12T2 (x, y ) + ε 2T3 (x, y , z )

(4)

+ε1ε 2T4 (x, y , z ) + ε 22T5 (x, y , z ) + ⋅⋅⋅

The eigenstrain µijε1ε 2 (x) (i, j = 1, 2,3) depends on the constitutive behaviour of microscale

and

mesoscale

phases.

In

this

work,

we

will

assume

that

µijε ε (x) = µijf (x, y, z) , and they will be discussed in the following sections for details. 1 2

Let

A1 = −

∂ ∂ (aijkl (y , z ) ), ∂x j ∂xl

A2 = −

∂ ∂ ∂ ∂ (aijkl (y , z ) ) − (aijkl (y , z ) ), ∂y j ∂xl ∂x j ∂yl

A3 = −

∂ ∂ (aijkl (y , z ) ), ∂y j ∂yl

∂ ∂ ∂ ∂ A4 = − (aijkl (y , z ) ) − (aijkl (y , z ) ), ∂x j ∂zl ∂z j ∂xl A5 = −

∂ ∂ ∂ ∂ (aijkl (y , z ) ) − (aijkl (y , z ) ), ∂y j ∂zl ∂z j ∂yl

A6 = −

∂ ∂ (aijkl (y , z ) ), ∂z j ∂zl

(5)

and B1 = −

∂ ∂ (kij (y , z ) ), ∂xi ∂x j

B2 = −

∂ ∂ ∂ ∂ (kij (y , z ) )− (kij (y , z ) ), ∂yi ∂x j ∂xi ∂y j

B3 = −

∂ ∂ (kij (y , z ) ), ∂yi ∂y j

B4 = −

∂ ∂ ∂ ∂ (kij (y , z ) )− (kij (y , z ) ), ∂xi ∂z j ∂zi ∂x j

B5 = −

∂ ∂ ∂ ∂ ∂ ∂ (kij (y , z ) )− (kij (y , z ) ), B6 = − (kij (y , z ) ). ∂yi ∂z j ∂zi ∂y j ∂zi ∂z j

(6)

The three-scale differentiation operator is defined as: 8

∂ ∂ 1 ∂ 1 ∂ → + + . ∂x j ∂x j ε1 ∂y j ε 2 ∂z j Taking (3) and (7) into (1), and we have the following equality

(7)

∂ β ε1ε 2 (x)(T (x) − T ) − ∂ (a ε1ε 2 (x) ∂uε1ε 2 (x) ) + ∂ a ε1ε 2 (x) µ ε1ε 2 (x) ) kl ε1ε 2 ∂x j ij ∂x j ijkl ∂xl ∂x j ( ijkl

(

)

(

= ( ∂ + 1 ∂ + 1 ∂ ) βijε1ε 2 (x)(Tε1ε 2 (x) − T ) ∂x j ε1 ∂y j ε 2 ∂z j

)

−( A1uε1ε 2 (x) + 1 A2uε1ε 2 (x) + 12 A3uε1ε 2 (x)

ε1

(8)

ε1

+ 1 A4uε1ε 2 (x) + 1 1 A5uε1ε 2 (x) + 12 A6uε1ε 2 (x))

ε2

ε1 ε 2

ε2

ε1ε 2 + ( ∂ + 1 ∂ + 1 ∂ ) ( aijkl (x) µklε1ε 2 (x) ) = fi (x) ∂x j ε1 ∂y j ε 2 ∂z j

Also, substituting (4) and (7) into the second equation of (1) yields ∂Tε ε (x) − ∂ (kijε1ε 2 (x) 1 2 ) ∂xi ∂x j ∂Tε ε (x) 1 ∂Tε1ε 2 (x) 1 ∂Tε1ε 2 (x)   = −( ∂ + 1 ∂ + 1 ∂ )  kijε1ε 2 (x)( 1 2 + + ) ∂x j ε1 ∂y j ε 2 ∂z j  ∂x j ε 1 ∂y j ε 2 ∂z j  = − B T ( x ) + 1 B T ( x) + 1 B T ( x)

(

ε1

1 ε1ε 2

2 ε1ε 2

ε12

3 ε1ε 2

)

+ 1 B4Tε1ε 2 (x) + 1 1 B5Tε1ε 2 (x) + 12 B6Tε1ε 2 (x) = g (x)

ε2

ε1 ε 2

ε2

Then, the three-scale correction terms can be established in the following u1 (x, y ) = N p (y )

∂u 0 + P1 (y )(T0 − T ) + ∫ h(y , y% )µ1f (x, y% )dY% , Y ∂x p

u 2 (x, y ) = N pq (y )

∂ 2u 0 ∂T ∂µ f (x, y% ) % + Pα1 (y ) 0 + ∫ H p (y , y% ) 1 dY , ∂x p ∂xq ∂xα1 Y ∂x p

∂u ∂P (y ) ∂N p (y ) ∂u 0 u 3 ( x, y , z ) = M ( z ) + M h (z ) 0 + M h (z ) 1 (T0 − T ) ∂yh ∂x p ∂xh ∂yh h

+ P0 (z )T0 + M h (z ) ∂ ∂yh u 4 ( x, y , z ) = M h ( z )

∫ h(y, y% )µ Y

f 1

(x, y% )dY% + ∫ h 0 (y , z, z% )µ f (x, y , z% )dZ% , Z

2 ∂Pα1 (y ) ∂T0 ∂ 2u 0 ∂N pq (y ) ∂ u 0 + M h (z ) N p (y ) + M h (z ) ∂yh ∂x p ∂xq ∂xh ∂x p ∂yh ∂xα1

∂T ∂T +M (z )P1 (y ) 0 + P0 (z ) N α1 (y ) 0 + M h (z ) ∂ ∂xh ∂xα1 ∂yh h

+M h (z ) ∂ ∂xh

∫ h(y, y% )µ Y

f 1

(x, y% )dY% ,

9

∫

Y

H p (y , y% )

∂µ1f (x, y% ) % dY ∂x p

(9)

∂ 2u 0 ∂ 2 N pq (y ) ∂ 2u 0 hr u 5 ( x, y , z ) = M ( y , z ) + M (y, z ) ∂yh ∂yr ∂x p ∂xq ∂xh ∂xr hr

+ M hr (y , z )

∂N p (y ) ∂ 2u 0 ∂N p (y ) ∂ 2u 0 + M hr (y , z ) ∂yh ∂xr ∂x p ∂yr ∂xh ∂x p

+ M hr (y , z )(

2 ∂P1 (y ) ∂T0 ∂P1 (y ) ∂T0 ∂ Pα1 (y ) ∂T0 + + ) ∂yr ∂xh ∂yh ∂xr ∂yr ∂yh ∂xα1

+ M (y , z ) ∂ ∂yr

∂ µ f (x, y% )dY% + M hr (y , z ) ∂ 2 % h ( y , y ) ∫Y ∂xh 1 ∂yr ∂yh

hr

+ D h (y, z ) ∂ ∂yh

∫

Y

H p (y , y% )

+ D h (y, z ) N r (y ) + Eh (y , z )(

∫

Y

H p (y , y% )

∂µ1f (x, y% ) % dY ∂x p

∂µ1f (x, y% ) % dY + D h (y , z ) ∫ h(y , y% ) ∂ µ1f (x, y% )dY% Y ∂x p ∂xh

2 ∂ 2u 0 ∂N pq (y ) ∂ u 0 + D h (y, z ) ∂xh ∂xr ∂yh ∂x p ∂xq

∂Pα (y ) ∂T0 ∂T0 ∂N k (y ) ∂T0 ∂T + ) + G h (y , z )( 1 + P1 (y ) 0 ) ∂xh ∂yh ∂xk ∂yh ∂xα1 ∂xh

+ ∫ H 0h (y , z, z% ) ∂ µ f (x, y , z% )dZ% , Z ∂xh (10) where u 0 = u 0 (x) are the zeroth-order homogenization displacement. N (y ) , p

N pq (y ) , P1 (y ) , Pα1 (y ) , M h (z ) , P0 (z ) , M hr (y , z ) , Eh (y , z ) , G h (y , z ) and D h (y , z ) ( p, q, h, r =1,2,3) are the unit cell functions on the mesoscale and microscale domain, respectively. h 0 (y , z, z% ) , H 0h (y , z, z% ) , h(y , y% ) and H p (y , y% ) are so-called microscale and mesoscale transformation influence functions for the eigenstrain. For example, the physical meaning of

∫ h(y, y% )µ

f 1

Y

(x, y% )dY% is that the eigenstrain

µ1f ( x, y% ) provides elastic deformation of the magnitude of h(y , y% ) owing to volume or shape changing at an infinitesimal neighborhood of a point y% ∈ Θ (Yuan and Fish, 2009; Fish et al.(2019); Yang et al. (2019)). In addition, we have that

10

T1 (x, y ) = Nα1 (y )

∂T0 ∂ 2T0 , T2 (x, y ) = Nα1α 2 (y ) , ∂xα1 ∂xα1 ∂xα 2

T3 (x, y , z ) = M k (z )

∂Nα1 (y ) ∂T0 ∂T + M k (z ) 0 , ∂yk ∂xα1 ∂xk

T4 (x, y , z ) = M k (z )

∂Nα1α 2 (y ) ∂ 2T0 ∂ 2T0 , + M k (z ) Nα1 (y ) ∂yk ∂xα1 ∂xα 2 ∂xk ∂xα1

(11)

∂ 2 Nα1α 2 (y ) ∂ 2T0 ∂ 2T0 T5 (x, y , z ) = M kl (y , z ) + M kl (y , z ) ∂yk ∂yl ∂xα1 ∂xα 2 ∂xk ∂xl + M kl (y , z )

∂Nα1 (y ) ∂ 2T0 ∂Nα1 (y ) ∂ 2T0 ∂ 2T0 + M kl (y , z ) + Pk (y , z ) N l (y ) ∂yk ∂xα1 ∂xl ∂yl ∂xk ∂xα1 ∂xk ∂xl

∂Nα1α 2 (y ) ∂ 2T0 , + Pk (y , z ) ∂yk ∂xα1 ∂xα 2 where T0 = T0 (x) are the homogenized temperature solutions. Nα1 (y ) , Nα1α 2 ( y ) ,

M k (z ) , M kl (y , z ) , and Pk (y , z ) ( α1 , α 2 , k , l =1,2,3) denote the local cell functions on the mesoscale and microscale domains, respectively. And, they all can be solved in the following, successively. Taking

(10)

into

(1)

and

considering

the

different

values

of

ε 1mε 2n (m, n = −2, −1, 0,1, 2) yield to the following equalities: ∂ β ε1ε 2 (x)(T (x) − T ) − ∂ (a (y , z ) 1 ( ∂uε1ε 2 k (x) + ∂uε1ε 2l ( x) )) ε1ε 2 ∂x j ij ∂x j ijkl ∂xl ∂xk 2

(

)

+ ∂ ( aijkl (y , z ) µklε1ε 2 (x) ) = 12 A6u 0 + 1 1 A5u 0 + 12 A3u 0 + ε 1 12 A6u1 + ε12 12 A6u 2 ε1 ε 2 ∂x j ε2 ε1 ε2 ε2 + 1 ( A4u 0 + A5u1 + A6u3 + ∂ ( β ijε1ε 2 (x)T0 (x))) + 1 ∂ (aijkl (y , z ) µklf ) ε2 ε 2 ∂z j ∂z j +

ε1 ∂ ( ( β ε ε ( x)T1 ( x, y )) + A4u1 + A5u 2 + A6u 4 ) ε 2 ∂z j ij 1 2

(12)

+ 1 ( A2u 0 + A3u1 + A5u3 + ∂ ( βijε1ε 2 ( x)T0 ( x))) + 1 ∂ ( aijkl ( y , z ) µklf ) ε1 ε1 ∂y j ∂y j + ∂ ( β ijε1ε 2 ( x)(T0 ( x) − T ) ) + ∂ ( β ijε1ε 2 ( x)T1 ( x, y ) ) + ∂ ( β ijε1ε 2 ( x)T3 ( x, y, z ) ) ∂x j ∂y j ∂z j + ( A1u 0 + A2u1 + A3u 2 + A4u3 + A5u 4 + A6u5 ) + ∂ ( aijkl (y , z ) µklf ) + ⋅⋅⋅ = fi (x) ∂x j Obviously, from A6u 0 =0, A5u 0 =0, A3u 0 =0, A6u1 =0 and A6u 2 =0, we can find that u 0 (x) are independent of y and z , u1 ( x, y ) and u 2 ( x, y ) are independent of z . Further, from (12), the following equalities can be established:

O(ε 2−1 ) :

A4u 0 + A5u1 + A6u3 + ∂ ( β ijε1ε 2 (x)T0 (x)) + ∂ (aijkl (y , z ) µ klf ) = 0, (13) ∂z j ∂z j 11

O(ε1ε 2−1 ) :

A4u1 + A5u 2 + A6u 4 + ∂ ( βijε1ε 2 (x)T1 (x, y )) = 0, ∂z j

(14)

A2u 0 + A3u1 + A5u3 + ∂ ( βijε1ε 2 (x)T0 (x)) + ∂ (aijkl (y , z ) µ klf ) = 0, (15) ∂y j ∂y j

O(ε1−1 ) :

O(ε10ε 20 ) : A1u 0 + A2u1 + A3u 2 + A4u 3 + A5u 4 + A6u 5 + ∂ (aijkl (y , z ) µklf ) ∂x j + ∂ ( β ijε1ε 2 (x)(T0 − T ) ) + ∂ ( β ijε1ε 2 (x)T1 ) + ∂ ( β ijε1ε 2 (x)T3 ) = fi (x). ∂x j ∂y j ∂z j

(16)

Thus, (13)-(16) will be analyzed in detail to get the multiscale local functions, homogenized parameters and homogenized equations, respectively. At first, taking into account (10), (13) can be rewritten as:

(

∂aijnh (y , z ) ∂z j

+

p ∂M lnh ∂u0 n ∂N nm ∂u0 m ∂P1n ∂ (aijkl (y , z ) ))( + + (T0 − T ) ∂z j ∂zk ∂xk ∂yk ∂x p ∂yk

h   f % ) +  ∂ (a (y , z ) ∂P0 k ) − ∂β ij (y , z )  (T − T ) % % h ( y , y ) µ ( x , y ) dY 1 ijkl 0 ∫Y  ∂z j ∂zl ∂z j   ∂ ∂ ∂ + (aijkl (y , z ) (aijkl (y , z ) µklf ) = 0 h 0 (y , z, z% )µ f (x, y , z% )dZ% ) − ∫ Z ∂z j ∂zk ∂z j

+

∂ ∂yk

Then, by virtue of (17), the following equations are obtained:  ∂ ∂M lnh (z ) ∂aijnh (y , z ) , in Z , − ∂z (aijkl (y, z ) ∂z ) = ∂z j j k   h z ∈ ∂Z .  M ln (z ) = 0, and

 ∂  ∂P0 k (z )  ∂β ij (y, z ) , in Z ,  ∂z  aijkl (y, z ) ∂z  = ∂z l j   i  z ∈ ∂Z ,  P0 k (z ) = 0,

(17)

(18)

(19)

From Lax-Milgram theorem and Korn’s inequality, (18) and (19) have the unique solutions M lnh (z ) and P0 k (z ) (h, k = 1, 2,3) . Also, we have the equality for the eigenstrain functions given by

∂ ∂ (aijkl (y, z )( ∂z j ∂zk

∫

Z

h 0 (y, z, z% )µ f (x, y , z% )dZ% − µklf )) = 0 ,

(20)

and h0 (y , z, z% ) also satisfy the Dirichlet boundary conditions. Next, considering (10) and (14), the following equations which describe the

ε1ε 2−1 -term are formulated by:

12

pq  ∂aijnh (y , z ) ∂ ∂M lnh  p ∂ 2u0 m ∂N nm ∂ 2 u0 m ∂T + + + P1 0 ( a ( y , z ) ) ( N   nm ijkl ∂z j ∂z j ∂zk  ∂xh ∂x p ∂yh ∂x p ∂xq ∂xn  ∂ ∂µ1f (x, y% ) % f %+ ∂ % % % + h ( y , y ) µ ( x , y ) dY H ( y , y ) dY ) (21) 1 p ∂xh ∫Y ∂yh ∫Y ∂x p

 ∂ ∂β ij (y , z )  ∂P ∂T + (aijkl (y , z ) 0 k ) −  Nα1 0 = 0.  ∂z ∂xα1 ∂zl ∂z j   j Clearly, form (18) and (19), we can find that the equality (21) always hold. By analogy, and thanks to (10) and (15), the following equations can be obtained from the ε1−1 -term: ∂M lmp ∂u0 m ∂ (aijmp (y , z ) + (aijkl (y , z ) )) ∂y j ∂zk ∂x p +

p ∂M lnh ∂N nm ∂u ∂ ((aijnh (y , z ) + aijkl (y , z ) ) ) 0m ∂y j ∂zk ∂yh ∂x p

+

∂ ∂y j

 ∂M lnh ∂ ( ( y , z ) ( y , z ) ) a + a  ijnh ijkl ∂zk ∂yh 

∫ h(y, y% )µ

f 1

Y

∂ (aijkl (y , z ) M ln ) ∂ ∂ ∂ − (aijkl (y , z ) µ klf ) + ( ∂y j ∂zk ∂yh ∂y j

 (x, y% )dY%  

h

+

∂aijkl (y , z ) ∂ ∂y j ∂zl

∫

Z

∫ h(y, y% )µ Y

f 1

(x, y% )dY% )

∂ ∂ (aijkl (y , z ) h 0 (y , z, z% )µ f (x, y , z% )dZ% + ∂z j ∂yl

∫

Z

h 0 (y , z, z% )µ f (x, y , z% )dZ% )

p ∂M lnh ∂P1n  ∂  ∂ ∂  h ∂N nm  ∂u0 m + + ( M ln ) )  aijkl (y , z )  (aijnh (y , z ) + aijkl (y , z )  (T0 − T ) ∂z j  ∂yk ∂yh  ∂x p ∂y j  ∂zk ∂yh  ∂βij (y , z ) ∂aijkl (y , z ) ∂ ∂P1 (y ) ∂P (z ) ∂aijkl (y , z ) − (T0 − T ) + ( )(T0 − T ) + 0 k (T0 − T ) = 0 ∂y j ∂zl ∂yk ∂y j ∂zl ∂y j (22) Moreover, g (z ) represents an integrable function on microscale region

Z ∈ R n (n = 1, 2,3) and “ − ” represent the averaging operator defined by

1 g (z )dz . (23) Z ∫Z Appling the operator (23) to both sides of equality (22), and the homogenization g=

equations are gotten on the mesoscopic unit cell Y as follows:

13

1  aijmp (y ) ∂ ∂bij1 (y )  ∂N lmp  ∂u0 m  ∂ ∂P1k 1 1 + ( a ( y ) ) + ( a ( y ) ) +  (T − T )  ∂y j ∂y j ijkl ∂yk  ∂x p  ∂y j ijkl ∂yl ∂y j  0   1 + ∂ (aijkl (y ) ∂ ∫ h(y , y% )µ1f (x, y% )dY% ) (24) ∂y j ∂yl Y

+ ∂ ( ∫ aijkl (y , z )( ∂ ∂y j Z ∂zl

∫

Z

h 0 (y , z, z% )µ f (x, y , z% )dZ% − µ f ))dZ = 0

where 1 aijkl (y ) = 1 Z

bij1 (y ) = 1* Z

∫

Z

∫

Z

*

(aijkl (y , z ) + aijpq (y , z )

∂M plk )dz , ∂zq

(25)

∂P0 k (z )    aijkl (y, z ) ∂z − β ij (y, z ) dz , l  

(26)

1 and aijkl (y ) , bij1 (y ) (i, j , k , l = 1, 2, 3) are termed as the mesoscopic homogenized

parameters. Then, by virtue of (24) and the arbitrary homogenized solutions

∂u0m , ∂x p

N lmp (y ) and P1k (y ) should satisfy the following equations: 1  ∂ aijmp (y ) ∂N lmp (y ) 1 − ( a ( y ) ) = , in Y ,  ∂y ijkl ∂yk ∂y j j   p is Y- periodic,  N lm (y )

(27)

1  ∂  1 ∂P1k (y )  ∂bij (y ) a ( y ) = , in Y ,   ijkl ∂yl  ∂y j  ∂y j   P (y ) is Y- periodic,  1k

(28)

and h( y , y% ) denote the periodical function and satisfy the equalities described by ∂ (a 1 ( y ) ∂ ∂y j ijkl ∂yl

∫ h(y, y% )µ

f 1

Y

+ ∂ ( ∫ aijkl (y , z )( ∂ ∂y j Z ∂zl

∫

Z

(x, y% )dY% )

h 0 (y , z, z% )µ (x, y , z% )dZ% − µ f ))dZ = 0

(29)

f

Analogously, utilizing the operator (23) to both sides of (A.1), and we can obtain the following equation:

14

−a

∂ 2 u0 m ∂ 2 u0 m 1 p ∂ − (y ) ( a (y) Nlm ) ∂xk ∂x p ∂x j ∂x p ∂y j ijkl

−a

p pq ∂N nm ∂ 2 u0 m ∂N nm ∂ 2u0 m 1 ∂ − (y ) (a (y ) ) ∂yh ∂x j ∂x p ∂y j ijnh ∂yh ∂x p ∂xq

1 ijmp

1 ijnh

−bij1 (y )

∂T0 ∂P ∂T ∂T 1 1 − aijkl (y ) 1k 0 − ∂ (aijkl (y ) P1k ) 0 ∂x j ∂yl ∂x j ∂y j ∂xl

∂Pα k (y ) ∂T0 ∂T 1 − ∂ (aijkl − ∂ (bij1 (y ) Nα1 ) 0 (y ) 1 ) ∂y j ∂yl ∂xα1 ∂y j ∂x p − ∫ (aijkl (y , z )( ∂ Z ∂zl 1 − aijkl (y ) ∂ ∂yl

∫

Z

h 0 (y , z, z% ) ∂ µ f (x, y , z% )dZ% − ∂ µ klf ))dZ ∂x j ∂x j

∫ h(y, y% ) ∂∂x Y

1 − ∂ (aijkl (y ) ∂ ∂y j ∂yl

∫

Y

(30)

j

1 µ1f (x, y% )dY% − ∂ (aijkl (y ) ∫ h(y , y% ) ∂ µ1f (x, y% )dY% ) Y ∂y j ∂xl

H p (y , y% )

∂µ1f (x, y% ) % dY ) = fi (x). ∂x p

Similarly to (23), the integrable function g (y ) can be defined on mesoscopic domain Y ∈ R n (n = 1, 2,3) by

1 g (y )dy . (31) Y ∫Y Therefore, applying the operator (31) to both sides of equalities (30), and the g=

homogenized equations defined on macroscopic domain Ω are formulated as follows:

 ∂ 2 ∂  2 1 ∂u0 k (x) ∂u0 l (x)   ∂x ( bij (T0 (x) − T ) ) − ∂x  aijkl 2 ( ∂x + ∂x )  j  l k   j − ((a (y , z )( ∂ ∫ h 0 (y , z, z% ) ∂ µ f (x, y , z% )dZ% − ∂ µklf ))dZ )dY  ∫Y ∫Z ijkl ∂zl Z ∂x j ∂x j  % − a 1 (y ) ∂ h(y , y% ) ∂ µ f (x, y% )dYdY  ∫Y ijkl ∂yl ∫Y ∂x j 1  x ∈ Ω, = f i (x),  x ∈ Γu , u 0 (x) = u (x),  2 ∂u (x) = pi (x), x ∈ Γσ , v j aijkl 0 k ∂xl 

(32)

where 2 aijkl = 1 Y

bij2 = 1* Y

∫

Y

∫

Y

1 1 (aijkl (y ) + aijpq (y )

*

1 (bij1 (y ) + aijpq (y )

∂N plk (y ) )dy , ∂yq

(33)

∂P1 p (y ) )dy , ∂yq

(34)

2 2 , bij (i, j , k , l = 1, 2, 3) are defined as the macroscale homogenization and aijkl

15

parameters. Further, f i (x) (i = 1, 2,3) are substituted by (32), and (30) can be rewritten by: 1 − aijmp (y )

∂ 2 u0 m ∂ 2 u0 m 1 − ∂ ( aijkl (y ) N lmp ) ∂x j ∂x p ∂y j ∂xk ∂x p

1 − aijnh (y )

p ∂N nm ∂ 2 u0 m ∂N pq ∂ 2u0 m 1 − ∂ (aijnh (y ) nm ) ∂yh ∂x j ∂x p ∂y j ∂yh ∂x p ∂xq

−bij1 (y )

∂T0 ∂P ∂T ∂T 1 1 − aijkl (y ) 1k 0 − ∂ (aijkl (y ) P1k ) 0 ∂x j ∂yl ∂x j ∂y j ∂xl

∂Pα k (y ) ∂T0 ∂T 1 − ∂ (aijkl − ∂ (bij1 (y ) Nα1 ) 0 (y ) 1 ) ∂y j ∂yl ∂xα1 ∂y j ∂x p − ∫ (aijkl (y , z )( ∂ Z ∂zl 1 − aijkl (y ) ∂ ∂yl

= −a

Z

h 0 (y , z, z% ) ∂ µ f (x, y , z% )dZ% − ∂ µ klf ))dZ ∂x j ∂x j

∫ h(y, y% ) ∂∂x Y

1 (y ) ∂ − ∂ (aijkl ∂y j ∂yl 2 ijkl

∫

∫

Y

j

1 µ1f (x, y% )dY% − ∂ (aijkl (y ) ∫ h(y , y% ) ∂ µ1f (x, y% )dY% ) Y ∂y j ∂xl

H p (y , y% )

∂µ1f (x, y% ) % dY ) ∂x p

∂ 2u0 l ∂T + bij2 0 − ∫ ∫ ((aijkl (y , z )( ∂ ∂x j ∂xk ∂x j Y Z ∂zl

1 − ∂ µklf ))dZ )dY − ∫ aijkl (y ) ∂ Y ∂x j ∂yl

∫

Z

h 0 (y , z, z% ) ∂ µ f (x, y , z% )dZ% ∂x j (35)

f % ∫Y h(y, y% ) ∂∂x j µ1 (x, y% )dYdY

∂ 2u0m ∂T0 Owing to the arbitrariness of and in (35), the mesoscopic functions ∂xk ∂x p ∂x j N lmpq (y ) , Pα1k (y ) and H p (y , y% ) should satisfy the following problems:  ∂  1 ∂N lmpq (y )  2 1 y a ( )   ijkl  = aiqmp − aiqmp (y ) ∂ y ∂ y k   j  ∂N lmp (y )  ∂ 1 p 1 − ∂y ( aijlq (y ) N lm (y ) ) − aiqlk (y ) ∂y , in Y , j k  pq  N (y ) is Y- periodic.  lm  ∂  1 ∂Pα1k (y )  ∂ 1 aijk1 α1 (y ) P1k (y )   aijkl (y )  = −biα1 (y ) − ∂ y ∂ y ∂ y l j   j  ∂P1k (y ) ∂  1 1 2 in Y , − aiklα1 (y ) ∂y − ∂y bij (y ) Nα1 (y ) − biα1 , l j   Pα1k (y ) is Y- periodic.

(

(

(36)

)

)

16

(37)

∂µ1f (x, y% ) % ∫Y H p (y, y% ) ∂x p dY )

1 − ∂ (aijkl (y ) ∂ ∂y j ∂yl

= −∫

Y

∫

Z

((aijkl (y , z )( ∂ ∂zl

1 − ∫ aijkl (y ) ∂ Y ∂yl

∫

Z

h 0 (y , z, z% ) ∂ µ f (x, y , z% )dZ% − ∂ µklf ))dZ )dY ∂x j ∂x j

∫ h(y, y% ) ∂∂x Y

+ ∫ (aijkl (y , z )( ∂ Z ∂zl

∫

Z

% µ1f (x, y% )dYdY

(38)

j

h 0 (y , z, z% ) ∂ µ f (x, y , z% )dZ% − ∂ µ klf ))dZ ∂x j ∂x j

1 1 (y ) ∂ ∫ h(y , y% ) ∂ µ1f (x, y% )dY% + ∂ (aijkl (y ) ∫ h(y , y% ) ∂ µ1f (x, y% )dY% ) + aijkl Y ∂yl Y ∂x j ∂y j ∂xl Finally, the high-order three-scale homogenization solutions of thermo-mechanical

coupling problems (1) can be defined in the following:

u1ε1ε 2 (x) = u 0 (x) + ε1 ( N p (y )

∂u 0 + P1 (y )(T0 − T ) + ∫ h(y , y% )µ1f (x, y% )dY% ), Y ∂x p

uε21ε 2 (x) = u1ε1ε 2 (x) + ε12 ( N pq (y )

∂ 2u 0 ∂T ∂µ f (x, y% ) % + Pα1 (y ) 0 + ∫ H p (y , y% ) 1 dY ), (39) ∂x p ∂xq ∂xα1 Y ∂x p

u3ε1ε 2 (x) = uε21ε 2 (x) + ε 2u3 (x, y , z ) + ε 1ε 2u 4 (x, y , z ) + ε 22u5 (x, y , z ), T1ε1ε 2 (x) = T0 (x) + ε1 Nα1 (y )

∂T0 , ∂xα1

T2ε1ε 2 (x) = T1ε1ε 2 (x) + ε12 Nα1α 2 (y )

∂ 2T0 , ∂xα1 ∂xα 2

(40)

T3ε1ε 2 (x) = T2ε1ε 2 (x) + ε 2T3 (x, y , z ) + ε1ε 2T4 (x, y , z ) + ε 22T5 (x, y , z ), where u1ε1ε 2 (x) , T1ε1ε 2 (x) , uε21ε 2 (x) , T2ε1ε 2 (x) , uε31ε 2 (x) and T3ε1ε 2 (x) are termed as the first-order, the second-order two-scale solutions and high-order three-scale solutions, respectively. Besides, all the unit cell solutions are decided by equations (18) -(20), (27)-(29), (36)-(38), (A.3)-(A.8), (B.2), (B.3) and (B.7)-(B.9) in Appendix A and B. In

addition,

from

(39)

and

(40),

the

high-order

strain

fields

eijε1ε 2 (uε31ε 2 (x)) (i, j = 1, 2,3) and heat flux qiε1ε 2 (x) (i = 1, 2, 3) can be obtained by: eijε1ε 2 (uε31ε 2 (x)) = ε1ε 2

qi

ε1ε 2

(x) = −kij

1 ∂u3,ε1ε 2i (x) ∂u3,ε1ε 2 j (x) ( + ), 2 ∂x j ∂xi ∂T ε1ε 2 (x) ( x) 3 . ∂x j

4. Three-scale formulae of reduced order model 4.1. Formulae of the reduced unit cell models 17

(41)

First of all, considering (20), a reduced local cell function can be constructed by discretizing eigenstrains in the following (Yuan and Fish, 2009) Mz

µijf (x, y, z ) = ∑ N (α ) (z ) µij(α ) (x, y ), z

z

α z =1

µ

(α z ) ij

(x, y ) =

1 Z (α z )

Mz

Nα ∑ α (

z)

(z ) = 1,

z =1

(42)

µ (x, y , z )dZ ),

(∫

f ij

Z (α z )

and the partitions of unity condition for eigenstrain shape functions satisfy the piecewise constant approximation given by (α ) 1, z ∈ Z z , N (α z ) (z ) =  (43) (α ) 0, z ∉ Z z , in which the whole volumes of the microscale local cell are partitioned into M z

nonoverlapping subdomains described by Z (α z ) . Then, taking (42) and (43) into (20) gives the following equality

{a

ijkl

mic(α z ) mic(α z ) (y , z )[ Pklmn (y , z ) − I klmn (y , z )]}

,z j

= 0, ∀α z , y ∈ Y , z ∈ Z

(44)

where mic(α z ) Pijmn (y, z ) = ( ∫

Z (αz )

h0mni (y, z, z% )dZ% ), z j ,

(45)  I klmn , z ∈ Z (α z ) , I (y, z ) =  elsewhere. 0, Similarly, combining (42)-(45), a reduced mesoscale problem is gotten by discretizing mic(α z ) klmn

the eigenstrains for (29)

{a

1 ijkl

meso(α y ,α z )

(y )[ Pklmn

meso(α y )

(y ) − I klst

}

meso(α z ) (y ) S stmn (y )]

,yj

= 0, y ∈ Y ,

(46)

where meso(α y ,α z )

Pklmn S

(y ) = ∫

Y

meso(α z ) klmn

(α y )

(α ) meso(α z ) hijmy, yn (y, y% )Sijkl (y% )dY% ,

(y ) = −a (y )−1 ( ∫ 1 ijkl

Z (α z )

mic(α z ) aijkl (y , z )S mnkl (y, z )dZ ),

mic(α z ) mic(α z ) mic(α z ) Sklmn ( y, z ) = Pklmn (y, z ) − I klmn (y , z ),

I

meso(α y ) klmn

(47)

 I klmn , y ∈ Y (α y ) , (y ) =  elsewhere, 0,

and

µ1f (x, y ) =

µ

mic(α z ) mn

Mz

∑S α z =1

( x, y ) =

meso(α z ) klmn My

mic(α z ) (y ) µmn (x, y ),

(48)

∑N

(α y )

α y =1

18

(y ) µ

mic(α y ,α z ) mn

(x),

1, y ∈ Y (α y ) , N (y ) =  (α ) 0, y ∉ Y y , 1 mic(α ,α ) mic(α z ) µmn y z (x) = (α ) ( ∫ (α y ) µmn (x, y )dY ), y Y Y (α y )

(49) (50)

are all periodic functions. Further, from (32), (42)-(48), the order reduction equations for the macroscale stress are given as follows: My

Mz

2 σ ij (x) = aijkl ekl (x) + ∑ ∑ Aijkl

(α y ,α z )

α y =1 α z =1

mic(α y ,α z )

µkl

(x) − bij2 (T0 (x) − T ) ,

where

(51)

}

meso(α ,α ) meso(α ) 1 meso(α z ) = 1 ∫ {aijkl (y )[ Pklmn y z (y ) − I ijst y (y ) S stmn (y )] dY . (52) Y Y Thus, combining (33), (34), (51) and (52) yields to the homogenized equations in the (α ,α z )

Aijmny

following:

− ∂ σ (x) = f (x), x ∈ Ω i  ∂x j ij  v σ (x) = p (x), x ∈ Γ , i σ  j ij  u0 i (x) = ui (x), x ∈ Γu . (53) mic(α ,α ) In addition, owing to the arbitrariness of macroscale value µmn , x p y z (x) , and considering (42)-(48) and (52), the high-order influence functions (38) can be rewritten by:

−

∂ mn (α ,α ) meso(α z ) 1 1 (aijkl (y )Π kpl y z (y )) = −aipkl (y ) Sklmn (y ) ∂y j meso(α y ,α z ) klmn

+ a (y ) P 1 ipkl

∂ meso(α ,α ) mn (α ,α ) 1 (y ) + (aijkp (y )Θ kmn y z (y )) + ηip y z , ∂y j

(54)

where mn (α y ,α z )

Π kpl

(α ) meso(α z ) (y ) = ∂ ( ∫ (α y ) H kpijy (y , y% )S mnij (y% )dY% ), ∂yl Y

meso(α y ,α z )

Θ kmn

Y

mn (α y ,α z )

ηip

−∫

Y

(y ) = ∫

= −∫

Y

(α y )

(α y )

∫

Z

(α ) meso(α z ) hijk y (y , y% )S mnij (y% )dY% ,

meso(α y ,α z )

(α y )

1 aipkl (y ) Pklmn

}

mic(α z ) mic(α z ) a (y , z )[ Pklmn (y , z ) − I klmn (y , z )] dZdY ( α z ) { ipkl

(55)

(y )dY ,

and they also satisfy the periodic boundary conditions. Similarly to (44), (46) and (54), combining (42)-(48), (A.8) can be rewritten as follows:

19

∂ h (α z ) meso(α z ) 1 aijkl (y , z ) Frmkl (y , z ) = aihkl (y ) S klrm (y ) ∂z j

(

)

+ aihkl (y , z ) S

mic(α z ) klrm

∂ mic(α z ) (y , z )+ (aijkr (y , z )Qhmk (y , z )), ∂z j

(56)

where Fkh (α z ) (y , z ) = mic(α z ) imn

Q

∂ ∂zk

∫

Z (α z )

H 0h (y , z, z% )dZ% , (57)

(y , z ) = ∫ (α z ) h (y , z, z% )dZ% . mn 0i

Z

4.2. Order reduction system for the HTRA approach In this subsection, the order reduction systems of the thermo-mechanical problems can be constructed to simulate eigenstrains at each iteration for the nonlinear macroscale computation. From (41), the HTRA solutions of strain fields are given by

eklmic (x, y , z ) = eklmac (x) + u1, yl (x, y ) + u 3, zl (x, y , z ) + ε 1u1, xl (x, y ) + ε1u 4, zl (x, y , z ) +ε1u 2, yl (x, y ) + ε 2u 3, xl (x, y , z ) + ε 2u 4, yl (x, y , z ) + ε 2u5, zl (x, y , z ) +

ε2 ε2 u 3, y (x, y , z ) + 2 u 5, y (x, y , z ) + ε1ε 2u 4, x (x, y , z ) ε1 ε1 l

l

(58)

l

+ε12u 2, xl (x, y ) + ε 22u5, xl (x, y , z ), where eijmac (x) = 1 2 (ui 0, x j (x) + u j 0, xi (x)) . Further, combining (10), (45)-(48), and 1

applying the homogenized operator mic( β z ) ij

e

kl ( β z ) ij

( x, y ) = A

Z ( βz )

∫

Z ( βz )

⋅ dZ on both sides of (58) yields

Mz

meso kl

( y )e

mic( β z ~α z ) (x, y ) + ∑ Pijkl (y ) µklmic(α z ) (x, y )

α z =1

+ε1u1, x j (x, y ) + ε1u (4,βzz j) (x, y ) + ε1u 2, y j (x, y ) + ε 2u3,( βxz j) (x, y ) + ε 2u (4,βyz j) (x, y ) +ε 2 u

( βz ) 5, zl

ε 2 (β ) ε 22 ( β ) (x, y ) + u3, y (x, y ) + u 5, y (x, y ) + ε1ε 2u (4,βx ) (x, y ) ε1 ε1 z

z

j

(59)

z

j

j

+ε12u 2, x j (x, y ) + ε 22u5,( βxz j) (x, y ), where β z ~ α z denotes that partitions β z and α z belong to the same microscopic unit cell, and Aijkl ( β z ) ( y ) =

1 Z

( βz )

∫

Z

( βz )

( I ijkl ( y , z ) + M ijk, zl ( z ))dZ .

Further, we also have the following equality

20

(60)

eijmeso (x, y ) = eijmac (x) + u1, y j (x, y ) My

=A

meso ijkl

mac kl

(y )e

Mz

(x) + ∑ ∑ Pijkl

meso(α y ,α z )

α y =1 α z =1

mic(α y ,α z )

(y ) µkl

(61)

(x),

and meso Aijkl ( y ) = I ijkl ( y ) + N ikl, y j (y ).

(62)

Then, taking (61) and (62) into (59) gives My

Mz

meso(α y ,α z )

meso eijmic( β z ) (x, y ) = Aijkl ( β z ) (y ) Aijkl (y )eklmac (x) + Aijkl ( β z ) (y ) ∑ ∑ Pijkl

α y =1 α z =1

Mz

(α y )

mic( β z ~α z ) + ∑ Pijkl (y ) N

α z =1

mic(α y ,α z )

(y ) µ kl

mic(α y ,α z )

(y ) µkl

( x)

(x) + ε1u1, x j (x, y ) + ε1u (4,βzz j) (x, y ) + ε1u 2, y j (x, y )

(63)

ε 2 (β ) ε2 u3, y (x, y ) + 2 u5,( βy ) (x, y ) ε1 ε1

+ε 2u3,( βxz j) (x, y ) + ε 2u (4,βyz j) (x, y ) + ε 2u5,( βzzl ) (x, y ) +

z

z

j

j

+ε1ε 2u (4,βxz j) (x, y ) + ε12u 2, x j (x, y ) + ε 22u5,( βxz j) (x, y ) where

 M z mic( β z ~α z ) (y ), for β z ~ α z ,  ∑ Pijkl (α ) mic( β z ~α z ) Pijkl (y ) N y (y ) = α z =1 ∑ α z =1 0, otherwise,  Mz

(64)

and mic(α z )  1 ( y , z )dZ , β z , α z ∈ same Z ,  ( β z ) ∫Z ( β z ) Pijkl P (y ) =  Z (65) 0, else. (β ) Meanwhile, integrating the equation (63) over partition domain Y y (where β z mic( β z ~α z ) ijkl

partition is embedded in β y ) yields to the equations for the following partitioned microscale strain mic( β y , β z )

eij

( x) = ∫

Y

(βy )

My

kl ( β z ) meso Amn (y ) Amnkl (y )dYeijmac (x) Mz

+∫

kl ( β z ) (y ) ∑ ∑ Pmnkl ( β y ) Amn

+∫

(βy )

Y

Y

meso(α y ,α z )

α y =1 α z =1

Mz

∑P α z =1

mic( β z ~α z ) ijkl

(y )N

( β ,β z )

(β )

+ε1u1, xyj (x) + ε1u 4, zyj ( β ,βz )

+ε 2u5, zyl

(β )

( x) +

(α y )

mic(α y ,α z )

(y )dY µij

mic(α y ,α z )

(y )dY µ kl

( x)

(β )

( β ,βz )

(x) + ε 1u 2, yy j (x) + ε 2u3, xyj

( x)

(66) ( β ,βz )

(x) + ε 2u 4, yy j

ε 2 ( β ,β ) ε 2 ( β ,β ) ( β ,β ) u3, y (x) + 2 u 5, y (x) + ε1ε 2u 4, x (x) ε1 ε1 y

z

y

j

( β ,β z )

+ε12u 2, xyj (x) + ε 22u5, xyj

z

j

( x)

where 21

y

z

j

( x)

(β )

u1, xyj (x) =

1

Y

(β )

u 2, yy j (x) =

(β y )

1

Y

(βy )

( β ,βz )

∫

(βy )

u1, x j (x, y )dY , u 4, zyj

∫

(βy )

u 2, y j (x, y )dY , u3, xyj

Y

( β ,βz )

Y

( β ,βz )

u 4, xy j

( x) =

(β )

u 2, xy j (x) = ( β ,βz )

u5, zyj

( β ,βz )

u 4, yy j

( β ,βz )

u3, yy j

( β ,βz )

u5, xyj

Y

( x) = ( x) =

( β ,βz )

Y 1

( x) =

u5, yy j

1

( x) =

( x) =

1

(β y )

Z ( βz )

∫

(βy )

Y

1

Y

(β y )

1

Y Y

(βy )

1

Z ( βz ) 1

Z ( βz )

1

1 Z ( βz )

1

Y

1

Y

1

(β y )

(β y )

Y

(βy )

1

( x) =

Y

(βy )

∫ ∫ Y

(βy )

Y ( βz )

1

Y

( βz )

1

Z ( βz )

∫ ∫ Y

(βy )

Y ( βz )

u 4, z j (x, y , z )dZdY ,

∫ ∫ Y

(βy )

Z ( βz )

u 3, x j (x, y , z )dZdY ,

u 4, x j (x, y , z )dZdY ,

u 2, x j (x, y )dY ,

(β y )

(β y )

1

( x) =

Z ( βz ) 1 Z ( βz )

∫ ∫

u5, z j (x, y , z )dZdY ,

∫ ∫

u 4, y j (x, y , z )dZ dY ,

∫ ∫

u3, y j (x, y , z )dZdY ,

Y

Y

Y

(βy )

(βy )

(βy )

Z ( βz )

Z ( βz )

Z ( βz )

∫ ∫

u5, y j (x, y , z )dZdY ,

∫ ∫

u5, x j (x, y , z )dZdY .

Y

Y

(βy )

(βy )

Z ( βz )

Z ( βz )

(67)

Therefore, the order reduction constitutive relations for the eigenstrains have the following forms mic(α y ,α z )

µij

mic(α ,α )

mic(α ,α )

(x) = f (eij y z ( x), σ ij y z (x)). (68) Obviously, equations (66) and (68) are composed of the order reduction systems for the independent unknowns values

mic(α y ,α z )

eij

(x) (i, j = 1, 2,3) . The partitioned

eigenstrains (66) are gotten by computing a system of equations (66)-(68) within each macro increment. Once the partitioned eigenstrain is obtained, (51) and (53) are used to update the macroscopic stress and displacement. 5. Finite element (FE) algorithms for the HTRA approach The HTRA algorithms mainly include three steps: (i) pre-processing, (ii) predictor-corrector nonlinear iterative procedure, and (iii) post-processing. Step 1: Pre-processing 1) Verify the distributions of inclusions and geometric structure at the microscale, mesoscale and macroscale area. Further, the different areas are partitioned into FE meshes. 22

2) Compute the FE solutions of M k (z ) , M h (z ) and P0 (z ) (k , h = 1, 2, 3) from the microscale unit cell problems (B.2), (18) and (19) on Z . Also, the mesoscale

homogenized

parameters

kij1 (y )

,

1 aijkl (y )

and

bij1 (y ) (i, j , k , l = 1, 2,3) are obtained by solving the equations (B.4), (25) and (26). 3) Evaluate Nα1 (y ) , N p (y ) and P1 (y ) (α1 , p = 1, 2, 3) from the mesoscale unit cell problems (B.3), (27) and (28) on Y . Then, the macroscale 2 homogenized parameters kij2 , aijkl and bij2 (i, j , k , l = 1, 2,3) are solved by

the equalities (B.6), (33) and (34). 4) Compute the high-order local unit cell problems (36), (37), (A.3)-(A.7) and (B.7)-(B.9) with the FE methods. 5) Establish the weak forms for the partitioned eigenstrain transformation influence functions (44), (46), (54) and (56), and obtain their finite element solutions. 6) Compute the FE solutions T0 (x) by the homogenized problems (B.5) on entire domain Ω . Step 2: Predictor-corrector iterative HTRA algorithms Loop over all macroscale load increments and iterations at each load increment until the systems converge. From the converged solutions obtained by the previous load increment, then we implement the following steps: 1) (predictor) Considering the zeroth-order reduced homogenization problems (equations (51)-(53)), we can construct the FE forms of the homogenized problems (53), and evaluate the macroscale nodal displacement ∆u ( x ) and strain ∆ekl (x) increments from the former converged solutions. 2) (predictor) Giving the predictive results of ∆ekl (x) to evaluate the incremental reduced-zeroth order unit cell equations (66)-(68) for ε1 = 0 , ε 2 = 0 and mic(α y ,α z )

evaluate the increment of eigenstrains ∆µkl 3)

mic(α y ,α z )

∆ekl (x) and ∆µkl

( x) .

(x) are interpolated from quadrature points to each FE

nodes at macroscale meshes. Compute the first-order derivatives ∆ekl , xm ( x) and 23

mic(α y ,α z )

∆µkl , xm

(x) by the derivatives of the FE shape functions. Further, ∆ekl , x ( x) m

mic(α y ,α z )

and ∆µkl , xm

(x) are interpolated from the quadrature points to each FE nodes

for the macroscale meshes. Also, evaluate the second-order derivatives mic(α ,α z )

∆ekl , xm , xn ( x ) and ∆µkl , xm , xyn

∂ 2T0 Similarly, ∂xk ∂xα1

(x) by the derivatives of FE shape functions.

can be also obtained by the derivatives of FE shape

functions. 4) (corrector) From macroscale solutions obtained from predictor phase, evaluate the HTRA solutions, equations (63)-(68), and compute the macroscale residual. If the systems are converged, go to the next load increment. Otherwise, go to the next iteration. Step 3: Post-processing mic(α y ,α z )

1) From the obtained solutions ( ekl , ekl , xm , ekl , xm , xn , µkl

mic(α y ,α z )

, µkl , xm

mic(α ,α z )

, µkl , xm , xyn

) and

∂ 2T0 , we can interpolate them to each of microscopic and mesoscopic cell FE ∂xk ∂xα1 mesh. 2) The first-, second-order and HTRA solution over whole domain Ω can be assembled by equations (39) and (40), respectively. 6. Numerical examples In this section, the HTRA approach proposed in this work is validated by some typical numerical examples. Since it is very difficult to obtain theoretical solutions of the thermo-mechanical coupling problems (1) for the periodic materials, we should replace uε1ε 2 (x) and T ε1ε 2 (x) by direct numerical simulation (DNS) with a fine mesh. Both elastic and hyperelastic materials are listed to confirm the presented approach. For all problems, one partition per phase is adopted. 6.1. Hyperelastic materials A heterogeneous structure Ω is composed of entire periodical unit cells which is illustrated in Fig.2 (a), and the mesoscopic cell Y and microscopic cell Z are given in Figs.2 (b) and (c) , respectively. The length and width of the rectangular inclusions in Y are 1/3, and the radius of the circular inclusions in Z is 0.3. Also, ε 1 = 1 6

24

and ε 2 = 1 36 . Set f (x) = (0, −400)T , g (x) = 100 J/m3, pi (x) = 0 , u(x) = (0, 0)T and T = T% = 0 K . Material3

Material1

(a) Fig.2. (a) Macroscale domain Ω ;

Material2

(b) (c) (b) Mesoscale cell Y ; (c) Microscale cell Z

Herein, the linear triangular element is taken for the original problems (1) by a fined mesh and the relevant homogenization problems with coarse meshes. The materials are assumed to be hyperelastic with the following constitutive relations (Ming, 2002): εε εε εε σ ijε ε (x) = aijkl ( x)eklε ε (x) + Bijklmn (x)eklε ε (x)emn (x) − β ijε ε (x)(Tε ε ( x) − T ) , 1 2

1 2

1 2

1 2

1 2

1 2

1 2

1 2

ε1ε 2

where Bijklmn (x) are the higher-order coefficients of elasticity tensor. The Young’s moduli and Poisson's ratios for the three partitions are E(1) = 10 , E(2) = 100 ,

E(3) = 1000 and ν (1) = 0.3 , ν (2) = 0.2 , ν (3) = 0.2 . The heat conduction and heat modulus

are

β (3) = 0.1. The

k (1) = 10 ,

k (2) = 0.01 ,

quadratic

forms

k (3) = 0.1

for

the

β (1) = 10 ,

and

matrix

of

the

β (2) = 1 ,

materials

are:

(1) (1) 2 σ 11(1) = σ 22(1) = ((e11(1) ) 2 + 2e11(1) e22 + (e22 ) ), σ 12(1) = 0 .

In this work, TDNS and uDNS are adopted as the reference solutions by the DNS for the nonlinear coupled problems. For the DNS, we have the FE meshes comprising 86400 degrees-of-freedom. The mesoscale cell and microscale mesh consist of 790 and 832 degrees-of freedom. For the homogenized models, four different FE meshes, consisting of 5000, 800, 200 and 72 degrees-of-freedom are chosen. The numerical errors evaluated by various multiscale models about the DNS, i.e. u DNS − u 0 (x), u DNS − uε21ε 2 (x), u DNS − uε31ε 2 (x) , in L2 and H 1 -norms are defined as:

v

2

L2 ( Ω )

= (∫ v dx)1 2 , Ω

v

= ( ∫ ( ∇v )dx)1 2 , 2

H1 (Ω) 2

Ω

1

and normalized with regard to the related L or H -norm of u DNS . As an example, 25

the relative errors in the L2 -norm for the zeroth-order, second-order and HTRA εε approaches are defined as u DNS − u I1 2 (x)

L2

uDNS

L2

with I=0, 2, 3 for different

multiscale models. Thus, the relative errors in the L2 and H 1 -norms for the three distinct models are introduced in Fig. 3 and Fig. 4, respectively.

(a) (b) 2 1 Fig. 3. Relative errors of temperature distribution in the L and H -norm for different models

(a) (b) 2 1 Fig. 4. Relative errors of displacement distribution in the L and H -norm for different models

Fig. 5 describes the numerical results of T0 (x) , T2ε1ε 2 (x) , T3ε1ε 2 (x) and TDNS for this example.

26

(a)

(b)

(c)

(d) ε1ε 2

Fig.5. Temperature distributions (a) T0 ( x) ; (b) T2

(x) ; (c) T3ε1ε 2 (x) ; (d) TDNS

Fig. 6 describes the numerical results of u 0 (x) , uε21ε 2 (x) , uε31ε 2 (x) and u DNS for this example.

(a)

(b)

27

(c) (d) εε εε Fig. 6. The first parts of displacement distributions: (a) u0 (x) ; (b) u 21 2 (x) ; (c) u31 2 (x) ; (d)

uDNS Further, the relative errors of the homogenization, second-order and HTRA solutions in the L2 and H 1 -norms are clearly illustrated in Figs. 7 and 8 as a function of ε 1 .

(a)

(b)

Fig.7. Relative error for the temperature distributions in L2 and H1-norms as a function

ε1 .

(a) (b) 2 Fig.8. Relative error for the displacement distributions in L and H1-norms as a function

ε1 .

28

According to the numerical results given in Figs. 3-8, the zeroth-order homogenization approach and second-order reduced asymptotic (SORA) approach have not good correctness, but the HTRA approach demonstrates excellently and clearly provides a better approximation of the FE solutions than the results gotten by zeroth-order homogenization approach and SORA approach. In addition, it should be noted that the SORA approach can also effectively catch the mesoscale oscillating information for the nonlinear thermo-mechanical problems, however, the microscale effect of the heterogeneous materials cannot be captured accurately. In short, all the computational results shown previously prove that the HTRA algorithms are efficient and valid in simulating the nonlinear problems with multiple spatial scales. 6.2. L-shape elastic materials An L-shape region with a multiple periodic microstructure is considered in Fig. 9. The microscale and mesoscale structures are the same as introduced in the previous example 1. The body force is f (x) = (0, −400)T , g (x) = 100 J/m3, pi (x) = 0 and T = T% = 0 K . Also, ε 1 = 1 8 and ε 2 = 1 48 , and the boundaries of the L-shape region are constrained.

Fig. 9. L-Shape domain Ω

The Young’s moduli and Poisson's ratios for three partitions are E(1) = 10 , E(2) = 100 ,

E(3) = 1000 and ν (1) = 0.3 , ν (2) = 0.2 , ν (3) = 0.2 , respectively. The heat conductivity and heat modulus are k (1) = 10 , k (2) = 0.05 , k (3) = 1 and β (1) = 10 ,

β (2) = 1 , β (3) = 0.1. Also, for the DNS, we use a FE mesh comprising 115200 degrees-of-freedom. Similarly, we consider four FE meshes for the macroscale models, consisting of 4360, 1508, 246 and 102 degrees-of-freedoms. The relative 29

errors in the L2 and H 1 -norms for the three different models are depicted in Fig. 10 and Fig. 11.

(a) (b) 2 1 Fig. 10. Relative errors of temperature distribution in the L and H -norm for different models

(a) (b) 2 1 Fig. 11. Relative errors of displacement distribution in the L and H -norm for different models

Fig. 12 describes the computational results for T0 (x) , T2ε1ε 2 (x) , T3ε1ε 2 (x) and TDNS in this example.

30

(a)

(b)

(d)

(d) εε

εε

Fig. 12. Temperature distribution (a) T0 ( x) ; (b) T2 1 2 ( x) ; (c) T3 1 2 ( x) ; (d) TDNS

Fig. 13 describes the distribution of the first displacement components evaluated by the different approaches.

(a)

(b)

31

(c) (d) εε Fig. 13. The first components of the displacement distributions: (a) u 0 ( x) ; (b) u 21 2 (x) ; εε

(c) u31 2 (x) ; (d) u DNS

Obviously, the HTRA solutions can again give the excellent approximations for the DNS solutions. 6.3. 3D heterogeneous materials We consider a 3D heterogeneous structures with multiple spatial scales shown in Fig.14(a), and the mesoscopic and microscopic regions are described in Fig.14(b) and Fig.14(c), respectively. The size of the long axe of the ellipsoidal inclusions is chosen as 0.5, middle axe 0.36 and short axe 0.2. The length, width and height of the rectangular inclusions in Y are 0.5, and the radius of the spherical inclusions in Z is 0.3. Besides, the linear tetrahedral element is chosen in this example for the nonlinear multiscale problems, and the materials parameters for the problems are same with that listed in Section 6.2.

32

(a) Fig. 14. (a) Domain Ω ;

(b) (b) Mesoscale cell Y ;

(c) (c) Microscale cell Z

The mesoscopic and macroscopic effective heat conductivity parameters are proposed as follows:

−1.10178e− 5 4.70109 e− 5   8.4468   k (y ) =  −1.10178e− 5 8.4468 3.02193e− 5   4.70109 e− 5 3.02193e− 5 8.4467   1 ij

and

−3.77252 e− 5 −3.49964 e− 5   7.1571   k =  −3.77252 e− 5 7.1573 −5.1741e− 5   −3.49964 e− 5 −5.1741e− 5 7.1573   Then, the macroscopic effective heat conductivity parameters are provided for 2 ij

the ellipsoidal inclusions listed by

2.2475e− 6 2.7491e− 6   8.2843   k =  2.2475e− 6 8.1916 1.4421e− 5   2.7491e− 6 1.4421e− 5 8.2666   In addition, the mesoscopic and macroscopic effective elastic parameters are given 2 ij

in the following: 0 0 0   16.121 6.4197 6.4197   0 0 0   6.4197 16.121 6.4199  6.4197 6.4199 16.118 0 0 0  1 aijkl (y ) =   0 0 4.7851 0 0   0  0 0 0 0 4.7846 0    0 0 0 0 4.7848   0 33

and 0 0 0   20.688 5.71438 5.71369   0 0 0   5.71438 20.687 5.71272  5.71369 5.71272 20.689 0 0 0  2 aijkl =  0 0 6.890 0 0   0  0 0 0 0 6.891 0    0 0 0 0 6.889   0 Then, the corresponding mesoscopic and macroscopic effective Young’s moduli are E1= E2=E3=12.4688 and E1= E2=E3=18.2149. The macroscopic effective elastic parameters are evaluated for the ellipsoidal inclusions by 0 0 0   15.8769 5.12921 5.14957   0 0 0   5.12921 15.3724 5.12497  5.14957 5.12497 15.62 0 0 0  2 aijkl =  0 0 0 5.16494 0 0    0 0 0 0 5.27528 0    0 0 0 0 0 5.15618   And the macroscale effective Young’s moduli obtained for the ellipsoidal inclusions are E1=13.3156, E2=12.8535 MPa and E3=13.0719 MPa. Besides, the macroscopic effective Young’s moduli are different in the diagonal owing to the ellipsoidal inclusions containing at the mesoscopic cell Y . From the numerical results introduced in this example, the HTRA solutions are efficient to predict the macroscale properties of the heterogeneous materials with multiple spatial configurations. Especially, the different inclusions at the mesoscopic region can also be effectively predicted by the HTRA approach discussed in this work. 7. Conclusions A novel HTRA algorithms have been established for predicting the nonlinear thermo-mechanical properties of a composite with multiple configurations. The key features of the presented approaches are: (i) an efficient reduced form for analyzing high-order nonlinear local cell functions with less computational time and (ii) a new high-order nonlinear homogenization solution that does not require higher-order continuity of coarse-scale solutions. Finally, the correctness of the three-scale model and the availability of the HTRA approach are confirmed with the DNS solutions. From the computational results, it is well illustrated that the HTRA model is efficient to simulate the nonlinear 34

thermo-mechanical coupling problems by solving the new high-order corrections and homogenized parameters. Numerical results clearly show that the results gotten by the HTRA approach agree well with those by the DNS solutions, and are quite accurate than the zeroth-order and the SORA approaches. Acknowledgements This work is supported by the National Natural Science Foundation of China (11701123) and also supported by the Fundamental Research Funds for the Central Universities (HIT.NSRIF.2020017). Appendix A: High-order formulae for the HTRA approach Substituting (10) into (16) and considering (7) yield to the following equality

−

∂u ∂N lmp ∂u0 m ∂N lmpq ∂ 2u0 m ∂ ∂ ∂ (aijkl 0l ) − (aijkl )− (aijkl ) ∂x j ∂xk ∂x j ∂yk ∂x p ∂y j ∂yk ∂x p ∂xq

2 p p 2 ∂M lnh ∂N nm ∂u0 m ∂ ∂ ∂ p ∂ u0 m h ∂N nm ∂ u0 m − (aijkl N lm )− (aijkl )− (aijkl M ln ) ∂y j ∂x p ∂xk ∂x j ∂zk ∂yh ∂x p ∂z j ∂yh ∂xk ∂x p

−

pq ∂M lmh ∂u0 m ∂ 2u0 m ∂M lnh ∂N nm ∂ 2u0 m ∂ ∂ ∂ (aijkl )− (aijkl M lmh )− (aijkl ) ∂x j ∂zk ∂xh ∂z j ∂xk ∂xh ∂y j ∂zk ∂yh ∂x p ∂xq

pq  ∂ 2 u0 m  ∂ ∂M lnh p ∂ 2u0 m ∂ h ∂N nm a ( M ) − ( a N nm )  ijkl  ln ijkl y y x x y z x x ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ k h p q  j k h p  pq ∂ 2 u0 m  ∂ ∂M lnhr ∂ 2 N nm ∂ 2u0 m ∂  ∂ p − M lnh N nm (aijkl )  aijkl  − ( ) ∂zk ∂yh ∂yr ∂x p ∂xq ∂z j  ∂yk ∂xh ∂x p  ∂z j p p ∂M lnhr ∂N nm ∂ 2 u0 m ∂M lnhr ∂N nm ∂ 2u0 m ∂ ∂ − (aijkl )− (aijkl ) ∂z j ∂zk ∂yh ∂xr ∂x p ∂z j ∂zk ∂yr ∂xh ∂x p

−

∂ ∂z j

35

pq ∂Dlnh p ∂ 2u0 m ∂Dlnh ∂N nm ∂ 2 u0 m ∂M lmhr ∂ 2u0 m ∂ ∂ ∂ − (aijkl N nm )− (aijkl )− (aijkl ) ∂z j ∂zk ∂xh ∂x p ∂z j ∂zk ∂yh ∂x p ∂xq ∂z j ∂zk ∂xh ∂xr

−

∂T ∂P ∂T ∂P ∂T ∂ ∂ (aijkl P1k 0 ) − 1k (aijkl 0 ) − (aijkl M lnh 1k 0 ) ∂y j ∂xl ∂y j ∂xl ∂z j ∂yh ∂xl

∂M lnh ∂P1k ∂T0 ∂P0 k ∂T0 ∂ ( P0 k aijkl ) ∂T0 ∂M lnh ∂ (aijkl P1k ) ∂T0 − aijkl − aijkl − − ∂zl ∂yh ∂x j ∂zl ∂x j ∂z j ∂xl ∂zl ∂y j ∂xh −

∂Pα1 ∂T0 ∂ (aijkl M lnh ) ∂P1k ∂T0 ∂P0 k ∂ (aijkl Nα1 ) ∂T0 ∂ (aijkl ) − − ∂y j ∂yl ∂xα1 ∂z j ∂yl ∂xh ∂zl ∂y j ∂xα1

∂ 2 Pα1 ∂T0 ∂M lnh ∂ ∂Pα1 ∂T0 ∂ (aijkl P0 k ) ∂Nα1 ∂T0 ∂ h − (aijkl M ln ) − (aijkl ) − ∂z j ∂yl ∂yh ∂xα1 ∂z j ∂yl ∂yh ∂xα1 ∂z j ∂yl ∂xα1 +

∂T ∂ ∂ ∂ βij (T0 − T ) ) + ( β ij Nα1 0 ) + ( ( βijT3 ) ∂x j ∂y j ∂xα1 ∂z j

−

∂ ∂ (aijkl ∂y j ∂xl

∫ h(y, y% )µ

f 1

Y

∂

∫ h(y, y% ) ∂x µ Y

(x, y% )dY% (aijkl )

f 1

l

∂M ∂ f ∂ ∫Y h(y, y% ) ∂xn µ1 (x, y% )dY% ) − ∂z j aijkl ∂yh

∂ ∂ (aijkl M lnh − ∂z j ∂yh − aijkl

∂ (x, y% )dY% ) − ∂y j

h ln

∂µ1f (x, y% ) % ∫Y h(y, y% ) ∂x j dY

∂ ∂ f ∂ ∂ f ( ∫ h 0 (y , z, z% ) µ (x, y , z% )dZ% ) − (aijkl ∫ h 0 (y , z, z% ) µ (x, y , z% )dZ% ) Z Z ∂zl ∂x j ∂z j ∂xl

−

∂M lnh ∂ ∂ f ∂ ∂ (aijkn ∫ h(y , y% ) µ1 (x, y% )dY% ) − (aijkl Y ∂z j ∂yl ∂xh ∂y j ∂yl

−

∂ (aijkn M lnh ) ∂ ∂zl ∂y j

−

∂M lnh ∂ ∂ (aijkn ∂z j ∂yl ∂yh

∫

Y

h(y , y% )

∫

Y

∫

Y

H p (y , y% )

∂ f ∂ ∂2 (aijkn M lnh ) µ1 (x, y% )dY% − ∂xh ∂z j ∂yl ∂yh

H p (y , y% )

∫

Y

∂µ1f (x, y% ) % dY ) ∂x p

H p (y , y% )

∂µ1f (x, y% ) % dY ∂x p

∂µ1f (x, y% ) % ∂ dY ) + (aijkl (y , z ) µklf ) + A6u5 = f i (x). ∂x p ∂x j (A.1)

Further, taking (30) into (A.1) and after some derivations and combinations, the following equalities are obtained

36

1  ∂ ∂aijnh (y ) ∂aijnh (y , z ) ∂ ∂Dlnh ∂M lnh − ( a ( y , z ) ) + − − ( a ( y , z ) )  ∂z ijkl ∂zk ∂y j ∂y j ∂y j ijkl ∂zk j  2 pq ∂M lnh   ∂N nm ∂ 2u0 m p ∂ u0 m − ∂ (aijkl (y , z ) )  + N nm + ∂ ∂z j ∂yk   ∂yh ∂x p ∂xq ∂xh ∂x p ∂yh

∫

Y

H p (y , y% )

∂µ1f (x, y% ) % dY ∂x p

∂M lm  1 + ∫ h(y , y% ) ∂ µ1f (x, y% )dY%  +  aihmr (y ) − ∂ (aijkl (y , z ) ) − aihmr (y , z ) Y ∂xh ∂z j ∂zk  hr

− aihlk (y , z ) +

pq p  ∂ 2u0 m ∂ 2 N mn ∂M lmr ∂ 2u0 n ∂N mn ∂ 2 u0 n − ∂ (aijlh (y , z ) M lmr )   + + ∂zk ∂z j ∂yh ∂xr ∂x p  ∂xh ∂xr ∂yh ∂yr ∂x p ∂xq

p ∂Pα ∂T ∂N mn ∂ 2u0 n ∂P1m ∂T0 ∂P1m ∂T0 + + + ∂ ( 1) 0 ∂yr ∂xh ∂x p ∂yr ∂xh ∂yh ∂xr ∂yr ∂yh ∂xα1

2 + ∂ ∂yh ∂yr

∂µ1f (x, y% ) % ∂ ∫Y H p (y, y% ) ∂x p dY + ∂yr

∂ µ f (x, y% )dY%  1  h 

∫ h(y, y% ) ∂x Y

1  ∂M lnl ∂aijkm ∂aijkl ∂aijkl (y ) ∂ ∂Glnh  + − − + − (aijkl ) ∂y j ∂y j ∂z j ∂zk   ∂z j ∂ym ∂ ( P0 k aijkl )  ∂Pα1k ∂T0 ∂T0   ∂ 1  ∂y ∂x + P1k ∂x  +  ∂z β ij M α1 + biα1 + β iα1 − ∂z α1 l  j  j  l

(

− aiα1kl

)

1 ∂P0 k ∂Elnh  ∂T0 ∂N k ∂T0  ∂bij ∂β ij ∂P0 k ∂aijkl ∂ − (a ) ( + )+ + − ∂zl ∂z j ijkl ∂zk  ∂xα1 ∂yα1 ∂xk ∂y ∂y j ∂zl ∂y j  j

∂F h  ∂T     − ∂ (aijkl ln )  Nα1 0 +  ∂  aijkl (y , z )  ∂ ( ∫ H 0h (y , z, z% ) ∂ µ f (x, y , z% )dZ% )   ∂z j ∂zk  ∂xα1  ∂z j  ∂xh  ∂zk Z  − aijkl (y , z ) ∂ ( ∫ h 0 (y , z, z% ) ∂ µ f (x, y , z% )dZ% ) − ∂ (aijkl ∫ h 0 (y , z, z% ) ∂ µ f (x, y , z% )dZ% ) Z ∂zl Z ∂x j ∂z j ∂xl + ∫ (aijkl (y , z )( ∂ Z ∂zl

∫

Z

 h 0 (y , z, z% ) ∂ µ f (x, y , z% )dZ% − ∂ µklf ))dZ + ∂ (aijkl (y , z ) µkfl )  = 0 ∂x j ∂x j ∂x j 

(A.2) Therefore, to make sure (A.2) always hold, D ( y , z ) , M (y , z ) , Fl (y , z ) , h

Glnh (y , z ) ,

Elnh (y , z )

and

hr

H 0h (y , z ) should satisfy the following equations,

respectively. 1  ∂ ∂aijnh ∂aijnh ∂Dlnh ∂M lnh ∂M lnh ∂ ∂  ∂z (aijkl ∂z ) = ∂y − ∂y − ∂y (aijkl ∂z ) − ∂z (aijkl ∂y ), in Z , k j j j k j k  j  h z ∈ ∂Z , (A.3)  Dln = 0,

 ∂ ∂M lmhr ∂M lmr 1 ( a ) = a − a − a − ∂ ( aijlh M lmr ) , in Z ,  ∂z ijkl ∂z ihmr ihmr ihlk ∂ z ∂z j k k  j  M hr = 0, z ∈ ∂Z ,  lm

37

(A.4)

 ∂ ∂bij1 ∂β ij ∂P0 k ∂aijkl ∂Fl ( a ) = + − , in Z ,  ijkl ∂zk ∂y j ∂y j ∂zl ∂y j  ∂z j  F = 0, z ∈ ∂Z ,  l

(A.5)

1  ∂ ∂Glnh ∂M lnh ∂aijnl ∂aijnh ∂aijnh (y ) ,in Z ,  ∂z (aijkl ∂z ) = − ∂z ∂y − ∂y + ∂y k j l j j  j  h z ∈ ∂Z , Gln = 0,

(A.6)

∂ ( P0 n aijnh )  ∂ ∂Elnh ∂P h 1 ∂ − aihnl 0 n ,in Z , ( a  ∂z ijkl ∂z ) = ∂z ( β ij M ln ) + bih (y ) + β ih − ∂z j ∂zl k j  j  E h = 0, z ∈ ∂Z . (A.7)  ln   ∂  ∂ f ( ∫ H 0h (y , z, z% ) µ (x, y , z% )dZ% )    aijkl   Z ∂xh  ∂zk   ∂ ∂ f ∂ f = − ∫ (aijkl ( ∫ h 0 (y , z, z% ) µ (x, y , z% )dZ% − µ kl ))dZ Z ∂zl Z ∂x j ∂x j

∂ ∂z j

+ aijkl +

∂ ∂ f ∂ ( ∫ h 0 (y , z, z% ) µ (x, y , z% )dZ% ) − (aijkl (y , z ) µklf ) Z ∂zl ∂x j ∂x j

(A.8)

∂ ∂ f (aijkl ∫ h 0 (y , z, z% ) µ (x, y , z% )dZ% ) Z ∂z j ∂xl

Appendix B: High-order three-scale formulae for the temperature fields Referring to [50], we have the following equality ∂Tε ε (x)  1  − ∂  kijε1ε 2 (x) 1 2 = 2 B6T0 + 1 1 B5T0 + 12 B3T0 + ε1 12 B6T1 + ε12 12 B6T2  ∂xi  ∂x j  ε 2 ε1 ε 2 ε1 ε2 ε2 + 1 ( B4T0 + B5T1 + B6T3 ) + ε1 1 ( B4T1 + B5T2 + B6T4 ) + 1 ( B2T0 + B3T1 + B5T3 ) (B.1)

ε2

ε2

ε1

+ ( B1T0 + B2T1 + B3T2 + B4T3 + B5T4 + B6T5 ) + ⋅⋅⋅ = g (x) Then, considering the coefficients of ε 2−1 -term in (B.1), M k (z ) (k = 1, 2, 3) should satisfy the following equations:  ∂  ∂M k (z )  ∂kik (y , z ) = 0, in Z , −  kij (y , z ) = ∂z j  ∂zi  ∂zi   z ∈ ∂Z . (B.2)  M k (z ) = 0, Similarly to (27), Nα1 (y ) are the solutions of the following equations:  ∂  1 ∂Nα1 (y )  ∂kiα1 1 (y ) in Y , −  kij (y ) ∂y  = ∂y , j i  ∂yi    N (y ) is Y- periodic,  α1

(B.3)

where

kij1 (y ) = 1 Z

∂M j (z )   kij (y, z ) + kip (y, z ) dz ,  Z ∂z p  

∫

38

(B.4)

and kij1 (y ) (i, j = 1, 2,3) denote the mesoscale homogenized coefficients. Also, the homogenized equations defined on macrostructure Ω can be obtained as follows:  ∂ 2 ∂T0 ( x) − ∂x (kij ∂x ) = f (x), x ∈ Ω, i j   x ∈ ΓT , T0 (x) = T% (x),  v k 2 ∂T0 (x) = 0, x ∈ Γq ,  i ij ∂x j

(B.5)

where kij2 = 1 Y

∫

Y

(kij1 (y ) + kik1 (y )

∂N j (y ) ) dy , ∂yk

(B.6)

and kij2 (i, j = 1, 2,3) denote the macroscale heat conductivity tensors. At

last,

the

higher-order

local

functions

Nα1α 2 (y ) ,

Pk (y , z )

and

M kl (y , z ) (k , l = 1, 2, 3) should satisfy the following equations, respectively.

 ∂  1 ∂Nα1α 2 (y )  ∂ 2 kiα1 2 (y ) Nα1 (y ) − kα11α 2 (y )   kij (y )  = kα1α 2 − ∂ y ∂ y ∂ y  i j i   ∂Nα1 (y )  1 in Y , − kα 2 j (y ) ∂y , j   Nα α ( y ) is Y- periodic,  12

(

)

(B.7)

 ∂  ∂Pk (y , z )  ∂kik (y , z ) ∂kik1 (y ) ∂ ∂M k (z ) − + (kij (y , z ) ) − ∂z  kij (y , z ) ∂z  = ∂y ∂yi ∂yi ∂z j i  j i   ∂M k (z )  ∂ in Z , + ∂z (kij (y , z ) ∂y ), i j  z ∈ ∂Z , (B.8)  Pk (y , z ) = 0,  ∂  ∂M kl (y , z )  ∂M l (z ) = kkj (y , z ) − kkl1 (y ) − ∂z  kij (y , z )  ∂ z ∂ z i  j j    ∂ in Z , + ∂zi ( kik (y , z ) M l (z ) ) + kkl (y , z ),  z ∈ ∂Z . (B.9)  M kl (y , z ) = 0, Remark Following the ideas proposed in (Yang et al., 2015; Yang et al., 2018), the existence and uniqueness of (B.2), (B.3) and (B.7)-(B.9) can be easily established. References Abdulle, A., Nonnenmacher, A., 2011. Adaptive finite element heterogeneous multiscale method for homogenization problems. Comput. Method Appl. M. 200, 2710-2726. 39

Allaire, G., 2003. Homogenization and two-scale convergence. SIAM J. Math. Anal. 23(6), 1482-1518. Allaire, G., Habibi, Z., 2013. Second order corrector in the homogenization of a conductive-radiative heat transfer problem. Discret. Contin. Dyn.-B 18(1), 1-36. Abdulle, A., Bai, Y., 2015. Fully discrete analysis of the heterogeneous multiscale method for elliptic problems with multiple scales. IMA J. Numer. Anal. 35(1), 133-160. Almqvist, A., Essel, E.K., Fabricius, J., Wall, P., 2008. Reiterated homogenization applied in hydrodynamic lubrication. PI Mech Eng J-J Eng 222, 827-841. Allaire, G., Briane, M., 1996. Multiscale convergence and reiterated homogenization. P. Roy. Soc. Edinb. A 126, 297-342. Babuska, I., 1976. Homogenization and its applications, mathematical and computational problems. in: B. Hubbard (Ed.), Numerical Solutions of Partial Differential Equations-III, Academic Press, New York. Babuska, I., Caloz, G., Osborn, J., 1974. Special finite element methods for a class of second order elliptic problems with rough coefficients. SIAM J. Numer. Anal. 31, 945-981. Bourgat, J.F., 1977. Numerical experiments of the homogenization method for operators with periodic coefficients. Computing methods in applied sciences and engineering, I, pp. 330–356, Lecture Notes in Math., 704, Springer, Berlin. Bacigalupo, A., 2014. Second-order homogenization of periodic materials based on asymptotic approximation of the strain energy: formulation and validity limits. Meccanica 49, 1407-1425. Bhattacharjee, S., Matouš, K., 2016. A nonlinear manifold-based reduced order model for multiscale analysis of heterogeneous hyperelastic materials. J. Comput. Phys. 313, 635-653. Boutin, C., 1996. Microstructural effects in elastic composites. Int. J. Solids Struct. 33(7), 1023-1051. Bensoussan, A., Lions, J.L., Papanicolaou, G., 2011. Asymptotic Analysis for Periodic Structures. American Mathematical Society, Rhode Island. Cioranescu, D., Donato, P., 1999. An Introduction to Homogenization. Oxford University Press. Chen, Q., Zhu, H.H., Yan, Z.G., Woody, J. J., Jiang, Z.W., Wang, Y.Q., 2016. A multiphase micromechanical model for hybrid fiber reinforced concrete considering the aggregate and ITZ effects. Constr. Build. Mater. 114, 839-850. Cao, L.Q., 2005. Iterated two-scale asymptotic method and numerical algorithm for the elastic structures of composite materials. Comput. Method Appl. M. 194, 2899-2926. Dong, H., Cui, J.Z., Nie, Y.F., Yang, Z.H., Wang, Z.Q., 2018. High-order three-scale computational method for heat conduction problems of axisymmetric composite structures with multiple spatial scales. Adv. Eng. Softw. 121, 1-12. Dong, H., Zheng, X.J., Cui, J.Z., Yang, Z.Q., Yang, Z.H., 2019. High-order three-scale computational method for dynamic thermo-mechanical problems of composite structures with multiple spatial scales. Int. J. Solids Struct. 169, 95-121. Dvorak, G.J., 1992. Transformation field analysis of inelastic composite materials. P. Royal. Soc. A, Math. Phys. Sci. 437, 311–327. E, W.N., Engquist, B., 2003. The heterogenous multiscale methods. Commun. Math. Sci. 1, 87-132. Farhat, C., Harari, I., Franca, L.P., 2001. The discontinuous enrichment method. Comput. Method Appl. Mech. Eng. 190, 6455–6479. 40

Fish, J., Yuan, Z., 2007. Multiscale enrichment based on partition of unity for nonperiodic fields and nonlinear problems. Comput. Mech. 40, 249–259. Fish, J., Kuznetsov, S., 2010. Computational continua. Int. J. Numer. Meth. Eng. 84, 774–802. Fish, J., Yang, Z.Q., Yuan, Z.F., 2019. A second-order reduced asymptotic homogenization approach for nonlinear periodic heterogeneous materials. Int. J. Numer. Meth. Eng. 119, 469–489. Fish, J., Yu, Q., Shek, K., 1999. Computational damage mechanics for composite materials based on mathematical homogenization. Int. J. Numer. Meth. Eng. 45, 1657-1679. Guan, X.F., Liu, X., Jia, X., Yuan, Y., Cui, J.Z., Mang, H.A., 2015. A stochastic multiscale model for predicting mechanical properties of fiber reinforced concrete. Int. J. Solids Struct. 56-57, 280-289. Guan, X.F., Liu, X., Jia, X., Yuan, Y., Cui, J.Z., Mang, H.A., 2015. A stochastic multiscale model for predicting mechanical properties of fiber reinforced concrete. Int. J. Solids Struct. 56-57, 280-289. Geers, M.G.D., Kouznetsova, V.G., Brekelmans. W.A.M., 2010. Multi-scale computational homogenization: Trends and challenges. J. Comput. Appl. Math. 234(7), 2174-2182. Geers, M.G.D., Kouznetsova, V., Brekelmans, W.A.M., 2001. Gradient-enhanced computational homogenization for the micro–macro scale transition. J. de Physique IV 11, 145–152. Holmbom, A., Svanstedt, N., Wellander, N., 2005. Multiscale convergence and reiterated homogenization of parabolic problems. Appl Math-CZECH 50, 131-151. Hughes, T.J.R., 1995. Multiscale phenomena: Green's functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods. Comput. Method Appl. M. 127, 387-401. Hou, T.Y., Wu, X.H., 1997. A multiscale finite element method for elliptic problems in composite materials and porous media. J. Comput. Phys. 134, 169-189. Kouznetsova, V.G., Geers, M.G.D., Brekelmans, W.A.M., 2004. Multi-scale second-order computational homogenization of multi-phase materials: a nested finite element solution strategy. Comput. Method Appl. M. 193(48-51), 5525-5550. Kuznetsov, S., Fish, J., 2012. Mathematical homogenization theory for electroactive continuum. Int. J. Numer. Meth. Eng. 91,1199-1226. Li, Z.H., Ma, Q., Cui, J.Z., 2016. Second-order two-scale finite element algorithm for dynamic thermo-mechanical coupling problem in symmetric structure. J. Comput. Phys. 314, 712-748. Miehe, C., 2002. Strain-driven homogenization of inelastic microstructures and composites based on an incremental variational formulation. Int. J. Numer. Meth. Eng. 55, 1285–1322. Mosby, M., Matous, K., 2016. Computational homogenization at extreme scales. Extreme Mech. Lett. 6, 68-74. Moulinec, H., Suquet, P., 1998. A numerical method for computing the overall response of nonlinear composites with complex microstructure. Comput. Method Appl. M. 157, 69–94. Ming, X., 2002. A study on nonlinear constitutive law and FEA of rubber-like hyperelastic materials. Beihang University. Mahnken R, Dammann C. 2016a. A three-scale framework for fibre-reinforced-polymer curing Part I: Microscopic modeling and mesoscopic effective properties. Int J Solids Struct 2016; 100-101: 341-355 Mahnken R, Dammann C. 2016b. A three-scale framework for fìbre-reinforced-polymer curing 41

part II: Mesoscopic modeling and macroscopic effective properties. Int J Solids Struct 2016; 100-101: 356-375 Nascimento, E.S., Cruz, M.E., Bravo-Castillero, J., 2017. Calculation of the effective thermal conductivity of multiscale ordered arrays based on reiterated homogenization theory and analytical formulae. Int. J. Eng. Sci. 119, 205-216 Oleinik, O.A., Shamaev, A.S., Yosifian, G.A., 1992. Mathematical Problems in Elasticity and Homogenization. North-Holland, Amsterdam. Oliver, J., Caicedo, M., Huespe, A.E., Hernandez, J.A., Roubin, E., 2017. Reduced order modeling strategies for computational multiscale facture. Comput. Method Appl. M. 313(1), 560-595. Oskay, C., Fish, J., 2007. Eigendeformation-based reduced order homogenization for failure analysis of heterogeneous materials. Comput. Method Appl. M. 196, 1216-1243. Patil, R.U., Mishra, B.K., Singh, I.V., 2017. A new multiscale XFEM for the elastic properties evaluation of heterogeneous materials. Int. J. Mech. Sci. 122, 277-287. Rodríguez, E.I., Cruz, M.E., Bravo-Castillero, J., 2016. Reiterated homogenization applied to heat conduction in heterogeneous media with multiple spatial scales and perfect thermal contact between the phases. J. Braz. Soc. Mech. Sci. 38, 1333-1343. Ramírez-Torres, A., Penta, R., Rodríguez-Ramos, R., Merodio, J., Sabina, F.J., Bravo-Castilleroc, J., Guinovart-Díazc, R., Preziosi, L., Grilloa, A., 2018. Three scales asymptotic homogenization and its application to layered hierarchical hard tissues. Int. J. Solids Struct. 130-131, 190-198. Ramírez-Torres, A., Penta, R., Rodríguez-Ramos, R. Alfio Grillo, A., 2019. Effective properties of hierarchical fiber-reinforced composites via a three-scale asymptotic homogenization approach Math. Mech. Solids DOI: 10.1177/1081286519847687. Smyshlyaev, V.P., Cherednichenko, K.D., 2000. On rigorous derivation of strain gradient effects in the overall behavior of periodic heterogeneous media. J. Mech. Phys. Solids 48, 1325-1357. Temizer, I., Wriggers, P., 2011. Homogenization in finite thermoelasticity. J. Mech. Phys. Solids 59, 344-372. Terada, K., Kurumatani, M., Ushida, T., Kikuchi, N., 2010. A method of two-scale thermo-mechanical analysis for porous solids with micro-scale heat transfer. Comp. Mech. 46(2), 269-285. Yu, X.G., Cui, J.Z., 2007. The prediction on mechanical properties of 4-step braided composites via two-scale method. Compos. Sci. Technol. 67, 471-480. Yang, Z.Q., Cui, J.Z., Nie, Y.F., Ma, Q., 2012. The second-order two-scale method for heat transfer performances of periodic porous materials with interior surface radiation. CMES: Comp. Model. Eng. 88(5), 419-442. Yang, Z.Q., Cui, J.Z., Sun, Y., Ge, J.R., 2015. Multiscale computation for transient heat conduction problem with radiation boundary condition in porous materials. Finite Elem. Anal. Des. 102-103, 7-18. Yang, Z.H., Cui, J.Z., 2013. The statistical second-order two-scale analysis for dynamic thermo-mechanical performances of the composite structure with consistent random distribution of particles. Comp. Mater. Sci. 69, 359-373. Yu, Q., Fish, J., 2002. Multiscale asymptotic homogenization for multiphysics problems with multiple spatial and temporal scales: a coupled thermo-viscoelastic example problem. Int. J. Solids Struct. 39, 6429-6452. Yang, Z.H., Zhang, Y., Dong, H., Cui, J.Z., Guan, X.F., Yang, Z.Q., 2017. High-order three-scale 42

method for mechanical behavior analysis of composite structures with multiple periodic configurations. Compos. Sci.Technol. 152, 198-210. Yang, Z.Q., Sun, Y., Cui, J.Z., Yang, Z.H., Guan, T.Y., 2018. A three-scale homogenization algorithm for coupled conduction-radiation problems in porous materials with multiple configurations. Int. J. Heat Mass. Tran. 125, 1196-1211. Yang, Z.Q., Sun, Y., Cui, J.Z., Ge, J.R., 2019. A three-scale asymptotic expansion for predicting viscoelastic properties of composites with multiple configuration. Eur. J. Mech.- A/Solids 76, 235-246. Yuan, Z., Fish, J., 2009. Multiple scale eigendeformation-based reduced order homogenization, Comput. Method Appl. M. 198, 2016-2038. Yuan, Z., Fish, J., 2009. Hierarchical model reduction at multiple scales, Int. J. Numer. Meth. Eng .79, 314-339. Yang, Z.Q., Hao, Z.W., Sun, Y., Liu, Y.Z., Dong, H., 2019. Thermo-mechanical analysis of nonlinear heterogeneous materials by second-order reduced asymptotic expansion approach. Int. J. Solids Struct. 178-179, 91-107. Zhang, J.L., Liu, X., Yuan, Y., Mang, H.A., 2015. Multiscale modeling of the effect of the interfacial transition zone on the modulus of elasticity of fiber-reinforced fine concrete. Comput. Mech. 55, 37-55. Zhang, S., Oskay, C., 2016. Reduced order variational multiscale enrichment method for elasto-viscoplastic problems. Comput. Method Appl. M. 300, 199-224. Zhang, S., Oskay, C., 2017. Reduced order variational multiscale enrichment method for thermo-mechanical problems. Comput. Mech. 59(6), 887-907. Zhang, H.W., Wu, J.K., Fu, Z.D., 2010. Extended multiscale finite element method for elasto-plastic analysis of 2D periodic lattice truss materials. Comput. Mech. 45, 623-635. Zabaras, N., Ganapathysubramanian, B., 2009. A stochastic multiscale framework for modeling flow through random heterogeneous porous media. J. Comput. Phys. 228, 591-618.

43




An effective high-order three-scale reduced asymptotic homogenization is proposed.




Nonlinear heterogeneous materials with multiple configurations are considered.




A novel model reduction scheme is given for solving high-order nonlinear unit cell problems.




Thermo-mechanical coupling behavior of the heterogenous materials is analyzed.

None declared.