Scholarly article on topic 'Second-order two-scale finite element algorithm for dynamic thermo–mechanical coupling problem in symmetric structure'

Second-order two-scale finite element algorithm for dynamic thermo–mechanical coupling problem in symmetric structure Academic research paper on "Materials engineering"

CC BY
0
0
Share paper
Academic journal
Journal of Computational Physics
OECD Field of science
Keywords
{"Second-order two-scale asymptotic expansion" / "Finite-element algorithm" / "Dynamic thermo–mechanical coupling problem" / "Axisymmetric and spherical symmetric structure" / "Periodic configuration"}

Abstract of research paper on Materials engineering, author of scientific article — Zhi-Hui Li, Qiang Ma, Junzhi Cui

Abstract The new second-order two-scale (SOTS) finite element algorithm is developed for the dynamic thermo–mechanical coupling problems in axisymmetric and spherical symmetric structures made of composite materials. The axisymmetric structure considered is periodic in both radial and axial directions and homogeneous in circumferential direction. The spherical symmetric structure is only periodic in radial direction. The dynamic thermo–mechanical coupling model is presented and the equivalent compact form is derived. Then, the cell problems, effective material coefficients and the homogenized thermo–mechanical coupling problem are obtained successively by the second-order asymptotic expansion of the temperature increment and displacement. The homogenized material obtained is manifested with the anisotropic property in the circumferential direction. The explicit expressions of the homogenized coefficients in the plane axisymmetric and spherical symmetric cases are given and both the derivation of the analytical solutions of the cell functions and the quasi-static thermoelasticity problems are discussed. Based on the SOTS method, the corresponding finite-element procedure is presented and the unconditionally stable implicit algorithm is established. Some numerical examples are solved and the mutual interaction between the temperature and displacement field is studied under the condition of structural vibration. The computational results demonstrate that the second-order asymptotic analysis finite-element algorithm is feasible and effective in simulating and predicting the dynamic thermo–mechanical behaviors of the composite materials with small periodic configurations in axisymmetric and spherical symmetric structures. This may provide a vital computational tool for analyzing composite material internal temperature distribution and structural deformation induced by the dynamic thermo–mechanical coupling response under strong aerothermodynamic environment.

Academic research paper on topic "Second-order two-scale finite element algorithm for dynamic thermo–mechanical coupling problem in symmetric structure"

Accepted Manuscript

Second-order two-scale finite element algorithm for dynamic thermo-mechanical coupling problem in symmetric structure

Zhihui Li, Qiang Ma, Junzhi Cui

PII: S0021-9991(16)00186-8

DOI: http://dx.doi.org/10.1016/jjcp.2016.03.034

Reference: YJCPH 6489

Journal of

Computational

Physics

To appear in: Journal of Computational Physics

Received date: 26 June 2015 Revised date: 4 February 2016 Accepted date: 15 March 2016

Please cite this article in press as: Z. Li et al., Second-order two-scale finite element algorithm for dynamic thermo-mechanical coupling problem in symmetric structure, J. Comput. Phys. (2016), http://dx.doi.org/10.1016/jjcp.2016.03.034

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. 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.

Highlights

• Second-order two-scale (SOTS) asymptotic expansions of the dynamic thermo-mechanical coupling problem in axisymmetric and spherical symmetric structure are obtained.

• SOTS finite element algorithm is proposed and unconditional stable implicit scheme is established.

• Anisotropic material is obtained by the homogenization with enhanced strength and minor thermal expansion in circumferential direction.

• The mutual interaction and simultaneous vibration of the temperature and displacement field are simulated and discussed.

• A new computational tool is presented for analyzing composite material internal temperature distribution and structural deformation under strong aerothermodynamic environment.

Second-order two-scale finite element algorithm for dynamic thermo-mechanical coupling problem in symmetric structure

Zhihui Li12, *, Qiang Ma12, Junzhi Cui3

1. Hypervelocity Aerodynamics Institute, China Aerodynamics Research and Development Center, Mianyang

621000, China

2. National Laboratory for Computational Fluid Dynamics, BUAA, Beijing 100191, China 3. LSEC, ICMSEC, The Academy of Mathematics and Systems Science, CAS, Beijing, 100190, China

Abstract: The new second-order two-scale (SOTS) finite element algorithm is developed for the dynamic thermo-mechanical coupling problems in axisymmetric and spherical symmetric structures made of composite materials. The axisymmetric structure considered is periodic in both radial and axial directions and homogeneous in circumferential direction. The spherical symmetric structure is only periodic in radial direction. The dynamic thermo-mechanical coupling model is presented and the equivalent compact form is derived. Then, the cell problems, effective material coefficients and the homogenized thermo-mechanical coupling problem are obtained successively by the second-order asymptotic expansion of the temperature increment and displacement. The homogenized material obtained is manifested with the anisotropic property in the circumferential direction. The explicit expressions of the homogenized coefficients in the plane axisymmetric and spherical symmetric cases are given and both the derivation of the analytical solutions of the cell functions and the quasi-static thermoelasticity problems are discussed. Based on the SOTS method, the corresponding finite-element procedure is presented and the unconditionally stable implicit algorithm is established. Some numerical examples are solved and the mutual interaction between the temperature and displacement field is studied under the condition of structural vibration. The computational results demonstrate that the second-order asymptotic analysis finite-element algorithm is feasible and effective in simulating and predicting the dynamic thermo-mechanical behaviors of the composite materials with small periodic configurations in axisymmetric and spherical symmetric structures. This may provide a vital computational tool for analyzing composite material internal temperature distribution and structural deformation induced by the dynamic thermo-mechanical coupling response under strong aerothermodynamic environment.

Keywords: Second-order two-scale asymptotic expansion; finite-element algorithm; dynamic thermo-mechanical coupling problem; axisymmetric and spherical symmetric structure; periodic configuration.

1. Introduction

*Electronic mail: zhli0097@x263.net. Tel: +86-10-8233-0957, Fax: +86-10-8231-7341.

With the rapid development of aeronautic and aerospace engineering, composite materials are widely used for manufacturing aircraft wing, satellite antenna and its supporting structures, solar wing and shell, and the thermal protection system for the space vehicle etc. The composites are mostly heterogeneous with complex micro periodic structures and with the appearance of complex and extremely severe environments; composite structures usually work under the coupled thermo-mechanical circumstances of spacecraft re-entry. The practical application of composites depends on an effective method to simulate and predict the coupled thermal and mechanical performances. Generally, the equations governing the coupled heat conduction and deforming processes for the composites possess rapidly varying coefficients, and the direct numerical simulations may be realized at cost of huge computational time. While some various numerical modeling approaches [1-3] are developed to obtain the effective thermal and mechanical properties, the homogenization methods [4-9] are introduced to reduce the equations with rapidly oscillating coefficients to the equations for the materials with homogeneous properties. By the asymptotic expansion of the corresponding physical variable, temperature or displacement for example, the homogenization technique is applied to various physical problems [4,9], including the coupled heat conduction-radiation problem [10], the coupled thermoelasticity problems [11-13], the large deformation problem [14] and the elastic problem for the hollow cylinder with discontinuous coefficients [15], and then the corresponding multi-scale finite element methods have also been established [16-18]. Based on this research, Cui & Cao [19-21] developed a Second-Order Two-Scale (SOTS) method to predict physical and mechanical behaviors of composite and perforated materials by introducing the second-order correctors, while the previous homogenization researches involve only the first-order asymptotic expansions. Su & Ma [22,23] extended the SOTS method to heat conduction, linear elasticity and thermo-elastic problems for quasi-periodic structures. Wang [24] developed the SOTS method for bending behavior analysis of composite plate with 3-D periodic configuration. Yang & Zhang [25,26] studied the second-order asymptotic solutions for integrated heat transfer problems with conduction, convection and radiation in periodic composite or porous materials. Feng, Wan & Yang [27-29] obtained the multi-scale second-order solutions for the quasi-static and dynamic thermo-elastic problems. By using the statistical SOTS method, Li [30] studied the mechanical properties of materials with a random distribution of grains and Yang [31] analyzed the thermo-elastic behaviors for the dynamic thermo-mechanical coupled response of random particulate composites.

Up to now, some researches have been performed on dynamic thermo-mechanical problems of periodic composites only in Cartesian coordinates. However, it is known that the axisymmetric and spherical symmetric

structures are often used in the engineering applications such as the axisymmetric spacecraft, cylindrical pressure vessel, culvert pipe, tube of the artillery, or the high pressure sphere tank. Multilayered or periodic composite configurations are often utilized to increase the bearing capacity of the structure. To analyze and simulate the thermal and elastic behaviors of this kind of structures, we should generalize the SOTS method to the thermo-mechanical coupling problems in cylindrical and spherical coordinate system. Chatzigeorgiou [15] analyzed the homogenization for the hollow cylinder with discontinuous properties in the cylindrical coordinate system, but the first or higher order correctors of the displacement haven't been discussed. In this work, we are essentially dealing with developing the finite-element algorithm for the dynamic thermo-mechanical coupling problems. In the earlier paper[32], the SOTS asymptotic analysis has been applied to the static elastic problem in the axisymmetric and spherical symmetric structure. In this paper, the asymptotic expansion of the temperature and displacement of the dynamic thermo-mechanical system in the axisymmetric and spherical symmetric structures will be presented, and the corresponding finite-element procedures will be implemented on the basis of the two-scale asymptotic analysis. Then, the unconditionally stable implicit algorithm is developed to simulate the thermo-mechanical coupled system, in which the explicit expressions of the homogenized coefficients in the plane axisymmetric and spherical symmetric cases are given. Some numerical examples are performed to demonstrate the accuracy and effectiveness of the present proposed algorithm and to simulate the vibration effects of the temperature increment and displacement.

For this purpose, the remainder of this paper is outlined as follows. The governing equations of the dynamic thermo-mechanical coupling problems for the axisymmetric and spherical symmetric structures are presented and the equivalent compact formulations are derived for the convenience of the SOTS analysis in Section 2. The SOTS expansions of the temperature increment and displacement are formally defined and the homogenized dynamic thermo-mechanical coupling problems, the homogenized coefficients and microscopic cell problems, are deduced and discussed in Section 3. The finite-element procedure is presented in Section 4. Some numerical examples are shown in Section 5, followed by the conclusions in Section 6. Throughout this paper, convention of summation on repeated indices is used and the letters in bold represent the matrix or vector functions or variables in the formulation. By O(£k ), k e N, we denote that there exists a constant c independent of £ and | O(£k) |< c£k .

2. Governing equations of thermo-mechanical coupling problems

2.1. Axisymmetric form

Consider two typical axisymmetric structures made of composite materials depicted in Fig.1. The infinite hollow cylinder shown in Fig.1(a) is multilayered periodically with two different materials along the radial direction, while the reinforced fibers are arranged in both radial and axial directions of the cylinder in Fig.1(b), both of which are homogeneous in circumferential direction. Also, by considering that the thermo-mechanical performance of such structures occurs only in the radial and axial directions, we can simplify the modal analysis only in the axisymmetric plane A. For convenience of our asymptotic analysis, set the indicial notation (X1, x2, ff) as the cylinder coordinate instead of the conventional (r, z, ff) . Neglecting the damping effect, the dynamic thermo-mechanical coupling equations are presented by

(a) Multilayered hollow cylinder (b) 3/4 axisymmetric structure with periodicity in two directions

Fig.1. Axisymmetric composite structures.

F F dff 1 d , F dff. ^ F del , . . —_

pcF —---—(aFxi—) + = h in Ax[0, t ],

dt x1 dxi dxi dt

F d2uf 1 d , F 3o1f2 OÎ _ . , —,

pF—2---—(O)- —^+ ^ = f in Ax[0, t ], (1)

dt x1 dx1 dx2 x1

F d2u2 1 d . F \ dof _ . . —.

pF—t--—( xO) -—^ = f in A x [0,t ],

dt x dx2

in which ff (x, t) = T £(x, t) — T0 denotes the temperature increment with T £(x, t) and T0 as the absolute temperature and reference temperature, respectively. u£ (x, t) = [u£ u^] denotes the displacement vector with u£ and u^ being the radial and axial displacement, respectively. The superscript "T" denotes the transpose operator and [0, t ] is the time interval. Suppose the plane A is composed of entire cells with periodicity £, i.e.,

A = U £ (Q + z),

where Q = [0,1]2 is the normalized cell domain and I = {z = (z1, z2) e Z2 | £(Q + z) e A} is the index

set. Assume the materials are homogeneous and isotropic, and p£, c£, a£ and /£ denote density, specific

heat, heat conductivity and thermal modulus, respectively. By the periodic configuration, we have

p£ (X) = p(x / £) = p( y ), c £ (x) = c(x / £) = c(y), a£ (x) = a( x / £ ) = a( y), / (x) = /(x / £ ) = /(y),

in which p(y) , c(y) , a(y) and /(y) are 1-periodic functions, and y e Q . Here, there exists the

problem of "Two Scale" which represents the macroscopic scale x that describes the macro geometry of a

structure and the microscopic scale y , which characterizes the micro periodic cell domain inside the structure,

and they are connected by the periodicity as y = x .

Define the stress vector C£ by

o£ = [of o2e

XT 012 ]

and the strain vector ef by

f f f f ir № du 22 uf duf +du 22 lT

e = le, e e.n I = I—— —— —- —— +--—I .

L 1 2 q> 12J L -.J

dx1 dx2 x1 dx2 dx1

e,, is the volume strain and

e = e + e2 + e;.

The constitutive relation of this coupled system reads as

C£ = D£e£ - J/0£

with D£ being the constitutive matrix,

Af + 2^e Xe Xf 0

Xe + 2^ e Xe 0

Xe Xe Xe + 2jue 0

0 0 0 ue

where Xe and je are the Lamé constants of the materials. Generally, Xe, je and /3e can be expressed by the Young's modulus Ee, Possion's ratio Ve and thermal expansion ae as

u =■

(1 + Ve )(1 - 2V e y 2(1 + v e )

1 - 2v e

respectively. From the periodic configuration of the plane A , we have similarly

A£ (x) = A( x / £ ) = A( y), 1£ (x) = n( x / £ ) = n( y), E£ (x) = E (x / £ ) = E (y), v£ (x) = v( x / £ ) = v( y), a£ (x) = a( x / £ ) = a( y).

The functions A(y) , ¡(y) , E(y) , v(y) and a(y) are also 1-periodic in Q. J denotes the column

vector

J = [1 1 1 0]T.

f = [f f2]T is the volume force, and h is the heat source term.

For the convenience of the SOTS analysis, define the following fourth-order elastic tensor

C£i (x) = A£(x)81}8kl + 21 £(x)8lk8fl, i, j, k,I = 1,2,

and the second-order stress tensor

v & dXi v x V и ^

where 8 is the Kronecker delta, and

öfj = o£ , o£22 = a 2 , a£2 = ah.

Substitute Eq.(3) into Eq.(l) and then we can define the axisymmetric dynamic thermo-mechanical coupling problem in terms of 0£ (x, t) and u£ (x, t) , as follows

£ £ д0£ д , £ д0£. „ П£ д2u2 n£ du£ , .

p c£x1---(a£x1-) + T0xfl--+T0ß —L = hx1 in A X [0, t ],

dt dx dx dtdx dt

дu -_d_

dt2 dx.

pe^4£r-dr ( xCVi -dxr )(8 ЛХ)+d7 ( xiSS ßeoe)

u £ dU £

+ 81[(Л£ + 2ߣ )+ Л£ ^-ße0£ ] = fx in A X [0, "],

0£ ( x,0) = 0,u ; ( x, 0) = ( x,0) = 0 in A,

0£ (x, t) = 0 on Г X [0, "], a£ -—nt = q on Г£ X [0, "],

f (x, t) = 0 on r x [0, t ], ojn = pt on r2 x [0, t ]. In this model, dA denotes the boundary of A that is smooth and

dA = r U r2 = r1 U r2, r1 n r 2 = r n r 2 = 0.

n = [n1 n2]T is the outward unit vector normal to the corresponding boundary. The heat flux q is applied

on r2 and the surface force p = [px p2] is imposed on r2 . At the initial time, the system is static with no temperature increment, displacement or velocity.

Furthermore, if the temperature and displacement are independent of the axial variable X2, Eq.(4) reduces to the plane axisymmetric problem. In this paper, we consider the following one-dimensional problem presented in terms of temperature QE (x, t) and the radial displacement U (x, t) as

E E dOe d , edOe. ^ ? d 2u? ^ ? du? , . r , — -, pc?X—--—(xa? -—) + T0+ = hx in [r0,rjx[0, t ],

dt 3xj dxl

px ^[ x1(r+2^) ^ax )( xw

dt dx1 dx ox dx1

u ? du ?

+(a+2M?)^+a^-poe = fx in [ro,ri]x[0,"],

O? ( x,0) = 0, u? (Xi,0) = 0,-^,0) = 0 in [ro, rj,

O£(r0, t) = 0, a ? -— (ri, t) = q, u? (r0, t) = 0, of (ri, t) = -p in [0, t ], dx1

where r0 and r denote the inner and outer radius, respectively, and p is the pressure applied on the outer surface.

2.2. Spherical symmetric form

If we consider a composite hollow sphere with inner radius r0 and outer radius rl and the constituent

materials are periodically distributed in only radial direction and homogeneous in other two directions, the spherical symmetric dynamic thermo-mechanical coupling problem is similar with the plane axisymmetric case in (5). The only attention should be paid to the fact that the volume strain becomes

< = < + 2<,

and the constitutive relation now reads as

'a+2^ ? 2a dxf F ?

_ a 2(a uf F _

_ xi _

Again the thermo-mechanical coupling problem in the spherical symmetric structure is presented in terms of O ? ( X, t ) and u? ( X, t ) by

pVxf f(xy f) + r02x^df = hx2 in K, rj x [0,T],

C/i C/Ai Ul C/Ai C/'

px2 ^ [ x2(^ + 2^) |£] (2 xi^«^ ) + ^ ( x^f)

dt dx1 dx1 dx1 dx1

d« £

+ 4(£ )«£ + 2x1A£ -2x$£e£ = fX in [r0,r1]x[0, t ], dx1

$£ ( x1,0) = 0, u£ ( x1,0) = 0,-4 x1,0) = 0 in [r„, r1], dt

f(r0, t) = 0, a£-(r1, t) = q, u£ (r0, t) = 0,of(r1, t) = -p in [0, t ],

3. SOTS asymptotic analysis

3.1. SOTS expansion for space axisymmetric problem

In this section, we use the SOTS method to analyze the space axisymmetric problem described in Eq.(4). According to the procedure of the asymptotic analysis [4,9], firstly assume that Q£(x, t) and u£(x, t) can be formally expanded as follows.

$£ ( X, t ) = f ( X, t ) + f 1 ( X, y, t ) + £ 2 $2 ( X, y, t ) + O(£ 3 ), u£ (X, t) = u0 (X, t) + £U1 (X, y, t) + £2u2 (X, y, t) + O(£3),

where e0(X, t) and u0(X, t) are homogenized solutions independent of the microscopic variable y .

$1 (X, y, t) , $2 (X, y, t) , u1 (X, y, t) and u2 (X, y, t) are called as the first- and second-order correctors of the temperature increment and the displacement, respectively. Respecting

d d 1 d -->-+ £ -,

taking Eqs.(7) and (8) into the coupling equations in Eq.(4), rewriting these two equations above as the £ power series, equating the coefficients of the same power on both sides, and taking into account of the properties of 00(x,t) and u0(x,t), one obtains the first two systems of equations:

O( £-1):

-f (ax $-f (ax f) = 0,

dy, dy, dy, dx,

-( xC,

dul)( xCm du0) j «10)+± jf) = 0.

dx, dy

O(?0) :

pxf ■

--(ax1 —-)--(axl —1 )--(ax1 —1 )

dt 3y( 3y( 3y( dxt 3x( dyt

-f (axi O + T0 xftl- (+ dr) + F = hxi,

axi axi at ayi axi at

Px ^-A(xc ^uL-A(xc M)-A(xc ^uL)-A(xc

w dy dy dy dx, dx, dy dx, dx,

-A (jjOui)(^0) + A (SvxiPOi) {S^O)

+sn[a+2M) ^+aM+du0)-o=fxi.

xi m dxk

From Eq.(9), the first corrector Q1(X, y, t) and u\(X, y, t) are formally defined by

O( x, y, t ) = N_( y )

dOp( x, t ) dxa

uf ( x, y, t ) = Nam ( y ) + Mk ( y) ^^ - H ( y O ( x, t), a = 1,2,

1 dx„ x.

where Na (y) , Nakm(y), Mk (y) and Hk (y) are 1-periodic functions in Q, and satisfy the following four cell problems, respectively,

d , dNa, da . --(a--) =- in Q,

JX dy = 0,

- T- (Cjki ) in °

L^kJy = 0,

d dMk da . 3Vj dyl dyt \Mkdy = 0,

-f (Cj, fk)=f in o,

cfyj dyl àyi \oHkdy = 0.

According the formulations above, we have the following important relationships:

r c..ad-^dy = f a^mdy,

Jo « dy, y Jo dy

ir dH ( odNam,

ioj j = ¡0 P — tfy,

J =i pdM-y

10 (18)

>Q dy, Jq dy,

In fact, the weak forms for N^ (y) , Nakm(y), Mk(y) and Hk(y) can be easily obtained from Eqs.(13) (14) and (15) as follows.

f C dNakm dfr = - C dfr

jQCjk/ dy, dy dy = )qC,

JQ m~ dy, dy^ JQ" dy,

k ri _

dy = -i ÀÀ^^-dy

JQC- H H ^

JQ "k dy, dyj

where ( are arbitrary 1-periodic test functions. Making use of the symmetrical characteristic of c tl, Eq.(16) follows by taking ( = Mt in Eq.(19) and ( = Namm in Eq.(20), Eq.(17) follows by taking ( = H , in Eq.(19) and (= Nam in Eq.(21), and Eq.(18) follows by taking (= H, in Eq.(20) and (= M, in Eq.(21).

Next, we proceed to Eq.(10). By taking into account the expressions of 01(x, y, t) and u\(x,y, t) in Eq.(11), the homogenized thermo-mechanical coupling equations for 00(x,t) and u0(x,t) are obtained on the basis of an integral average over Q as

(pc)0 x $ -A. (a0 x, + T0 xfi.^«- + T0fi = hx, y ' 1 dt dxr 1 dx/ 0 1 j dx,dt dt 1

p d2«j d dul s d À 0 + d „1$

p x1~dr-T- (x1Cm^rk-)-(À,u1 )+ (x1PA)

dt dx, j dx, dx,

V «L + Ài -j- - ftL

x, j dx.

+Sn[(À + 2^f — + -Pie,] = fx,

\0 /-.0 p0 p1 p0 p1 n0

where (pc)0, p0, a0, fi, fi, fi fi j, (À + 2^)1, ÀÀ and ÀÀ are the so-called "homogenized" or "effective" coefficients defined by

(pc)0 I (pc - T0fidrL)dy,p° =JWt L pdy, a = I (a + a T^dv,

\Q\ J Q dyk \Q |J Q \Q |J Q dyk

fi = Q JQ (fi+ fi]bdyfi = Q ^Q (fi+ C- it*' fi = Q Jq d,fi = Q JQ fi+Àt )dy,

Ck, = Q I(Cmi + Cjma dNmL)dy,(À + = Q \ (À + 2M + AdM±)dy.

j \Q\JQ j dyai \Q\JQ dyk

À0 =-L r (SÀ+r-^dy,!1 =-L r À )dy

j \Q\J^ij ijU dy, ij \Q\JQ ij dyt 1

respectively, where | Q | denotes the measure of the domain Q. By Eqs.(16), (17) and (18), there are Now the homogenized dynamic thermo-mechanical coupling problem can be defined as follows

d2U0 d

( xfi"i $-t w( x/)

+ Sj[(A + + M = fx in A x [0, t ],

x1 dxk

d0( x,0) = 0, u 0( x,0) = 0,—'-(x ,0) = 0 in A, dt

00(x, t) = 0 on r1 x [0, t ], a0—n = q on T2 x [0,7],

u 0 (x,0) = 0 on r1 x [0, T],a°nj = pi on r2 x [0, t ]. and then the homogenized constitutive relation reads as

— D"e0 -fi"0o —

Cm C0 1122 4 C0 1 1112 du0 "Aû

C0 2211 C0 2222 C0 2212 dx2 322

401 4 (4 + 4 u0 3" H"

C0 1211 C0 1222 4 C0 1212 _ X1 3

du0 + du°

dx2 dxL

It can be seen that the entries of D° in the circumferential direction i.e., (A + 2p) and A® are

different from the entries in the other two directions, and the homogenized thermal modulus is also not

the same as 3° and //2. For example, we can choose the cell domain Q with good symmetry such that the homogenized coefficients have the following relationships.

C o — C o C o — C o C o — C o — C o — C o — 0 C1111 — C 2222, C1122 — C2211, C1112 — C1211 — C1222 — C2212 — 0,

Aji — À02 , 4.2 — 4 — 0, A" — 3/2 , A" — 0 •

However, there are obvious distinctions between C10111 and (A so are /I and , which

implies that homogenized material is anisotropic.

To capture the local oscillating behavior of the temperature, displacement and stress fields more precisely, we need the second-order correctors. For this purpose, turning back to Eq.(10) and defining the specific

expressions of the second-order correctors û2 (x, y,t) and Uk (x, y, t) as

x , x32a0(x,t) ^ , , 1 aa0(x,t) ,3a0(x,t) e2(x,y,t) = N (y) ' 7 + Ma(y)--^^-R(y)- 0

dxa dxa

x1 dxa

- T,s, y) - T„ (y,.1 d»10<x 0

dxa dt

2/ x lr / x, t) ,, / 1 x, t) w10( x, t)

ul(x,y,t) = N (y) . md ' + M„1t„(y)--^^-pk(y) •

1 2 dxadxa2 x1 dxa x1

- a (y)mx» - „,(y^- z„(y)^

dxa x, at

a1,a2 = 1,2.

where N,

aa Ma , R , , U, Naam, Mam, P , Vak , Wk and Z„1k are the 1-periodic cell

functions defined in Q. Taking Eqs.(25), (11) and (22) directly into Eq.(10), and the cell problems for these functions are obtained as

3 , Na,, 3 r M dNa + „ 0„ . -3-= a— (aNa) + + aSaa -a 5aa in Q,

3V' dy, dya 1 3y„ 12 ^ 2

J>a«2 dy = 0,

d dMa dNa 0

= a—— + a^1a, -a in Q,

dy, dy, dy1

iX^ = 0,

-d1 (a 1R) =^C - W^T - in Q,

dy, dy, dy'

\Rdy=a

-d_ (a ^ = $3^ + p8 in q,

dy, oy, vyk

iQSaa dy = 0

d . dU. ndM, 0 . „

-—(a—) = in Q,

dy, dy, dyk

d dN, ___ (C

aa2km ) = d (c N ) + C ^Na1km ~(Cjkl ^ ) = ^ (Cjka2 N akm) +C akl ^

ayj dy, dyj 721 2 dy,

+ Caa - C

\0Naa2km,dy = 0,

akm x ^ r 9 ^ N 1 ^ 3Mk

Ck, = ¿Ufc-(CjtaMt) + C,aki ~ + ]

% dy 3yj 1 'dy, 1 1

+ Cnu + Cama- Cla+dd- (¿Na:m)+ in Q,

dy, 1 dy, 1 dyt ¡„Mamdy = 0,

-d-Cj^) = d-(CmMk ) + + 8i[WMt+ W + 2ß) - (Я + 2^] in Q,

"V "V %

\QPkdy = 0,

-d-(CjUd-VaaL)=d-(ciJkaHk)++dd-jn,)m e,

J/akdV = a

d _ dW^ „ dHh , „ „ „о Э _ „ _ЭЯк

";Jki "

-T" (Cjkilp) = C;iki^ + 8iß-ß0 (SjXHJ -8Л(Л-^ + ßß in Q, oyj oy, oy, oyj ОУк

\wkdy=о

-^(Cj-d^) = (P-pVa in Q, oVj

JQ Z akdy = 0

Finally, the SOTS expressions of Q£ and u£ are presented by

0\ x, t ) = 0( x, t) + eN( y)

+ e\ Na y )

Эв(x,t) + M (y)_1 Э0о(x,t)

dx„ dx

дв (X t) dV(x,t)

-R( y)- To Sa ( y) ^ ) - ToU ( y)

а a ч a

i du°(x,t).

dx„ dt

X dxa -) + O(e3),

< ( x, t ) = u0( x, t ) + e( Nam ( y) ^^ + Mt ( y ) - Hk ( y)eo( x, t))

dxa, Xi

N ^ ( >&£!> + M,M ( y) -L ^

к ( y) u^ - a ( y) - w; ( y)0^ - Zak ( y)^ t)

dxa x i at

- Pk ( y)

■) + 0(£Ъ).

x dxa x

This is the general expansion of the axisymmetric dynamic thermo-mechanical coupling problem, and various simplified cases can be considered. In the practical engineering, the quasi-static thermo-elasticity problem is useful, and in this situation, both de and ue are independent of time variable t. Thus, the corresponding problem is formulated as follows.

Э . „ две. , .

--(a x-) = hxi in A,

dx, dx,

-^(xiCki d-f )-d- WW) + ^(xß8,в)

dxJ dx, dxJ dxJ

+ 8i[(W + 2ße ) ^ + r dUt-ß0e ] = fxi in A, x dx.

0e ( x) = 0° on ri, ae-n, = q on Г2

U (x) = 0 on Гх,ае,,п, = p, on Г2.

Based on Eq. (26), the corresponding SOTS expansion of 6k(x) and uek (x) are simplified as

r (X) = 0O(x) + eN (y) + e2(N** (y+m* (y)

1 dû0( x )

01 dxa aai dxadxa ° xi dxa

) + O(e3),

ue( X) = u°( x) + e( ( y)

du°( x )

aikm\s ' dx x

+Mk ( y)

u10( x)

H ( y)^o( x))

+e2( no( y)

52u:(x)

+Makm (y)

1 du0(x)

da axa xx dx

P ( y)

u0( x)

d^o( x)

x2 ■a1k ( y) dxa

Vak (y)-

W ( y)^ ) + O(e3).

3.2. SOTS expansion for plane axisymmetric problem

From the SOTS expansion deduced in the previous section, by directly deleting the terms in the axial direction x2 and y2, it is easy to obtain the SOTS expansions for the plane axisymmetric problem in Eq.(5).

Here, we directly give the results without further derivation. Since the macro domain is composed of the entire cells, we set all the cell functions to satisfy the zero boundary condition so that the SOTS approximate solution can exactly satisfy the Dirichlet boundary condition of the original problem. The SOTS expansions of

Qe (x1, t) and u1e (x1, t) for the problem (5) are defined by

x1, t ) = 60( x1, t ) + eN1( *) +e2( Nn( +y^

x1 dx1

- R( y1) ^dïA - To s ( y- T0U ( ) + O(e3),

dt dtdx1 x1 dt

— r)u0(x t) — u0(x t)

ue (x, t) = < (x, t) + e(N1 (y ) p ' + M(y ) - H(yje, (x1, t))

+e2( Nn( y1)

d2u°( x1, t )

+Mx( y)

1 du°( x1, t ) x1 dx1

- P( y1)

u1 ( x1,t)

-V ( X) ^ - W ( y)^ - Z ( y)^^ ) + 0(e3),

in which the first-order cell functions Nx, Nx, M and H satisfy the following problems

d dNy da .

-—(a—4 = — in [0,1],

dy1 dyx dyx N1(0) = N1(1) = 0,

-d [(A + 2M)dN] =d (A + 2p) in [0,1],

dy1 dyx dyx

L N1(0) = N1(1) = 0,

d „ dM dX -—[(X + 2^)—] =— m [0,1], dv1 dv1 dv1

m (0) = M (1) = o,

d . dH^ dB . -—[(X + 2ft)—] in [0,1],

dvi d>>i dvi

H (0) = H (1) = 0,

and the homogenized thermo-mechanical problem for 00 (X1, t) and u° (x1, t) is described as

(pc)0Xi ^-d— (Xi«0 f^) + ToxB ToB;^f = hxi in [10,x[0,F],

dt dx1 dx1 dtdx1 dt

p0 x1 a^ 3 [ x^+2D0 M.]OT)+dL (B

dt dx1 dx1 dx1 dx1

u 0 Thi0

+ (X + 2ft)1 + X0 ^ft = fa m [r>, rjx [0, 7],

$( x,0) = $, <(x,0) = 0,—4x,0) = 0 in [r0, rj,

$0^0, t) = 0, a °—(r1, t) = q, ->0, t) = 0,a'10(r1, t) = -p in [0, t ]. dx.

where the homogenized coefficients p0,(pc)°,a0,(X + 2ft)0,X0,(X + 2ft)1,ft0 and are defined by

p = £pdy,(pc)0 = £(pc - ^ft—) dy,

10 (a + a^N )dy, (X + 2ft)0 = £ [X + 2ft + (X + 2ft ^ ]dy,

X0 = J"1 [X + (X + 2ft) —]dy = £ (X + XdN)dy, (X + 2ft)1 = £ (X + 2ft + X—)dy,

0 dy1 0 dy1 0 dy1

, rK„ „ ,dH_ „„ „M , r^ „dH

b0=J>+/^ d=j>+a+D ^H /?=j>+Bdf- d=J0(B+^f

respectively. From Eq.(30), the homogenized constitutive relation becomes

= D°e0 -fi°d0 =

(X + 2ft)0 X0 X0 (X + 2ft)1

3x1 ft0 "

-10 ft _

_ x1 _

which also implies that the homogenized material is anisotropic.

Due to the simplicity of cell problems, the explicit expressions of the coefficients can be directly obtained

(p)0 = fW)c(t) - W-P(t)) ]dy, ) W À(t) + 2M(t)

a0 = 1/f dt,(A + 2u)° = 1/f1-1-dt, (33)

J° a(t) J° A(t) + 2^(t) 1 J

A0 = (A + 2M)0 f1 A(t) dy,(A + 2M)1 = 4fA)M) + M(t)dt + (A20, J° A(t) + 2ju(t) J° A(t) + 2^(t) (A + 2^)°

m=(A+2^)0 f—m)—d = fm(t )A+^ m) dt = r1 P"A(t)+2^(t m) dt, J0 A(t) + 2^(t) 9 J° A(t) + 2^(t) J° A(t) + 2^(t)

Based on the results above, for some simple cell structure, we can obtain the homogenized coefficients more

specifically. For example, consider that the cell domain [0,1] is composed of two different materials with the

occupation portion T and 1 - T, the coefficients of the two materials are expressed by subscript indices such

as a1 , a2 , etc., then we have

p =PT + P(1 -T),(pc)0 = plclT+p2c2{1 -T) +-T(1 -T)T°m -m)2

,(A+2m)0 =-

(A + 2m)(1 -t) + (A, + 2m2)t

(A + 2m)(a + 2m)

aj(l-t)+a2T (X + 2a)(1 -t) + (X2 + 2a2)t

Xo _ XX + 2a )(1 - T) + X(X + 2^2 )T (X + 2a)(1 -t) + (X2 + 2a2)t

(X+2a)1= (x + 2a)' + 4t(1-t)(-2I + a-a)(a-a) , (^ + 2^ )(1 - t) + (X2 + 2a2)t

g _32(X + 2a)(1 -T) + /(X + 2a)t g 0 + t(1 -t)2(a-a)(£ -p2) P (X + 2a)(1 -t) + (X + 2a)t * P (X + 2 A )(1 -t) + (X + 2a2)t

From these expressions, it is interesting to compare the coefficients (X + 2a)° with (X + 2 a)1, and also

ft' with . Following the idea of Chatzigeorgiou [15] that the Possion's ratio plays little role in the material

coefficients, then both X and A are positively correlated with the Young's modulus E due to the

expression Eq.(2), and it is concluded that (X + 2A) > (X + 2a)°. For the relationship between /0 and

, if the heat expansion a increases with bigger E, then > /0, otherwise < . Of course, for

the homogeneous material, T _ 1 or T _ 0, and there are (X + 2A)1 _ (X + 2a)' ,3* _ 3° .

For a more general case in Eq.(33), we can also have (X + 2a)1 > (X + 2a)' , the proof can be seen in Ma et al.[32].

Simultaneously the analytical solutions of N1 , N1 , M and H can be obtained by

N ( y) _-y1 + a0 fy — dt, N (y1) _-y1 + (X + 2a)' f*-1-dt,

1W * J' a(t) * K w J' X(t) + 2A(t)

x fy X0 -X(t) i „, , p 3°-3(t) 7

M,(y,)_ I -^^dt,H(y.)_ I ^-i--— dt,

11 Jo X(t) + 2A(t^ Jo X(t) + 2A(t)

respectively.

Generally, it is difficult or even impossible to solve the homogenized dynamic thermo-mechanical coupling problem in Eq.(30) analytically, but the corresponding quasi-static problem can be easily obtained explicitly. For example, consider the following model:

1 d , 0 d6>„ . --— X = 0 in [//,, /-J,

X 0X1 dx*1

1 du 0

(A + 2A)Qd?r + (A + -(A + 2m)1 ^ = m -K)~ + in Krj

30 a i„ \ a1 „.0 r\ _0 1

A>(ro) _ PAW _ a1,^) _ o,^10(/1) _ 0.

with the inner and outer temperature increment A° and A1, Ao is solved by

A _ A ln x + B, A _ (A1 - A0) / ln(r1 / ro), B _ (A0 ln r1 - A1 ln ro) / ln(r1 / ro). Then, the thermoelasticity equation becomes

(A + 2m)'

d 2u10

+ (A + 2 m) - (A + 2 m)—2 = — (C ln x1 + D)

c=(3° - m 4, D=( 4+5)mo - 5m,

The general solution for the homogenized problem is

u10 ( x1 ) = k1 x[ + k2 x1 Y + Ex1 ln x1 + Fx1,

where k1 and k2 are two undetermined constants, and

y= l(a + m e =-

l(A + 2m)0 ' (A + 2m)0 - (A + 2m)1 Then the strain solutions are obtained by

D[(A + 2M)0 - (A + 2M)1] - 2(A + 2M)0C [(A + 2m)0 - (A + 2M)1]2

with the stress component by

fe10 = k1 yxy-1 - k2 Yx-Y-1 + E ln x1 + E + F, I e0 = k1 xY-1 + k2x-Y-1 + E ln x + F,

00 = k1 [(A + 2m)0 Y+ A0 ]xY-1 + k2 [-(A + 2m)0 Y+ A°]x-r-1 + [(A + 2m)0 + A°]( E ln x1 + F ) + (A + 2m)0 E, = k1[(A + 2m)1 +AVK-1 + k2[(A + 2m)1 -AV]x-r-1 + [(A + 2m)1 + A0](E ln x1 + F ) + A0 E.

Combined with the boundary conditions, k1 and k1 are obtained as

1 V ((A + 2d)°Y+A°)

- r-r-l((X + 2d)0 + A0)(E ln r0 + F)] 1

[r0-Y-1 ((A + 2d)0 + A°)( E ln r + F) + (r0-Y-1 - r^-1)(A + 2d)0 E ')(E ln r, + F)],

[rY1((A + 2d)0 +A°)(E ln r0 + F) - (r0Y-1 - r1r-1)(A + 2d)0 E

V (-(A + 2d)°Y+A°) - r0Y-1((A + 2d)0 + A0 )(E ln r + F)]

V - rr-1r-r-1 - r-Y-1rY-1

r /0 /1 /0 '1

Based on the expressions of the homogenized coefficients in Eq.(31), the second-order correctors can also be further simplified. It immediately follows that M1 — 0 and the microscopic functions N11, R , S, U,

N11, M1, P , V , W and Z satisfy the following reduced cell problems, respectively,

d . dN11 d .

| («"7^) — ~T(«N1) in [0,1],

■j ay1 dy1 ay1 [ Nn(0) — Nn(1) — 0,

—pc-ToBdr-(pc)° 111 [0,1],

< dy1 dy1 dy1

Ir(0) — R(1) — 0,

-JL(adS-) —B-N-+B-B in [0,1],

dyl dyl dyl S (0) — S (1) — 0,

d . dU adM 0 ■ m n

-X" (a~r) — P^~ + P-P in [0,1],

dy1 dy1 dy1 U (0) — U (1) — 0,

)--f[(a + 2^^] —-f[(a + 2d)N1] in [0,1], < dy1 dy1 dy1

{Nu(0) — Nn(1) — 0,

Id . dM, d . — . dM d . dN. .

-—[(a + 2d)—1] — 3~[(a + 2d)M ] + (a + 2d)—+—(AN1)-A—1 in [0,1],

ayl dy1 ayl ayl ayl ayl

Mj(0)—Mj(1) — 0,

-—[(A+2d)—] — —(2dM) + A—+(A+2d)-(A+2d)1 in [0,1], dy1 dy1 dy1 dy1

lp(0) — P(1) — 0,

—[(A+2d)dV] —-f [(A+2d)H]+-f (BN1) in [0,1], dy1 dy1 dy1 dy1

IV (0) — V (1) — 1,

la[(A+2d)aW] —(AH)-B+B; in [0,1], J dy1 dy1 dy1 ay1

lw (0) — W (1) — 1,

—[(X + 2 a) ] _P-P in [0,1], dy1 dy1

[Z (0) _ Z (1) _ 1.

It is also possible to obtain the second-order cell functions from these problems above. As an example, considering the problem for N11. By an integration in y1, there is

dNu. N1 + C

in which C1 is a constant to be determined. Integrating the equation again over 0 to 1 and using the periodicity of N11 (y1), one obtains

C _ a0£ N1(t)dt _ a0[-2 + a0 J' J' 1

Integrating Eq.(36) in y1, we have

dsdt ].

nu( y1) _ - {; N(t )dt+C1 j; dt+c2.

Taking into account the boundary condition of N11(y1), then C2 _ 0 . Finally, from the expression of N1,

N11(y1) _1 y2 - a0 f dsdt - a0[1 - a0 f1 dsdt]iyidt.

11 1 r1 Jo Jo a(s) 2 a(s) ^ a(t)

By the same argument, the other cell problems can also be solved.

From the results above in this section, all the explicit expressions of the cell functions in the SOTS expansion of Eq.(29) can be solved successively and analytically. Moreover, although it is difficult to obtain the solution of the homogenized dynamic coupled problem described in Eq.(30), the quasi-static problem in Eq.(35) can be solved easily. Compared with the classic solutions, the static temperature field is trivial, while the thermal stress is quite different because of the anisotropic property. X is an important factor to affect the stress distribution and can be defined as the anisotropic ratio of the homogenized material. Of course for isotropic material, X _ 1.

Now that all the terms in the SOTS expansion for the quasi-static plane axisymmetric thermoelasticity problem can be expressed explicitly, it is convenient to apply this asymptotic model to predict the thermal and mechanical behaviors of the composite materials under some extreme environment of strong aerodynamic heating and force, and to design and optimize the composite structure in the engineering.

3.3. SOTS expansion for spherical symmetric problem

The SOTS expansion for the spherical symmetric problem described in Eq.(6) is similar with the one of the plane axisymmetric problem. Also, without derivation, we directly have the following SOTS expansion of

0e(xj, t) and u£ (xj, t) as

0e ( Xj, t ) = 0( xj, t ) + £ Nj( yj)

0( xj, t )

■ e 2( Nn( yj)

d 20o( xj, t ) 3xj2

- R( y1)

0( xj, t ) dt

32Uj°( Xj, t )

- ^ ( yj)^^. - ToU ( + 0(e3),

2 duj0( Xj, t )

u £ ( Xj, t) = uj ( Xj, t) + e ( Nj(yj)-

u°( Xj, t ) — .2u0( Xj, t )

-+ M ( yj)-

- H ( yj)0o( Xj, t ))

+e 2( Njj( yj) ^^X^+Mj( yj)^ ^u^ - P( yj)Zufel

Xj dxj

-K ( yj) - r ( yj)0^ - z ( yj) dX^ ) + ^

ox, X, dt

where N1, N1, M and H satisfy the same cell problems as the plane asymmetric case, respectively. The homogenized spherical symmetric thermo-mechanical coupling problem becomes:

op0x2dr-jr(x2fl0dr)+dr=hx2 in [ro,i]x[o,r],

Ui C/.A'I C/.A'I UI C/.A'I Ui

p x2 ^-dd- [ X2a + 2A)0 -id" (2 X") + ^ (x^)

dt dx2 ax1 ax1 dx2

+ 4(X + a)0«0 + 2^° ^-= fx in [ro,rjx[0, t ],

ro (x, 0) = 0, «0 (x, 0) = 0,—4x,0) = 0 in [r0, r2], dt

0o(r0,t) = 0,a" drL(r1,t) = q,«0(ro,t) = 0,O"10(r1,t) = -p in [0,7]. ox1

where the homogenized coefficient (X + A)° is defined by

,0 ^ fir,o \ ri3X(t)A(t) + 2ft2(tK (X0)2

(X + a)0 = I [(X + a) + X-]dy = I —; —i-±jL dt+—J~—,

JoLV dy, Jo X(t) + 2A(t) (X+2a)0

and the other coefficients are the same as the plane axisymmetric case in Eq.(33). We can also consider

obtaining the explicit formula for the simple structure in [0,1] as the axisymmetric case in Eq.(34). Under

this condition, there is

(X + A)0 =

:{iT (X+A)(X2+2m2) + (j -r)V(X+ A)(X+2Aj) ]2

(X+ 2Aj)(j -T) + (X + 2A2)T +t(J -r)[(V(X+A)(X + 2 a) -J (X +A2)(X2 + 2A2) )2 - (X-Xz)2]}.

The homogenized constitutive relationship becomes

= D°e0 -fi°Û0 =

(A + 2M)0 A0 dx1

2A° (A + M)0 _ u0 P _

_ x1 _

It is also concluded that the anisotropic material is obtained after homogenization.

Analogous to the problem in Eq.(35), it is also interesting to consider the following quasi-static spherical symmetric thermo-mechanical model

-d-(xy r = 0 in ft,l],

C'/vVi C/vVi

d^ + (A + 2M)0- ^ AW = e\eQ{rx) = e\0(rQ) = o,^) = 0.

(A + 2m)0^ut + (A + 2m)°— —u1— [2(A + M)°-A0]2^ = (3°-Pi)— + 3°in h,l],

2r n0 dr

The expression of A0 is given by

4 , ^(A-A) - = r -ror

G0 = —+ 5,4

Substituting A0 into the thermoelasticity equation, we have

d2y0 0

(A + 2m)0 + (A + 2M)

0 2É.u0 -[2(A + m)0 - A0]^^^ = C1 + D-1,

with C = 25(3° - P ),D = A(P" - 2P ).

Then the general solution of u10 is solved by

u0 = k1xj'1 + k2 x(2 + Ex1 + F,

v_-\^¡\+8t v =-1^/1+87 _= 2(A + M)0 -A0

/1 ~ ,Y2 ~ , /

, F = -

(A + 2M)0 D

2[(A+2m)0 -2(A+M)0 + A0] 2[-2(A+M)0 + A ]

with the strains given by

ff = kxYx x*~l + k2y2 xf2-1 + E,

fI = k1 xY-1 + k2 xY2 -1 + E + F,

and the stresses by

al = k1[(A + 2M)0 Y1 + 2A°]xp-1 + k2[(A + 2m)0 Y2 + 2A° ]xY -1 7T 2A"F

+ [(A + 2M)0 + 2A ]E +

al = k1[2(A + M)0 +A°Y1]x[1 -1 + k2[2(A + M)0 + AY2K2 -1 + [2(A + M)0 + A]E +2(A + M)0 F.

Taking the boundary conditions into account, the underdetermined constants k1 and k2 are obtained as

k = V ((A + Jw + A) [((A+2M)0 + 2A°)E (rY-1 - r? -1) + 2A0 F^lT],

V((A+2m) Y + 2A ) tor

k2 = =

V((A+2m) Y2 + 2A )

-[((A+2m)0 + 2A0)E (r1Y1-1 - r0Y1-1) + 2A0 F ^-

Where V = r0n-1rY-1 - ^-1rY-1.

For the second-order cell functions, only P is particularly posed, which satisfies the problem

d „ dP d „ — , dM „ „ 0

--[(A + 2m)—] =—[(-A + 2m)M] + 2A-+ 2(a + m)-2(A + M)0 in [0,1],

d_y1 dv1 dy1 dy1

P(0) = P(1) = 0,

while other periodic functions N11 , R , S, ¿7, N11 , M1 , F , JF and Z satisfy the same cell problems as the plane axisymmetric case.

4. The SOTS Finite-Element Algorithm

In this section, the Finite Element (FE) algorithm for SOTS analysis method of the space axisymmetric problem in Eq.(4) is given out. Firstly, we define the first-order two-scale (FOTS) approximation solutions lf1 and uf'1 by

"f = < +e(Naikm +- HA),

dxa1 x1

and SOTS approximate solutions Qe'2 and uek'2 by

A =lf,1 ^+M„ 11 - R l - TSa POr - W1

a_a1 dx dx a x dx dt aa dx dt X dt

a1 a2 1 a1 a2 1

uf,2 = uf,1+f2(N +M 1 ^ - p< - V in - r A - z

"k "k vva1a,km -s -s ^lvla,km -s 2 ''k ^a1k -s.2 /■

xa1 dxa2 1 x1 dxa1 x1 ^ x1 dt

respectively. The main computational procedures of our proposed method are presented as follows:

1. Set the material coefficients pe, ce, Xe, Ae (or Ee, Ve), ¡e (or ae ), and the reference temperature T0 . Determine the axisymmetric domain A and the cell domain Q . Make the spatial discretization of the two domains and obtain the computation meshes.

2. FE computation of the cell problems.

2. !. Solve the first-order cell functions N^ , N^i^ , Mk and Hk in Q. The weak form for Na is given by

i a —— dy = -f a-^^dy,

JQ dy i dy, ' JQ dyai *

in which ( is arbitrary ¡-periodic test functions in Q, and the weak forms for N^^, Mk and Hk are represented in Eqs.Q9), (20) and (21), respectively.

2.2. From Eq.(23), compute the homogenized coefficients (pc)0, p0 , a0, ¡¡° , ¡° , Cj^ ,

(X+ 2a)1 and X .

2.3. Solve the second-order ceU functions Nala2, Maj , R , Saía2, U , Na^km , Ma,km , Pk,

Vaik , Wk and Zaik in Q by the following weak forms, respectively,

f a^Na(y = f (a^ + aS - aS^ )( - f aN^dy,

Q dy . dy . JQ dya2 aa " JQ a 5y„2 ''

f dMa d( f dNa 0

)Qa~dyt didy = Jq^ + S -a Sa)(dy, fadR(y = f (pc-T05drL-(Pc)°)(dy,

jQ dyt dy, jQ dyk

f a^(y = f ¡dNaka +5S -¡a )(dy,

Jq dy, dy, JQ ^ dyk P aa Paa V y

f adU(y = f (¡¡M +5-5)(dy, Jq dy, dy, h"1 dyk P P W

f C dNggkm ^ ^^ ^ C dN0,hn + c _ ^ _ - c n ^

)q ijk dy dy. a dy '"2"°' ^'ama'WW JQSk«2ivakm dy

iQCk ^ ^dy = i/'" (C*„ M - l/'mjMk ^dy

+£ (Cflk, d-Nkm + Cima - Cla )frdy - £ M^ " M dy

- -I +

\pdrdy = I[^^dr^ + (A+ 2M) -a + 2ft)1]^1dy - joCjk^Mk Mdy - ^dy.

0 dy dy 3yk Jo dyj dyk

J0C- 3jt = Jo «0- it+*■*-£)** - ioCjkaHk f^ - Jo N f^

JC-^ ft=Jo (Ci-H+M-k)*# - ¡0 <* IH:- io^Hi ft

f C^^aLMdy = f (p-ptdy. •o j dy <V Jo 1

It is observed that all the cell problems of the temperature field are elliptic and can be easily solved by the standard finite-element method with the same global stiff matrix in solving the discretized algebraic equations and so are the cell problems of the displacement field. Therefore. it will not take too much extra time for computing the second-order cell functions.

3. FE computation of the homogenized problem

From Eq.(24), the weak form for 00 and u° is formulated by

f (-P£LxM.e dx + f ^x e^dx + f x^y ^eodx = f hx,0odx + f qXieods,

T0 1 dt 0 to 1 dx' dx' ^ ; dt 0 TO 10 Jr2 T0 1 0

£ px1 —U-Hidx + £ x1 (e0 )T D0e0dx - £ x1 (e° )T fi°60dx = £ x1 pfi^s.

where e0 and u.. i = 1.2 are arbitrary test functions in A that satisfy e0 = 0 on r1 and ui = 0 on r1.

e0 = [P 0 P 0 P 0 P 0 ]T = [ du1 du2 u1 du2 + du1 ]T

dx1 dx2 x1 dx1 dx2

The homogenized problem is strongly coupled by parabolic and hyperbolic equations. so the discretization in the spatial domain is carried out firstly. and then the time stepping scheme is used in the temporal domain

[0, t ]. According to this procedure, let the homogenized displacement u° and temperature lo be presented in terms of node values de and ? by

u _[N^)]^de(t),0o _[N'(x)\n?(t),

respectively in each element, where Nv denotes the number of nodes per element. The [N] and [N'] denote the shape functions of u° and lo, respectively. Moreover, the nodal strain and thermal gradient are

expressed by

*0 = [ BW d e, = [ b \Näe,

with [B] and [B'] denoting the gradient matrix for u° and do, respectively. Based on Eq.(39), the space discretized system for the dynamic thermo-mechanical problems is given by

x M?e (t)+K?e (t)+(Ce )T de (t) _ Ge (t),

^ Meude (t) + Kide (t) - Ce?e (t) _ Fe (t).

where Ne is the number of elements, the superscripts " •" and " ••" denotes the first- and second-order time derivatives, respectively, and

Ml _ jA,PV']T[N']x1dAe,Kl _¡^[B']T[B^,

Ml _ ¡A, p°[N]T[N]xxdAe, Kl _ J [B]T^0[B]x1dAe,Ce _ [B]T fi0[N']x1dAe,

Ge _ L [N']T ~X1dAe + f [N']T qx1dse,Fe _ j [N]T fxxdAe + J_ [N]Tpxxdse.

JA T' 2 ^ T o A 2^

By finite element assembling process, the global space discretized system is given by

|M?(t) + K?(t) + CTd(t) _ G(t), IMud(t) + Kud(t) - C?(t) _ F(t).

To deal with the derivative terms d(t), d(t) and ?(t), divide the temporal domain [0, t ] into Nt equidistance sub-internals,

o _ to < t1 < ••• < tNt _ t,

with the time step At _ 7/Nt, t _ nAt, n _ 0,1,2,...,Nt. The unknowns at the time tn are denoted by,

dn _d(tn),dn _d(tn),dn _d(tn),? _?(tn),

For the weak form of the heat problem, we use the " 8 "-scheme, that is

?n+1 -

+ 8K?n+1 + (1 -8) K?n + C1

dn+ - dn

_8Gn+1 + (1 -8) G"

At " At

where 0 < 8 < 1, and for the dynamic thermo elastic problem, we adopt the Newmark scheme

(-X M + K )dn+1 - C?n+1 _ Fn+1 + M (—L- dn +—dn + (— - 1)dn), VAt2 u u * u tfAt nAt 2n

where the acceleration dn+1 and velocity dn+1 are computed by

dn+1 (dn+1 - dn) —— dn - (— - 1) 'dn, dn+1 _ dn + [(1 -a>)dn + codn+1 ]At,

nAt nAt 2n

and the parameters 0 <C< 1, 0 1/2. Usually, when C> 0.5 and n^ 0.25(0.5 + c)2 , the scheme (41) is unconditionally stable. Rewrite Eqs.(40) and (41) together in matrix form as

Ml + S/AtKl C 1

M + K..

Ml- (1 -8)AtKl C 1

8AtGn+1 + (1 - 8) AtGn

Fn+1 + M (— dn + ^^-1)rfn) u nAt 2n

At the initial step, d0 is computed by

Mud0 + Kud0 - C? _ F0.

Therefore, the computing procedure of the homogenized problem is programmed as follows:

3.1. Verify t , Nt and the time step At _ t / Nt. Take n _ 0 as the initial time step, and set the parameters y, C and n

3.2. Assemble the matrices Mu, Ku, C , Ml and Kl, and form the global matrix in the computational scheme

Ml + 8AtKl -C

2 Mu + Ku

Ml- (1 - 8)AtKl C ' 1

Assemble the initial global load vector F0, determine the variables d0, d0 and ?, and compute the initial velocity d0 by

H0 _ Mu'(F0 - Kud0 + C?°). 3.3. Assemble the vector Fn+1 and Gn+1 and form the global load vector

SistG n+1 + (1 -S) AtGn

Fn+1 + Mu (— dn + (— - \)dn ) TjAt 2n

Apply the corresponding boundary conditions for the temperature increment and displacement and obtain dn and £n+l by Eq.(43).

3.4. Compute the velocity iln+1 and acceleration iln+1 with the expressions in Eq.(42).

3.5. Set n = n +1, if n > Nt, finish computing, otherwise return to Step 3.3.

4. Assembly of the FOTS and SOTS solutions of temperature, displacement, thermal gradient, strain and stress.

4.1 In each time step, assemble FOTS and SOTS solutions of de and uek by Eqs. (37) and (38) dd0

respectively, where the term - is obtained by the backward difference. From the Newmark scheme, the

du0 dv av

velocity —— and acceleration -is directly obtained, while -— is computed by the space

dt dt2 dxadt

derivative of the velocity

dp d@s,1

4.2. Compute the homogenized, FOTS and SOTS solutions of the thermal gradient, » , and

dxt dxt

__, where the partial differential operation of Eq.(8) is used to compute the derivatives of the first- and

second-order correctors.

4.3 Compute the homogenized, FOTS and SOTS solutions of strain by

0 r 0 e = [e 0 0 0 e2 eQ e 12 [ du° du° dx1 dx2 x1 du° + dx2 du20]T dx1

e II «,1 «,1 e2 ee ee,l]T = [ du,« 12 J L -, dx1 du1 dx2 «,1 u1 x1 du,1,1 —— + dx2 du ]T dx1

re II e,2 e,2 el22]T = [ df dx1 du1 «,2 u1 du,2 , duf]7

e2 ee dx2 X1 dx2 dx1

and the stress approximations by

_0 *a0 0 O0 n —£,1 T\£ £,1 T Q£ n£,1 —£,2 r\£ £,2 T/3£n£,2

7 = D e -pe0,O = De - J pe ,7 = De - Jpe . For the quasi-static thermo-mechanical problem (27), by deleting all the time derivative terms, the temperature increment 4 and displacement d are obtained by the following two steps

K£ = G, Kud = C£ + F .

5. Numerical Examples and Discussions

5.1 Space axisymmetric quasi-static thermo-mechanical problem

The quasi-static thermo-mechanical problem is solved by the computable modeling on the basis of Eq.(27) as the first example. Consider a composite infinite cylinder with inner radius 1 and outer radius 2, the microscopic cell domain is shown in Fig.2(a), where there are periodicities in both radial and axial directions. The configuration of unit square [1,2] X [0,1] in the axisymmetric plane is depicted in Fig.2(b) with £ = 1/8. On the outer surface r = 2, e£ = 2K , and on the inner surface r = 1, set e£ = 0K . There is a constant heat source h = 500 J /m3 throughout the materials, and no pressures are imposed on the surfaces.

1/4 3/4

(a) Cell domain Q

□ □ □ □ □ □ z z

□ □ □ z z z z z

□ □ □ z z z z z

□ □ □ z z z z z

□ □ □ z

□ □ □ z

□ □ □ z

□ □ □ □ □ z z z

(b) Macro periodic domain in unit square

(c) Periodic domain for FE computation. Fig. 2. Cell domain and the macroscopic periodic composite domain (£ = 1/8)

Since it is difficult to find the exact solutions of above problem, we have to take de and u£ to be their FE solutions in very fine mesh for comparison. The information of computational meshes with different £ are listed in Tab.1, and it is clear that the mesh partition number of SOTS approximate solutions is much less than that of fine original solutions and independent of £ . That means that the SOTS method can greatly save the memory and computational time, which is very important in engineering computation.

Table 1. The information of the computational meshes.

Original fine mesh

£ = 1/8 £ = 1/16 £ = 1/32

Cell mesh Homogenized mesh

Number of Nodes 9393 18753 37473 1203 585

Number of elements 18208 36416 72832 2276 1204

Mesh size 6.7451E-3 3.3726E-3 1.6863E-3 5.3960E-2 2.2100E-2

The material coefficients are listed in Tab.2 and by Eq.(23) the homogenized heat conductivity is computed as a0 = 70.2916W/mK, and the computed homogenized elastic constitutive matrix is

1.1216 0.5463 0.5500 0

0.5463 1.1216 0.5500 0

0.5500 0.5500 1.1517 0

0 0 0 0.2849

x10Pa,

thus, the homogenized Possion's ratio v0 and Young's modulus E0 in radial and axial directions are given by

C0 + C0

1111 1122

= 0.3275,E0 = 2C10212(1 + vu) - 7.5641x1010Pa,

and the Possion's ratio v* and Young's modulus E* in circumferential direction are

= 0.3232,El = [(X + 2ft)1 -X101](1 + V?) - 7.9617x1010Pa,

* (0 + 2ft)1 + 0

respectively. It is observed that in the circumferential direction, the Possion's ratio is smaller while the Young's modulus is bigger than the ones in the other two directions, showing that the homogenized anisotropic material has the biggest strength in the circumferential direction. The computed homogenized thermal modulus vector is

0 = [2.3836 2.3836 2.3682 0f x106Pa, which implies that the circumferential thermal modulus is smaller than that in other two directions. Related to

the value of E° and v°, it is also indicated that the homogenized circumferential thermal expansion is also

much smaller.

Table 2. Coefficients of the constituent materials

Material 1 ( k1 ) Material 2( k2 )

Young's modulus(GPa) 66.2

Possion's ratio 0.321

Heat conductivity(W/mK) 118.1

Thermal modulus(GPa/K) 2.491E-3

117.0 0.333 2.036 1.905E-3

By the symmetry of the cylinder, we only choose the domain with one periodicity in the axial direction for our computation, that is A = [1,2] X[0, £] shown in Fig.2(c). The axial displacement on z = 0 and £ is fixed to satisfy the symmetric condition.

Fig.3 shows the different approximate solutions and the present FE fine solution of the temperature increment. Due to the discontinuity of the material coefficients, there are humps in the embedded square in the temperature field, finely revealed in Fig.3(a). The homogenized temperature increment is shown in Fig.3(b), which is the exact solution of the following one-dimensional problem:

d —p -—(a0 x = hxi, Po(1) = 0, p (2) = 2,

UAi tiAi

where Q is solved by

„ , . h ,2^ 8a + 3h,

Q0(x1) = -—(x1 -1) + 0, o lnx1'

4a 4a ln 2

and the maximum of the temperature is calculated by

Q0max =Q(x1) « 2.2564, xl =

8a0 + 3h

• 1.7249.

2h ln2

It is also seen from Figs.3 (c) and (d) that it is difficult for the FOTS solution to describe the exact distribution, while there is hardly any difference between the FE fine solution in Fig.3(a) and the present SOTS solution in Fig.3(d).

1,4 1.6

(a) FE fine solution QE (K)

1.4 1,6

(b) Homogenized solution (K)

0.15 0.1 0.05 0

1 1.2 1.4 1.6 1.8 2

(c)FOTS solution 0e,l(K)

1.2 1.4 1.6 1.8 2

(d)SOTS solution 0e'2(K) Fig. 3. Comparison of the temperature increments with £ = 1/8.

> »01 |j4| n n y

(a) FE fine solution dd£ / 3x1 (K/m)

0.15 0.1

'0.05 0

(b) Homogenized solution d Q0 / 3x1 (K/m)

(c) FOTS solution dQ1,1 / dx (K/m)

1 1.2 1.4 1.6 1.8 2

(d) SOTS solution d0£,2 / dx1 (K/m) Fig. 4. Comparison of the thermal gradients in radial direction with £ = 1/8.

The approximate solutions of thermal gradient in the radial direction is shown in Fig.4 and it is more clear that the FOTS approximation dO£,:1 / 3x1 in Fig.4(c) is not sufficient to give a good result. The SOTS

approximation dO£'2 / 3x1 in Fig.4(d) is more acceptable in very good agreement with the present FE fine

solution of £ = 1/8 in Fig.4(a). After adding the second-order correctors, the thermal gradient can be simulated accurately and effectively, and sequentially, we can obtain the thermal flux with oscillating behavior between the interface of the two materials correctly. Although extra work is done for computing the second-order correctors, it is worthwhile as the real temperature field can be reproduced precisely.

By the feature of the temperature field, the corresponding thermo-elastic problem can also reduce to the plane axisymmetric problem with the one-dimensional homogenized equation as

0 d u1 0 1 du1 (l+i,i\1 u1 = (R0 R0\Q0 + R0 dQ0

C1111~TT + C1111——_— (A + 2M) ~r = (R11 -PV) — + R11-

dxX1 dx1 X1 X1 dx1

Since it is well-known that if the material is isotropic, C0111 = (Â + 2^)1 and RR =R°, the equation above reduces to the classic plane axisymmetric thermal-stress equation. Fig.5 shows the approximate solutions of the

radial displacement and it is observed that both the FOTS and SOTS asymptotic solutions can effectively give the detailed information of the displacement field.

0.15 0.1 0.05 0

0.15 0.1 "0.05 0

0.15 0.1 0.05 0

2.5E-05 3.1E-05 3.7E-05 4.3E-05 4.9E-05

1 1.2 1.4 1.6 1.8 2

(a) FE fine solution Uj (m)

(b) Homogenized solution Uj (m)

1 1.2 1.4 1.6 1.8 2

(c) FOTS solution Uj (m)

1 1.2 1.4 1.6 1.8 2

(c) SOTS solution U{ (m) Fig. 5. Comparison of the radial displacements with £ = 1/8.

To compare the solutions quantitatively, define the following error functions for 0£ and u£ :

0 n£ n 1 n£ n£,1 2 /}£ /}£,2 0 £ 0 1 ££,12 £ £,2

ee = 0 -0o,ee = 0 -0 , ,ee =0 -0 , ,eu = u -u ,eu = u -u , ,eu = u -u , , and the relative errors of the temperature increment and the displacement in norms are listed

with different £ in Tabs. 3 and 4, where the semi-norm (energy) is used in computing the H1 norm. Although the present homogenized, FOTS and SOTS solutions are in accordance with the FE fine solutions, by expansion to the second order, the errors are improved significantly and the local oscillation within the micro cell of the temperature and the displacement are captured more accurately. On the other hand, as £ is

decreasing, all the errors in L norm also decrease while the errors of the homogenized solution in H1

norm keep almost the same. This is because the homogenized solution is incapable to describe the microscopic information so that the local fluctuation makes the errors between the homogenized and the original solution exist all the time. By adding the correctors, the errors are decreasing as £ is becoming smaller.

Table 3. Computed relative errors of the temperature increments in norms with different £ .

£ \\eQ\\L2 \\Q£\\I} \\4\\L> \\Q£\\L2 \\eQ\\l2 \\Q£\\L2 \\eQ\H \\Q£\\H1 \\eQ\H \\Q£ \\H1 \\eQ\H \\Q£ \\H1

1/8 1.7680E-2 1.0739E-2 2.4449E-3 5.6213E-1 3.3741E-1 5.89491E-2

1/16 3.7065E-3 1.7831E-3 1.0256E-3 4.6979E-1 1.2783E-1 5.5565E-2

1/32 3.3757E-3 1.1174E-3 7.4971E-4. 4.8604E-1 1.0688E-1 5.0687E-2

Table 4. Computed relative errors of displacements in norms with different £ .

£ Kl* !!«% \\ eh \\L2 \\"£\\l2 \\ e2u \\L* \\"£\\l2 \\eU \\H! \\ u£ \\Hl \\eUu \\h1 \\u£ \\H1 \\eU\\H1 \\u£ \\H1

1/8 6.1039E-3 3.9232E-3 2.8263E-3 3.5162E-1 7.8192E-2 5.6394E-2

1/16 3.5822E-3 2.6748E-2 2.7378E-3 3.5108E-1 7.1161E-2 5.0634E-2

1/32 1.4415E-3 7.9828E-4 7.8631E-4 3.4858E-1 7.0251E-2 5.0008E-2

Figs. 6 and 7 show the approximate solutions of the radial and hoop stress, respectively, from which the discontinuous stress is reproduced by adding the first- and second-order correctors. From Eq.(44), it can be known that the displacement and stress are dependent on both the temperature increment and thermal gradient. Due to the inner heat source, generally, the stress inside the embedded square is bigger than that around the matrix, with the maximum lying in one embedded square, and around this square the maximum temperature increment is located.

0.15 |-0.1

0.05 0

(a) FE fine solution 0* (MPa)

(b) Homogenized solution 0 (MPa)

(c) FOTS solution 0t,j (MPa)

1.2 1.4 1.6

(d) SOTS solution 0£'2 (MPa) Fig.6. Comparison of radial stresses.

(a) FE fine solution 0£ (MPa)

0.15 0.1 0.05 0

i ■ !

1 1.2 1.4 1.6 1.8 2

(b) Homogenized solution (MPa)

■ tu ■ t I

(c) FOTS solution 0e/ (MPa)

(d) SOTS solution (MPa) Fig.7. Comparison of hoop stresses

Similarly, define the following error functions for the stress component

eO _£ _0 el _£ _£,1 e2 _£ _£,2

Oj =0j -0j,ài =0j -0j ,Ôj =0j -0j ,

¿2 = 02 = 02 -02 ¿2 = 02 -02 , ¿0 = 00 - 00,d0=00-00 ,d0 =00 -00 ,

and the computed relative errors of the stresses in L norm are listed in Tab.5, which demonstrate that there is

higher accuracy in the SOTS approximations for the stress computation than in the FOTS and the homogenized

solutions.

Table 5. Computed relative errors of the stresses in L norm with different £ .

IIO°|L IIÔIL II02|I

II 0j£ II

II 0j£ 11 L

II 0j£ 11 L

II Ô IIL2 II 0£ I I L2

II ÔII

2 I II2

II 0£ 11L 2

II 02% II 0£ 11 L 2

IIÔ0IIL2 II Ol II

0 IIL2

II 0£ 11L 2 II 0£ IIL

II02^ II 0EO 11 L 2

1/8 8.2514E-2 2.5832E-2 1.6082E-2 1.1601E-1 3.8705E-2 3.3252E-2 8.7738E-2 1.6715E-2 1.2329E-2 1/16 8.5760E-2 1.5187E-2 1.0988E-2 1.1427E-1 3.1495E-2 3.0004E-2 8.4709E-2 1.1656E-2 1.0217E-2 1/32 8.6415E-2 1.2012E-2 1.0230E-2 1.1271E-1 3.0392E-2 3.0000E-2 8.4371E-2 1.0101E-2 9.7311E-3

From Figs.3-7 and Tabs. 1,3 and 4 in this section, it can be concluded that the SOTS approximation solution is valid to simulate the coupled thermal and mechanical behavior of the composite material in the space axisymmetric structure. The difference between the FOTS and SOTS approximate solutions may be not very significant for the displacement field, but for the temperature field, the second-order correctors are indispensable. Based on this consideration, the present SOTS finite element algorithm is always implemented

in the actual engineering computation.

5.2. Axisymmetric dynamic thermo-mechanical coupling problem

(a) Symmetric layered cell domain (b) The cross section domain of the layered

material (£ = 1/4)

Fig.8. Multilayered infinite composite cylinder and the cell domain

The plane axisymmetric dynamic thermo-mechanical coupling problem is studied in this section. Consider the infinite multilayered cylinder with inner radius 1 and outer radius 2, but this time the cell domain is depicted in Fig.8(a) and the cross section of the cylinder is shown in Fig.8(b), which means that the materials are arranged periodically only in the radial direction. The non-dimensional computation is carried out and the material coefficients are listed in Tab.6. The cylinder is static initially at t = 0, for t > 0, there is a constant heat source h = 10 applied in the cylinder. The temperature increment is set as zero on both the inner and outer surface. The inner surface is fixed while the outer surface is set free. Based on the scheme in Eq.(43), we perform the finite element computation for the homogenization solution on [0,1] by 64 elements and 65 nodes, with the mesh size 0.015625. Set the time step At = 0.005 , the parameters of the time integration S = 0.5, C = 0.5 , n = 0.25, and 7 = 10. All the cell functions are obtained analytically.

Table 6. Non-dimensional material coefficients.

Material 1( k1) Material 2( k2) Young's modulus (E) 10 1

Possion's ratio (V) 0.25 0.3

Heat conductivity (a) 10 1

Thermal modulus ( ß ) 0.01

Mass density ( p) 10

Specific heat ( c ) 4

0.05 2 3

From Tab.6, the homogenized coefficients are easily obtained by the explicit expressions in Eq.(34) with T = 0.5, therefore,

a0 -1.1812,(A + 2^)0 -1.2746, A0 - 0.4856,(A + 2a)1 - 5.7930, / - 0.4788,/?° - 0.3586, a0 = 6.0,(pc)0 - 23.3156.

The homogenized Possion's ratio V° and the Young's modulus E0 in the radial direction are given by

V° = ■

(A + 2a)0 +A0

0.2759, E0 = [(A + 2a) - A0](1 + V0) « 1.0067,

while and E° in the circumferential direction are

0.0773,El = [(A + 2a)1 -A0](1 + V) « 5.7197,

(A + 2a)1 +A0

respectively. It is seen that there is a much smaller Possion's ratio and bigger Young's modulus in the circumferential direction, which is consistent with the case in Section 5.1. Set the reference temperature T0 = 50, and we can define the following thermo-mechanical coupling ratios in the radial and circumferential directions, respectively, as

(ß°)2 T i

(ß0)2 T

c (A + 2AyW c (1 + 2a)W

Since (A+ 2a)1 > (A + 2a)° and fff < J30, in the sequential, only r0 is considered. From Carter [33], the

thermal and elastic process has become significant for the material with rC° = °.5 , and rC° = ° corresponds

to the uncoupled thermo-elastic behavior, and in this case, rC° ~ °.3856 from Eq.(45) and 7°.

Figs. 9 and 1° show the evolutions of the approximate solutions of the temperature increment and radial displacement at various time with £ = 1/8 in different line types. By a fine mesh with elements 1024 and mesh size 9.765625E-4, the FE solutions are also obtained for comparison and represented by circle symbols. Due to the continuous heat source from the beginning to the time t = 1, the temperature rises throughout, forms the parabolic curve in the radial direction with the maximum lying around r = 1.5 and causes a

transient acceleration for the displacement field, which leads to the two-way movement along the radial direction from the center r = 1.5 to the inner and outer surface (in Fig.10(a)), respectively. As the time goes on, the temperature keeps on rising (in Figs.9(b)-(d)). Since the outer surface is in unconstraint, the cylinder is expanding and all of the radial displacements become positive (in Figs.10(b)-(d)). When the displacements around the inner surface are also going to become positive, the reflection of the elastic wave occurs due to the fixed inner surface and the inertia effect at about t = 7.0, which can be clearly observed in Figs.10(e)-(f). The fluctuation of the displacement field simultaneously causes the vibration of the temperature field. In Fig.9(g), it can be seen that at about t = 9.0, the maximum of the temperature falls a little. Then, both of the temperature and radial displacement keep on rising, and the displacements in Figs.10(f)-(h) are becoming positive again throughout. It can be predicted that at t > 10.0 there are continuous oscillation and mutual influences in both of the two physical fields.

It is also seen that the macroscopic behaviors are well captured at each time step by the homogenized solutions. Based on these results, the FOTS and SOTS solutions give good approximations. The SOTS solutions almost coincide with the FE solutions in fine mesh, while there are minor errors by the FOTS solutions, which can be observed in the detailed graphs of the locality in each of the figures. All the results demonstrate the reliability and effectiveness of the SOTS method to predict the dynamic thermo-mechanical behaviors of periodic composite materials in the plane axisymmetric structure.

(a) t = 1.0 (b) t = 2.0

1.4 1.6

1.4 1.6

(c) t = 3.0

(d) t = 5.0

1.4 1.6

0.5 j-

f/ 06 . \ "

0.4 - // 0.59

CD 0.58

0.3 - j 0.57 \

. 3 0.56 \

0.2 iT 1.66 1.68 1.7 o 1.72 e£

0.1 -u 00 e£,i \\

01 , 1 , , 1 , , e£, 2 ..... ■ . . . «

1.2 1.4 r 1.6 .8 2

(e) t = 7.0

(f) t = 8.0

1.2 1.4 1.6 1.8 2

(g) t = 9.0 (h) t = 10.0

Fig. 9. Time history of temperature increments (£ = 1 / 8) .

1.2 1.4

1 1.2 1.4

(a) t = 1.0

(b) t = 2.0

1 1.2 1.4 1.6

1 1.2 1.4 1.6

(c) t = 3.0

(d) t = 5.0

(e) t = 7.0

(f) t = 8.0

-0.005

1.4 1.6

- 0.1315

0.12 0.131 0.1305 V -----^

0.1 0.13 0.1295

0.08 - 1.90 1.92 1.94 1.96 1.98

0.06 -

0.04 - o u; ----" -------u; 1

0.02 - u; 2

0 : ....... 1 i .......

1.2 1.4 r 1.6 1.8 2

(g) t = 9.0 (h) t = 10.0

Fig.10. Time history of radial displacements (£ = 1/8)

Fig. 11 shows the approximate solutions of the temperature at t = 5.0 with £ = 1/2k , k = 2,3,4,5.

The errors of the FOTS approximations in Fig.11(a) are observed more clearly in £ = 1/4, while the SOTS solutions are still almost the same as the FE solutions in fine mesh. For the approximate solutions of the radial displacement shown in Fig.12, both the FOTS and SOTS solutions can not give the correct values on r = 2 since they do not rigorously satisfy the Neumann boundary condition, but as £ becomes smaller, the errors are decreasing. On the other hand, there are no errors on r = 1, since all of the asymptotic solutions satisfy the Dirichlet boundary condition by setting the zero boundary condition for all the cell functions. The detailed graph of the locality is also illustrated in the figures, from which it is seen that sometimes the FOTS approximate solutions may be closer to the original solution, but their derivates are not approximated correctly and this brings in the errors in computing the heat flux and the stresses.

1.2 1.4 1.6 1

(a) £ = 1/4 (b) £ = 1/8

(c) £ = 1/16 (d) £ = 1/32

Fig.11. Comparisons of temperature increments at t = 5.0 with different £.

(a) £ = 1/4

(b) £

1.4 1.6

(c) £ = 1/16 (d) £ = 1/32

Fig.12. Comparisons of displacements at t = 5.0 with different £.

The computed relative errors of the temperature increment and radial displacement with £ = 1/8 and

1/16 in L and H1 norms at different time are illustrated in Fig. 13 in different line types. For the homogenized solutions, the errors are bigger than the other two approximate solutions, and as £ becomes

smaller, only the errors in L norm become smaller (see Figs.l3(a)-(b) and Figs.l3(e)-(f)) while there is not much improvement of the errors in H1 norm. By comparing Fig.l3(c) with 13(d) and 13(g) with 13(h), the local varying of the thermal gradient and the stress of the coupled system can not be simulated only by the homogenized solutions. The errors reduce significantly both in norms by adding the first- and

second-order correctors and as the periodicity increases, the difference between these two approximations are becoming smaller. Another notation is that the evolution of the errors with time is becoming more stable with decreasing £.

time(t)

£ 0.035

á 0.025 .2

® 0.02 0.015 0.01 0.005 0

-0.005

1- 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 -

: - 80 - L2 Error H

r ------------- 8e' 1-L2 Error i

~ -----8C'2-L2 Error i

time(t)

3 0.06

(a) L error of Qe (£ = 1/8)

(b) L error of Qe (£ = 1/16)

<3 0.4 u]

■ 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1

-1 — !

. ! 1 I

- 1! - e0-H1 Error .

-s ------------- e£' 1-H1 Error -

- -----e£' 2-H1 Error -

- vV_________________

....................

time(t)

(c) H1 error of QE (£ = 1/8)

(d) H1 error of e£ (£ = 1/16)

time(t)

time(t)

(e) L2 error of u£ (£ = 1/8)

(f) L2 error of u£ (£ = 1/16)

lu 0.4

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 - 1 . . -

% - 1 -

! - li li -

- li - li - I 'i - li - - \\ ----------- - U0-H1 Error " - ui1-H1Error -

- I\ - w - u1'2-H1Error "

....................

time(t)

(g) H1 error of u£(£ = 1/8)

lu 0.4

- U0-H1Error

------------- ui1 - H1 Error

-----u1'2 - Hi Error

---------------------------------------------------

time(t)

(h) H1 error of u£(£ = 1/16)

Fig.13. Time history of computed relative errors of the temperature increment and displacement with £ = 1/8

and 1/16.

At a particular time t = 5, the computed relative errors of the temperature increment and the radial displacement in both L and H1 norms are listed in Tabs. 7 and 8, and the converging curves are plotted in Fig.14 with decreasing £ . Note that at this special time, the precisions of u^1 are better with smaller errors

than those of u£'2 in LL norm, as seen in Tab.8, while the errors of u^1 in H1 norm, which are more

important in the engineering, are still too large. Therefore, the second-order correctors are necessary to obtain the correct strain and stress distributions. On the other hand, by observing the numbers in the two tables and the curves in Fig.14, it can be indicated that in the LL norm, the converging order of the relative errors for the homogenized solutions is O(£) and O(£ ) for both FOTS and SOTS approximations, respectively. However, in H1 norm, the converging order is O(1) for the homogenized solutions as they keep almost unchanged, and O(£) for the FOTS approximations. For the SOTS approximate solutions in H1 norm, the converging order may be O(£) or higher.

Error estimation of the asymptotic analysis is of fundamental importance, as only if the error estimation is given, the approximate solution can be justified mathematically. However, it involves too much and quite a lot of work need to be done for the dynamic thermo-mechanical coupling problem discussed in this paper. Tedious theoretical manipulations need to be taken, and the computed relative errors shown in Figs. 13, 15 and Tabs.7, 8 may help obtain some inspirations.

Table 7. Computed relative errors of the temperature increments in L2 and H1 norms at t = 5 with

different £.

£ K0!!L2 110% !! 4 !!L2 !!0£!!I} !!0£!!L2 !|0£ llH1 I! 4!!^ l|0£ llH1 !! el!Ih1 l|0£ llH1

1/4 9.3782E-2 1.5047E-2 5.1262E-3 6.1054E-1 6.5283E-2 1.9215E-2

1/8 4.7045E-2 3.8336E-3 9.1295E-4 6.2728E-1 3.2910E-2 5.3504E-3

1/16 2.3530E-2 9.5347E-4 1.5841E-4 6.3173E-1 1.6553E-2 1.6640E-3

1/32 1.1757E-2 2.3599E-4 3.5320E-5 6.3286E-1 8.2847E-3 5.3110E-4

1/64 5.8226E-3 5.9027E-5 8.0048E-6 6.3314E-1 4.1237E-3 2.3896E-4

1/128 2.9413E-3 1.3831E-5 2.7411E-6 6.3321E-1 2.0169E-3 1.0638E-4

Table 8. Computed relative errors of radial displacements in L2 and H1 norms at t = 5 with different

£ II Ob 114% KIIL2 II u£ IIL2 II IIL2 II U£ IIL2 IkX II u£ IIH1 II ex KIIh> ikX> llux

1/4 6.8089E-2 1.2935E-2 1.6769E-2 7.3690E-1 8.3650E-2 2.7007E-2

1/8 3.3887E-2 3.1970E-3 4.1717E-3 7.3310E-1 4.3204E-2 6.9798E-3

1/16 1.6888E-2 7.8824E-4 1.0325E-3 7.3239E-1 2.1677E-2 1.9663E-3

1/32 8.4380E-3 1.9561E-4 2.5687E-4 7.3223E-1 1.0836E-2 3.8983E-4

1/64 4.2182E-3 4.8500E-5 6.4105E-5 7.3219E-1 5.3887E-3 1.0258E-4

1/128 2.1091E-3 1.1718E-5 1.5967E-5 7.3218E-1 2.6299E-3 4.5009E-5

(a) Temperature increment. (b) Radial displacement.

Fig. 14. Converging curves of the relative errors.

Fig.15 shows the approximate solutions of the radial and hoop stress with £ = 1/16 at t = 1.0 ,3.0,7.0 and 9.0, from which the discontinuity of the stresses on the interface of the two materials is clearly observed. As predicted in comparing the relative errors of radial displacement in Fig.12 and Tab.8, the FOTS approximate solution can not give the correct results, and the biggest errors lie around the two boundaries, due to the errors

of the radial displacement in H1 norm. By adding the second-order correctors, the satisfactory stress distributions are obtained at almost every point of the domain.

(a) af at t = 1.0

(b) al at t = 1.0

0.55 0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 0

(c) a£ at t = 3.0

(d) al at t = 3.0

(e) af at t = 7.0

(f) &L at t = 7.0

0.25 0.2 0.15 0.1 0.05 0

-0.05 -0.1 -0.15 -0.2

i •i

: \ zJl.

,. Ij I. »

!| Ï < t i

' ? || >1 ji !i

' Ï !i !i

1 ' Î (

_ „ e, 2

1.4 1.6

(g) p at t = 9.0

(h) OÎ at t = 9.0

Fig. 15. Approximate solutions of radial stresses and hoop stresses at different times with e = 1/16.

Finally, let us discuss the oscillating behavior of the system in different coupling ratio rc . Set T0 = 25,50,100 and we have r0 = 0.1941,0.3856,0.7609, respectively. Fig.16 shows the time history of the temperature increment and displacement at r = 1.5 with £ = 1/8. It is found that the curves converge faster and the whole system approaches to the steady state earlier with bigger r^. Actually, the dynamic

thermo-mechanical coupling of the structure is a two-way energy transformation process. The heat is transferred by conduction, part of which is converted into the elastic energy and causes the oscillation of the displacement if inertia force is considered. Conversely, part of the elastic energy is converted into heat and further transferred by conduction. With this continuously transferring process, the final steady state is reached with the energy balance between inside and outside of the system. In the example herein, the heat generated by the heat source is partly used to increase the temperature of the body and partly produces the mechanical oscillating. On the other hand, the straining of the structure produces extra heat which takes participate in the heat conduct process, causes the dissipation of the elastic energy and simultaneously suppresses the heat transferring. This bidirectional energy transformation leads to the lag effect of both the temperature and displacement field. As rc0 increases, the coupling effect is more remarkable and the system reaches static

more quickly. Conversely, with smaller r£°, the oscillating of the displacement will be kept for a longer time. In the extreme condition, if rc0 = 0 , which implies the energy transformation only occurs from heat to elasticity and without the damping effect, the thermal-induced oscillating will always exist and never dissipate.

5 10 15 20 25 30 35 40

time(t)

(a) Temperature increment

10 20 time(t)

(b) Radial displacement

Fig. 16. Time history of the 0£ and u£ at r = 1.5 ( £ = 1/8) with different rc°.

6. Concluding remarks

In this work, the governing equations of the dynamic thermo-mechanical coupling problems for the axisymmetric and spherical symmetric structures are presented. The SOTS expansions of the temperature increment and displacement are formally defined and the homogenized coefficients and microscopic cell problems are deduced and discussed. Then, the second-order finite element solutions of the temperature increment and displacement for the dynamic thermo-mechanical coupling problems in the axisymmetric and spherical symmetric structures are presented on the basis of the two-scale asymptotic analysis. Because of the existence of the extra terms in the circumferential direction in the coupled model, more cell functions and homogenized coefficients are obtained in comparison with the conventional coupling system in Cartesian coordinate [29]. Generally, it is indicated that the circumferential strength of the axisymmetric and spherical symmetric structure is enhanced with the bigger Young's modulus, smaller Possion's ratio, and smaller thermal expansion by the homogenization results. In the second part of this paper, we use the SOTS finite element algorithm to study the axisymmetric quasi-static and dynamic thermo-mechanical coupling problems where the results of the present SOTS finite-element algorithm match the FE fine solutions with hardly any difference and the accuracy and effectiveness of the present algorithm is verified.

It has been indicated from the present study that the computational meshes in the cell domain and the homogeneous domain are much coarser than the original periodic fine mesh. With the same heat and elastic stiff matrices in the left hand side in the algebraic equation for solving the cell problems, the cell functions are

obtained very fast. The homogenized solutions are also obtained quite quickly because of the homogenized coefficients, which make the iteration matrices well-conditioned by the comparison with the original matrices with rapid oscillating coefficients. Since the cell functions and the homogenized solutions are all independent of the periodicity £ , by the present SOTS asymptotic expansion in Eq.(26), the original thermo-mechanical coupling problem with arbitrary can be reproduced, which is much convenient for simulating and predicting the thermal and mechanical behaviors of the composite materials.

On the other hand, the computation scheme proposed in this paper gives out good performance to simulate the vibrations of the temperature increment and displacement in a long time. With the synchronization between the SOTS asymptotic analysis and the FE solution in the fine mesh, the very detailed information in the thermo-mechanical system can be captured well. It is also indicated that only by adding the second-order correctors, the stress distributions can be accurately obtained. The ratio r£° is an important factor to evaluate

the coupling effect of the two physical fields. With increasing rc°, the lag effects for the temperature and

displacement field become more significant.

As this work is only the beginning of the computable modeling of the dynamic thermo-mechanical coupling problems, in the future, the present method will be further developed and applied to the simulations of the spacecraft dynamic thermo-mechanical coupling response, structural deformation and failure under reentry aerothermodynamic environment.

Acknowledgements

The research is supported by National Key Basic Research and Development Program (2014CB744100) and National Natural Science Foundation of China (11325212,91530319), and also supported by China Postdoctoral Science Foundation (2014M562616). The first author would like to thank Francesca Gulinatti for her useful help on proofreading the problem of language of this paper.

References

[1] J. Auriault, Effective macroscopic description for heat conduction in periodic composites, Int. J. Heat. Mass. Tran. 26 (1983) 861-869.

[2] A. Norris, A differential scheme for the effective moduli of composites, Mech. Mater. 4 (1985) 1-16.

[3] M. Jiang, I. Jasiuk, M. Ostoja-Starzewski, Apparent thermal conductivity of periodic two-dimensional

composites, Comp. Mater. Sci. 25 (2002) 329-338.

[4] A. Bensoussan, J.L. Lions, G. Papanicolaou, Asymptotic analysis for periodic structures, North - Holland, Amsterdam, 1978.

[5] E. Sanchez-Palencia, Non-homogenous media and vibration theory, Lecture Notes in Physics, Vol. 127, Springer, Berlin, 1980.

[6] O.A. Oleinik, A.S. Shamaev, G.A. Yosifian, Mathematical problems in elasticity and homogenization, North-Holland, Amsterdam, 1992.

[7] G. Allaire, Homogenization and two-scale convergence, SIAM J. Math. Anal. 23(6) (1992) 1482-1518.

[8] VV Jikov, S.M. Kozlov, O.A. Oleinik, Homogenization of differential operators and integral functional, Springer Verlag, Berlin-Heidelberg, 1994.

[9] D. Cioranescu, P. Donato, An Introduction to Homogenization. Oxford University Press, New York, 1999.

[10] G. Allaire, K.E. Ganaoui, Homogenization of a conductive and radiative heat transfer problem, Multiscale Model. Simul. 7(3) (2009) 1148-1170.

[11] G. A. Francfort, Homogenization and linear thermoelasticity, SIAM. J. Math. Anal. 14 (1983) 696-708.

[12] W. J. Parnell, Coupled thermoelasticity in a composite half-space, J. Eng. Math. 56 (2006) 1-21.

[13] V. L. Savatorova, A.V Talonov, A. N. Vlasov, Homogenization of thermoelasticity processes in composite materials with periodic structure of heterogeneities, A. Angew. Math. Mech. 93(8) (2013) 575-596.

[14] N. Takano, Y. Ohnishi, M. Zako, K. Nishiyabu, The formulation of homogenization method applied to large deformation problem for composite materials, Int. J. Solids Struct. 37 (2000) 6517-6535.

[15] G. Chatzigeorgiou, N. Charalambakis, F. Murat, Homogenization problems of a hollow cylinder made of elastic materials with discontinuous properties, Int. J. Solids Struct. 45 (2008) 5165-5180.

[16] T.Y. Hou, X.H. Wu, A multiscale finite element method for elliptic problems in composite materials and porous media, J. Comput. Phys. 134 (1997) 169-189.

[17] W.N. E, B. Engquist, The heterogeneous multi-scale methods, Comm. Math. Sci. 1(1) (2003) 87-132.

[18] W.N. E, P. B. Ming, P. W. Zhang, Analysis of the heterogeneous multiscale method for elliptic homogenization problems, J. Amer. Math. Soc., 18(1) (2005) 121-156.

[19] J.Z. Cui, L.Q. Cao, Finite element method based on two-scale asymptotic analysis, Math. Numer. Sinica. 1 (1998) 89-102.

[20] J.Z. Cui, T.M. Shin, Y.L. Wang. The two-scale analysis method for bodies with small periodic configurations, Struct. Eng. Mech. 7(6) (1999) 601-604.

[21] L.Q. Cao, Multiscale asymptotic expansion and finite element methods for the mixed boundary value problems of second order elliptic equation in perforated domains, Numer. Math. 103(1) (2006) 11-45.

[22] F. Su, J.Z. Cui, Z. Xu, Q.L. Dong, A second-order and two-scale computational method for the quasi-periodic structures of composite materials, Finite Elem. Anal. Des. 46(2010) 320-327.

[23] Q. Ma, J.Z. Cui, Second-Order Two-Scale Analysis Method for the Quasi-Periodic Structure of Composite Materials under Condition of Coupled Thermo-Elasticity, Adva. Mater. Res. 629(2013) 160-164.

[24] Z.Q. Wang, J.Z. Cui, Second-order two-scale method for bending behavior analysis of composite plate with 3-D periodic configuration and its approximation, Sci. China Math. 57(8)(2014) 1713-1732.

[25] Z.Q. Yang, J.Z. Cui, Y. Sun, J.R. Ge, Multiscale computation for transient heat conduction problem with radiation boundary condition in porous materials, Finite Elem. Anal. Des. 102-103(2015) 7-18.

[26] Q.F. Zhang, J.Z. Cui. Existence theory for Rosseland Equation and its homogenized equation, Appl. Math. Mech. 33(12) (2012)1595-1612.

[27] Y.P. Feng, J.Z. Cui, Multi-scale analysis and FE computation for the structure of composite materials with small periodic configuration under condition of coupled thermo-elasticity, Int. J. Numer. Meth. Eng. 60 (2004) 241-269.

[28] J.J. Wan, Multi-scale Analysis Method for Dynamic Coupled Thermoelasticity of Composite Structures, PhD thesis. The Academy of Mathematics and Systems Science, CAS, China. 2007.

[29] Z.H. Yang, J.Z. Cui, Y.T. Wu, Z.Q. Wang, J.J. Wan, Second-order Two-scale analysis method for dynamic thermo-mechanical problems in periodic structure, Int. J. Numer. Anal. Model. 12(1) (2015) 144-161.

[30] Y.Y. Li, J.Z. Cui, The multi-scale computational method for the mechanics parameters of the materials with random distribution of multi-scale grains, Compos. Sci. Technol. 165 (2005) 1447-1458.

[31] Z.H. Yang, Y. Chen, Z.Q. Yang, Q. Ma, Dynamic thermo-mechanical coupled response of random particulate composites: A statistical two-scale method, Chin. Phys. B. 23(7) (2014) 605-616.

[32] Q. Ma, J.Z. Cui, Z.H. Li, Second-order two-scale analysis algorithm for axisymmetric and spherical symmetric elastic problems of periodic composite materials, Int. J. Solids Struct. 78-79(2016) 77-100.

[33] J.P. Carter, J.R. Booker, Finite element analysis of coupled thermoelasticity, Comput. Struct. 31 (1) (1989) 73-80.