Second-order asymptotic algorithm for heat conduction problems of periodic composite materials in curvilinear coordinates

Journal of Computational and Applied Mathematics 306 (2016) 87–115

Contents lists available at ScienceDirect

Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam

Second-order asymptotic algorithm for heat conduction problems of periodic composite materials in curvilinear coordinates Qiang Ma a,b , Junzhi Cui c , Zhihui Li a,b,∗ , Ziqiang Wang d a

Hypervelocity Aerodynamics Institute, China Aerodynamics Research and Development Center, Mianyang 621000, China

b

National Laboratory for Computational Fluid Dynamics, BUAA, Beijing 100191, China

c

LSEC, ICMSEC, Academy of Mathematics and Systems Science, CAS, Beijing 100190, China

d

College of Science, Guizhou Minzu University, Guiyang 550025, China

article

info

Article history: Received 6 May 2015 Received in revised form 21 February 2016 Keywords: Asymptotic analysis Heat conduction equation Quasi-periodic problem Coordinate transformation Finite element computation Error estimation

abstract A new second-order two-scale (SOTS) asymptotic analysis method is presented for the heat conduction problems concerning composite materials with periodic configuration under the coordinate transformation. The heat conduction problems are solved on the transformed regular domain with quasi-periodic structure in the general curvilinear coordinate system. By the asymptotic expansion, the cell problems, effective material coefficients and homogenized heat conduction problems are obtained successively. The main characteristic of the approximate model is that each cell problem defined on the microscopic cell domain is associated with the macroscopic coordinate. The error estimation of the asymptotic analysis method is established on some regularity hypothesis. Some common coordinate transformations are discussed and the reduced SOTS solutions are presented. Especially by considering the general one-dimensional problem, the explicit expressions of the SOTS solutions are derived and stronger error estimation is presented. Finally, the corresponding finite element algorithms are presented and numerical results are analyzed. The numerical errors presented agree well with the theoretical prediction, which demonstrate the effectiveness of the second-order asymptotic analysis method. By the coordinate transformation, the asymptotic analysis method can be extended to more general domain with periodic microscopic structures. © 2016 Published by Elsevier B.V.

1. Introduction Composite materials are made from two or more constituent materials with significantly different physical or chemical properties, that when combined, produce a material with characteristics different from the individual components. They can overcome the defects of a single material and expand the application scope. With the rapid development of aeronautic and aerospace engineering, composites are widely used for manufacturing aircraft wing, satellite antenna and its supporting structure, solar wing and shell, and the Thermal Protection System (TPS) for the space vehicle, etc., because of the good thermal stability, high specific stiffness and strength. The composites are mostly heterogeneous with complex micro structures. Developing an effective method to simulate and predict the thermal and mechanical performance is of importance and practical significance and is necessary for the design and manufacture of the composites.

∗

Corresponding author at: Hypervelocity Aerodynamics Institute, China Aerodynamics Research and Development Center, Mianyang 621000, China. E-mail addresses: [email protected] (Q. Ma), [email protected] (J. Cui), [email protected] (Z. Li), [email protected] (Z. Wang).

http://dx.doi.org/10.1016/j.cam.2016.04.007 0377-0427/© 2016 Published by Elsevier B.V.

88

Q. Ma et al. / Journal of Computational and Applied Mathematics 306 (2016) 87–115

The evaluation of the heat transferring for the composite materials mainly involves studying the properties of the second-order elliptic equations with rapidly oscillating coefficients and many related basic scientific problems, including multi-scale correlation model, and high-performance multi-scale algorithms should be considered. In the recent decades, the researches on the multi-scale algorithm have attracted the attention of many authors. Mathematically, Bensoussan et al. [1] introduced the asymptotic expansion, and applied this idea to various physical and mechanical problems, and gave the theoretical results of the converging behavior. Oleinik et al. [2] studied the homogenization problem and proved the well-known ε 1/2 error estimation. The homogenization method is investigated and extended to the steady or transient heat conduction problem, wave propagation [3], optimal control [4], and different differential operators and integral functions [5], and perforated structures [6]. The vibration theory of the heterogeneous media is considered by Sanchez-Palencia [7]. Allaire [8] considered the homogenization problem by introducing the well-known two-scale convergence method, and recently developed the asymptotic expansion of a coupled conduction–radiation heat transfer problem [9]. Based on this research, various multi-scale methods have been proposed [10,11], but only the first-order asymptotic expansions are considered. Recently, Cui and Cao et al. [12–17] introduced the Second-Order Two-Scale (SOTS) analysis method to predict physical and mechanical behaviors of composite and perforated materials more accurately by considering the second-order correctors. Especially by choosing proper representative cell domain, higher-order convergence and more accurate results were obtained [18]. Su and Ma [19–21] extended the SOTS method to heat conduction, linear elasticity and thermo-elastic problems for quasi-periodic structures. Wang [22] developed the SOTS method for bending behavior analysis of composite plate. Another remarkable work for the heat conduction of the periodic or porous composite material is the SOTS asymptotic expansion integrating conduction, convection and radiation is obtained by Yang, Zhang and Ma [23–26]. Incorporated with the thermo-elastic response of the composites, Feng and Wan [27,28] obtained the multi-scale second-order solutions for the quasi-static and dynamic thermo-elastic problems. Yang [29] used the statistical SOTS method to compute the dynamic thermo-mechanical coupled response of random particulate composites. According to most literatures, although the asymptotic studies of the heat conduction problem are extended to the porous materials [18], random distributed structure [29], the coupled conduction–radiation condition [9,24,25], etc., the numerical simulations were often tested and carried out only in the regular domains, rectangle or cuboid for example, and did not extend to more general and irregular domains. As we all know, the structures made of composites are not always regular, and the cylinder, ball, or fan-shaped structures, etc. are often used in the engineering. The hyperbolic plates and shells are applied in the TPS for the space vehicle with complicated microscopic structures and the SOTS method cannot be applied to these kinds of structures directly due to their curved profiles. Based on the idea of coordinate transformation that is often adopted in the finite difference computation, we try to find a proper transformation from the regular domain to this particular domain so as to apply the SOTS method in the regular domain and obtain the asymptotic solutions. Then, the heat transfer problem in the original domain becomes quasi-periodic problem after the transformation. Finally, we obtain the real asymptotic solutions by the inverse transformation. The mesh generation method can also be used to make the direct mapping between the two domains. By using this idea, the SOTS method can be extended to more arbitrary domains. The remaining part of this paper is outlined as follows: The heat conduction problem in the curvilinear coordinate system and the coordinate transformation are presented in Section 2. The SOTS analysis method is discussed in Section 3, including the cell problems, homogenized solutions and coefficients, the error estimation and some specific transformations. The finite element algorithm is proposed in Section 4. Three numerical examples are tested and discussed in Section 5, followed by conclusions and future work in Section 6. Throughout this paper, the common notations of Sobolev space and the convention of summation on repeated indices are used, and the letters in bold represent the matrix or vector functions in the formulation. By O(ε k ), k ∈ N, we denote there exists a constant c independent of ε such that |O(ε k )| ≤ c ε k . 2. Governing equation Without loss of generality, the coordinate transformation is illustrated in Fig. 1 in the two-dimensional domain. The structure occupying the domain Ω in Fig. 1(a) is quasi-periodically distributed in the conventional Cartesian coordinate x, which is to say that although we can choose a representative cell Υ shown in Fig. 1(c) to characterize the microscopic feature of the structure, the cells in the macroscopic domain may be bent or twisted such that all of them are different from each other. For this structure, the asymptotic expansion method [1–4] cannot be adopted directly, as it is only effective when the cells in the composite materials are periodically arranged along the axis. However, by choosing proper coordinate transformation,

ξ = L (x),

(1)

this quasi-periodic domain Ω is transformed into periodic domain Θ shown in Fig. 1(b) in the general curvilinear coordinate ξ . Then the asymptotic expansion method can be applied to analyze the heat conduction performance of the structure, and after obtaining the asymptotic expansion of the temperature, by the inverse transformation x = L −1 (ξ),

(2)

the asymptotic behavior of the temperature in the original domain Ω can be simulated and predicted. To make this transformation well defined, we assume the transformation (1) is proper, i.e. the function L (x) is monotonous bounded, continuous and its first partial derivative is also continuous, in which case the inverse transformation always exists.

Q. Ma et al. / Journal of Computational and Applied Mathematics 306 (2016) 87–115

89

Fig. 1. Coordinate transformation from a quasi-periodic domain Ω (a) to a regularly periodic domain Θ (b) with normalized cell domain Υ (c).

Based on this idea, firstly, assume that each of the materials that constitutes the composite is homogeneous and isotropic for simplicity, and consider the following static heat conduction problem in the domain Ω ∈ Rn in Cartesian coordinates,

   ∂ ∂ T ε (x)  ε  − a ( x ) = f (x) in Ω ,  ∂ xi ∂ xi  ∂ T ε (x) T ε (x) = 0 on Γ ,  aε (x) nxi = q(x) on Γ2 1 ∂ xi

(3)

in which n = 1, . . . , 3 denotes the dimension, aε (x) denotes the heat conductivity coefficient and ε is a small parameter related to the periodicity. The temperature is denoted by T ε (x) and the heat source term by f (x). The x = (xi )1≤i≤n is the coordinate variable. The boundary ∂ Ω is smooth and

∂ Ω = Γ1 ∪ Γ2 ,

Γ1 ∩ Γ2 = ∅.

nx = (nxi )1≤i≤n denotes the outward unit normal to Γ2 and q(x) is the imposed heat flux. There is some periodicity for the coefficients aε (x) in Ω and by the coordinate transformation (1), the problem (3) will be rewritten in the general curvilinear coordinate ξ = (ξi )1≤i≤n ∈ Θ . From Eq. (1), the partial derivative is changed by

∂ ∂ξj ∂ → . ∂ xi ∂ξj ∂ xi

(4)

90

Q. Ma et al. / Journal of Computational and Applied Mathematics 306 (2016) 87–115

Taking Eq. (4) into the equation in (3), we obtain the corresponding problem as

   1 ∂ ∂ T ε (ξ) ε   − J (ξ) ∂ξ Cij (ξ)a (ξ) ∂ξ J (ξ) = f (ξ) in Θ , i

j

  T ε (ξ) = 0 on Γ¯ 1 ,

E (ξ)Cij (ξ)aε (ξ)

∂ T ε (ξ) nξi = q(ξ) on Γ¯ 2 ∂ξj

(5)

where Cij is the adjoint metric tensor expressed as Cij (ξ) =

∂ξi ∂ξj , ∂ xk ∂ xk

and J is the Jacobian defined by

 ∂(x1 , . . . , xn )  . J (ξ) =  ∂(ξ1 , . . . , ξn )  



From the assumption on L (x), J > 0 in Θ . Moreover, assume that Cij is uniformly bounded, and L (x) ∈ (H 2 (Ω ))n , then we have the relationship 1 ∂ J ∂ξi

(Cij J ) = Dj (ξ) =

n  ∂ 2 ξj k=1

∂ x2k

.

The boundaries Γ¯ 1 and Γ¯ 2 are the transformed boundaries of Γ1 and Γ2 , respectively (see Fig. 1(b)), and also,

∂ Θ = Γ¯ 1 ∪ Γ¯ 2 ,

Γ¯ 1 ∩ Γ¯ 2 = ∅.

nξ = (nξi )1≤i≤n denotes the outward unit normal to Γ¯ 2 and E is the normalized factor with the expression as

   n   ∂ξk 2 2 nξk . E (ξ) = 1  ∂ xi i=1 For instance, in the three dimensional case, if ξ1 is constant on Γ¯ 2 , we have nξ = (±1, 0, 0),

E (ξ) = 1/|∇x ξ1 |,

where ∇x denotes the gradient operator with respect to x. Now the coefficient aε (ξ) in Θ is periodic and can be defined by

  ξ a (ξ) = a = a(η), ε ε

0 < a ≤ aε ≤ a¯ ,

(6)

where η = (ηi )1≤i≤n ∈ Υ is the microscopic variable opposed to ξ . The periodic cell domain is represented by Υ = [0, 1]n and a(η) is called a 1-periodic function with a, a¯ ∈ R as its lower and upper boundary. Moreover, suppose that Θ is composed of entire cells with periodicity ε , i.e.,

Θ=



ε(Υ + z ),

z ∈I

where I = {z = (zi )1≤i≤n ∈ Zn |ε(Υ + z ) ∈ Θ } is the index set. It is also assumed that Γ¯ 1 and Γ¯ 2 consist of the entire cell surfaces. The domain Θ then is regular such that the SOTS method can be applied and the SOTS asymptotic expansion for T ε (ξ) can be obtained. Finally, the SOTS approximation for real temperature T ε (x) is obtained by the inverse transformation equation (2). 3. Second-order two-scale asymptotic analysis 3.1. Second-order two-scale expansion of temperature Based on the heat problem (5) in curvilinear coordinate ξ , we use the SOTS method for the temperature field T ε (ξ). Firstly, assume T ε (ξ) is formally expanded as follows, T ε (ξ) = T0 (ξ) + ε T1 (ξ, η) + ε 2 T2 (ξ, η) + O(ε 3 ),

(7)

where T0 (ξ) is the homogenized solution that does not depend on the microscopic variable η. T1 (ξ, η) and T2 (ξ, η) are called the first- and second-order correctors. Respecting ξ and η as two independent variables, the partial derivative becomes

∂ ∂ ∂ → + ε −1 . ∂ξi ∂ξi ∂ηi

(8)

Q. Ma et al. / Journal of Computational and Applied Mathematics 306 (2016) 87–115

91

Take Eqs. (7) and (8) into the heat equation (5), rewriting the equation above as the ε power series, and equating the coefficients of the same power of ε , one obtains the first three equations, O(ε −2 ) : −Cij

∂ ∂ηi

 a

∂ T0 ∂ηj



= 0,

(9)

      ∂ ∂ ∂ ∂ T1 ∂ T0 ∂ T0 ∂ T0 O(ε ) : −Cij a − Cij a − Cij a − Dj a = 0, ∂ηi ∂ηj ∂ηi ∂ξj ∂ξi ∂ηj ∂ηj         ∂ ∂ ∂ ∂ ∂ T2 ∂ T1 ∂ T1 ∂ T0 O(ε 0 ) : −Cij a − Cij a − Cij a − Cij a ∂ηi ∂ηj ∂ηi ∂ξj ∂ξi ∂ηj ∂ξi ∂ξj   ∂ T0 ∂ T1 + = f. − Dj a ∂ηj ∂ξj −1

(10)

(11)

From the definition of T0 (ξ), it follows that Eq. (9) holds naturally. Using this property, Eq. (10) reduces to

−Cij

∂ ∂ηi

 a

∂ T1 ∂ηj

 = Cij

∂ a ∂ T0 . ∂ηi ∂ξj

Setting formally T1 (ξ, η) = Nα1 (ξ, η)

∂ T0 (ξ) , ∂ξα1

α1 = 1, . . . , n,

(12)

where Nα1 (ξ, η) is the 1-periodic function in η direction that satisfies the cell problem

   ∂ Nα1 (ξ) ∂ a(η) ∂   a (η) = Ckα1 (ξ) in Υ , − C (ξ)  ij ∂ηi ∂ηj ∂ηk     Nα1 (ξ, η) = 0.

(13)

Υ

1 1 Then, Nα1 (·, η) ∈ Hper (Υ ) with the function space Hper (Υ ) defined as



1 Hper (Υ ) =

ϕ|ϕ ∈ H 1 (Υ ), ϕ(η) 1 − periodic ,

 Υ

 ϕ dη = 0 .

From the cell problem Eq. (13), we have the following theorem. Theorem 1. For two different ξ , ξ ∈ Θ , let 0

Nα01 (ξ, η) = Nα1 (ξ , η0 ), 0

1

Nα11 (ξ, η) = Nα1 (ξ , η1 ), 1

and the following inequality holds

∥Nα01 − Nα11 ∥H 1 (Υ ) ≤ c |ξ 0 − ξ 1 |,

(14)

where the constant c is independent of ξ and ξ . 0

1

Proof. From the cell problem (13), the weak forms for Nα01 and Nα11 are



∂ Nα01 ∂ϕ dη = − ∂ηj ∂ηi



∂ Nα11 ∂ϕ 1 Cij (ξ )a dη = − ∂ηj ∂ηi



Cij (ξ )a 0

Υ

 Υ

with ϕ ∈

1 Hper

Cij (ξ ) 0

Ckα1 (ξ )a

∂ϕ dη, ∂ηk

Ckα1 (ξ )a

∂ϕ dη ∂ηk

0

Υ

1

Υ

(Υ ) being the test function. By subtraction, there is

 a Υ

∂(Nα01 − Nα11 ) ∂ϕ 0 1 dη = −[Cij (ξ ) − Cij (ξ )] ∂ηj ∂ηi

 a Υ

∂(Nα11 − ηα1 ) ∂ϕ dη. ∂ηj ∂ηi

92

Q. Ma et al. / Journal of Computational and Applied Mathematics 306 (2016) 87–115

From the continuity of Cij and boundedness of Cij and Nα11 , taking ϕ = Nα01 − Nα11 and using the Schwartz inequality gives c1 ∥Nα01 − Nα11 ∥2H 1 (Υ ) ≤ c2 |ξ − ξ | · ∥Nα01 − Nα11 ∥H 1 (Υ ) , 0

1

where the constants c1 and c2 are independent of ξ and ξ , finally we obtain 0

1

∥Nα01 (ξ 0 , η0 ) − Nα11 (ξ 1 , η1 )∥H 1 (Υ ) ≤ c |ξ 0 − ξ 1 |. From Theorem 1, it is seen that generally because of the existence of Cij , 1 (Υ ), Nα1 (ξ, η) ∈ H 1 (Θ ) × Hper

which means that each macroscopic point is associated with a particular cell function Nα1 . If Cij is constant or the macroscopic parameters in Eq. (13) can be separated, then Nα1 will be independent of ξ and can be more easily solved both analytically and numerically. Next, considering Eq. (11), we make an integral average on Υ . Since Cij is always the macroscopic term, referring the expression T1 in Eq. (12), the homogenized equation for T0 (ξ) is obtained as follows

−

1 ∂



J ∂ξi

∂ T0 J ∂ξk

Cij a0jk



= f,

(15)

with the so called ‘‘homogenized’’ or ‘‘effective’’ heat conductivity a0ij (ξ) defined by a0ij (ξ) =

1

|Υ |

  Υ

a(η)δij + a(η)

 ∂ Nj (ξ, η) dη, ∂ηi

(16)

where |Υ | denotes the measure of the domain Υ and we can immediately have a0ij (ξ) ∈ H 1 (Θ ). Based on Eq. (15), we defined the homogenized heat conduction problem for T0 (ξ) as

   ∂ T0 (ξ) 1 ∂  0  − C (ξ) a J (ξ) = f (ξ) in Θ , (ξ) ij  jk J (ξ) ∂ξi ∂ξk  ∂ T0 (ξ)  T0 (ξ) = 0 on Γ¯ 1 , E (ξ)Cij (ξ)a0jk (ξ) nξi = q on Γ¯ 2 . ∂ξk

(17)

Remark 1. Now we have obtained the first two expansions of T ε (ξ), i.e., T0 + ε T1 , if we take the truncated function T ε − (T0 + ε T1 ) into heat equation (5),

−

1 ∂ J ∂ξi



     ∂(T ε − T0 − ε T1 ) 1 ∂ ∂ T0 ∂ T1 Cij a(ξ) J = f + Cij a + J ∂ξj J ∂ξi ∂ξj ∂ηj     ∂ T1 1 ∂ ∂ T1 1 ∂ Cij a J +ε Cij a J . + J ∂ηi ∂ξj J ∂ξi ∂ξj

It can be noted that the residual on the right hand side is of order O(1) that does not equal to 0. In the practical engineering computation, it cannot be omitted for a constant ε , so engineers conclude that the first-order two-scale approximate solution T0 + ε T1 cannot be accepted and the micro-scale fluctuation of the temperature is far from being captured. This is the reason why it is necessary to consider the higher order expansions. Now we resolve to obtain the expression of the second-order corrector T2 (ξ, η). Taking Eqs. (12) and (15) into Eq. (11), we have

∂ −Cij ∂ηi



∂ T2 a ∂ηj



 ∂ Nα1 ∂ 2 T0 ∂ 0 = Ckα2 (aNα1 ) + Cα2 k a + Cα1 α2 − Cα2 k akα1 ∂ηk ∂ηk ∂ξα1 ∂ξα2        ∂ a0jα1 ∂ Nα1 ∂ 2 Nα1 ∂ Nα1 ∂ T0 ∂ 0 + Dj a a +a , + aδjα1 − ajα1 + Cij − ∂ηj ∂ηi ∂ξj ∂ξi ∂ηj ∂ξi ∂ξα1 

from which T2 (ξ, η) is defined formally as T2 (ξ, η) = Nα1 α2 (ξ, η)

∂ 2 T0 (ξ) ∂ T0 (ξ) + Mα1 (ξ, η) , ∂ξα1 ∂ξα2 ∂ξα1

α1 , α2 = 1, . . . , n,

(18)

Q. Ma et al. / Journal of Computational and Applied Mathematics 306 (2016) 87–115

93

where Nα1 α2 (ξ, η) and Mα1 (ξ, η) are the 1-periodic cell functions and satisfy the following cell problems, respectively.

   ∂ Nα1 (ξ, η) ∂ Nα1 α2 (ξ, η) ∂ ∂   − C (ξ) a (η) = Ckα2 (ξ) (a(η)Nα1 (ξ, η)) + Cα2 k (ξ)a(η)  ij  ∂ηi ∂ηj ∂ηk ∂ηk  +Cα1 α2 (ξ)a(η) − Cα2 k (ξ)a0kα1 (ξ) in Υ ,       Nα1 α2 (ξ, η) = 0.

(19)

     ∂ Mα1 (ξ, η) ∂ Nα1 (ξ, η) ∂  0  −Cij (ξ) a(η) = Dj (ξ) a(η) + a(η)δjα1 − ajα1 (ξ)   ∂ηi  ∂ηj ∂ηj       0  ∂ Nα1 (ξ, η) ∂ 2 Nα1 (ξ, η) ∂ ajα1 (ξ) ∂ +Cij (ξ) a(η) + a(η) − in Υ , ∂ηi ∂ξj ∂ξi ∂ηj ∂ξi          Mα1 (ξ, η) = 0.

(20)

Υ

Υ

1 By the same argument in the proof of Theorem 1, we can also prove that Nα1 α2 and Mα1 ∈ C (Θ ) × Hper (Υ ). ε In summary, the SOTS asymptotic expansion of T (ξ) is obtained as follows,

T ε (ξ) = T0 (ξ) + ε Nα1 (ξ, η)

  ∂ T0 (ξ) ∂ 2 T0 (ξ) ∂ T0 (ξ) + O(ε 3 ). + ε 2 Nα1 α2 (ξ, η) + Mα1 (ξ, η) ∂ξα1 ∂ξα1 ∂ξα2 ∂ξα1

(21)

3.2. Error estimation The error estimation of the expansion is presented in this section. First, let the first- and second-order approximate solutions be T ε,1 (ξ, η) = T0 (ξ) + ε T1 (ξ, η), T ε,2 (ξ, η) = T0 (ξ) + ε T1 (ξ, η) + ε 2 T2 (ξ, η)

(22)

respectively. From Oleinik et al. [2], we have the following two lemmas: Lemma 1. Let Ω be a bounded domain with a smooth boundary and Bε = {x ∈ Ω |ρ(x, ∂ Ω ) < ε, ε > 0}. Then, there exists δ0 such that for every ε ∈ (0, δ0 ) and v ∈ H 1 (Ω ), the following inequality holds

∥v∥L2 (Bε ) ≤ c ε 1/2 ∥v∥H 1 (Ω ) , where c is a constant independent of ε and v . Lemma 2. For w ∈ H 1 (Υ ) and

 ε(∂ Υ )

w 2 dη ≤ c ε

 εΥ



Υ

w dη = 0, the following inequality holds

|∇w|2 dη,

where ∂ Υ is the boundary of Υ and c is a constant independent of ε . Next we prove the following estimation: Theorem 2. Let T ε (ξ) satisfy the heat problem Eq. (5) and have the SOTS expansion in Eq. (21), T0 (ξ) ∈ H 3 (Θ ), Nα1 (ξ, η), Nα1 α2 (ξ, η), Mα1 (ξ, η) ∈ H 2 (Θ ) × H 1 (Υ ) the following inequality holds

∥T ε (ξ) − T ε,2 (ξ, η)∥H 1 (Θ ) ≤ c ε 1/2 ,

(23)

where the constant c is independent of ε . Proof. Throughout the proof, c denotes a constant independent of ε . Taking T ε −T ε,2 into the original equation and boundary condition, we obtain the following boundary value problem for T ε − T ε,2

   ∂(T ε − T ε,2 ) ε 1 ∂    C a J = F ε in Θ , − ij  J ∂ξ ∂ξj J i  ∂(T ε − T ε,2 ) ∂ T0 T ε − T ε,2 = −ε Gε on Γ¯ ,  ECij aε nξi = E α ik nξ − ε E H ε on Γ¯ 2  1 ∂ξj ∂ξk i

(24)

94

Q. Ma et al. / Journal of Computational and Applied Mathematics 306 (2016) 87–115

where

Fε =

 ∂ Nα1 ∂ Nα1 ∂ 2 T0 ∂ T0 J J + Cij a ∂ξj ∂ξα1 ∂ξj ∂ξα1 ∂ξi  ∂ 2 T0 ∂ 3 T0 ∂  + + Cij aNα1 J Cij aNα1 J ∂ξi ∂ξα1 ∂ξj ∂ξα1 ∂ξi ∂ξj   ∂ Nα1 α2 ∂ Nα1 α2 ∂ ∂ 2 T0 ∂ 3 T0 J + Cij a + Cij a J ∂ξi ∂ηj ∂ξα1 ∂ξα2 ∂ηj ∂ξα1 ∂ξα2 ∂ξi     ∂ Mα1 ∂ Mα1 ∂ 2 T0 ∂ Nα1 α2 ∂ 2 T0 ∂ ∂ T0 ∂ + J Cij a J + Cij a J +ε Cij a ∂ξi ∂ηj ∂ξ1 ∂ηj ∂ξ1 ∂ξi ∂ξi ∂ξj ∂ξα1 ∂ξα2     ∂ Mα1 ∂ T0 ∂ ∂ ∂ 3 T0 J +ε +ε Cij aNα1 α2 Cij a J ∂ξi ∂ξα1 ∂ξα2 ∂ξj ∂ξi ∂ξj ∂ξα1     ∂ Nα1 α2 ∂ ∂ 2 T0 ∂ ∂ 2 T0 +ε Cij aMα1 J +ε Cij a J ∂ξi ∂ξα1 ∂ξj ∂ξi ∂ξj ∂ξα1 ∂ξα2   3 ∂ Mα1 ∂ Mα1 ∂ 2 T0 ∂ Nα1 α2 ∂ T0 ∂ ∂ T0 J +ε Cij a J J , + ε Cij a + ε Cij a ∂ξj ∂ξα1 ∂ξα2 ∂ξi ∂ξi ∂ξj ∂ξα1 ∂ξj ∂ξα1 ∂ξi ∂ ∂ξi



Cij a

and on the boundaries

  ∂ 2 T0 ∂ T0 ∂ T0 G = Nα1 + ε Nα1 α2 + Mα1 , ∂ξα1 ∂ξα1 ∂ξα2 ∂ξα1   ∂ Nk , α ik = Cij a0jk − aδjk − a ∂ηj  ∂ Nα1 ∂ T0 ∂ Nα1 α2 ∂ 2 T0 ∂ Mα1 ∂ T0 ∂ 2 T0 H ε = −Cij aε + Nα1 + + ∂ξj ∂ξα1 ∂ξα1 ∂ξj ∂ηj ∂ξα1 ∂ξα2 ∂ηj ∂ξα1   2 3 ∂ Mα1 ∂ T0 ∂ Nα1 α2 ∂ T0 ∂ T0 ∂ 2 T0 + nξi . + Nα1 α2 + Mα1 +ε ∂ξj ∂ξα1 ∂ξα2 ∂ξα1 ∂ξα2 ∂ξj ∂ξj ∂ξα1 ∂ξα1 ∂ξj ε

To estimate the term ∥Gε ∥H 1/2 (Γ¯1 ) , we define the truncate function mε (Θ ) by

 mε ∈ C ∞ (Θ ),    mε = 1, if ρ(ξ, Γ¯ 1 ) ≤ ε, mε (Θ ) = mε = 0, if ρ(ξ, Γ¯ 1 ) ≥ 2ε,    ∥∇ mε ∥L∞ (Θ ) ≤ c . ε The support of mε (Θ ) is a neighborhood of Γ¯ 1 of thickness 2ε in Θ , which we denote by Uε , and set ψ ε = m ε Gε . One has that

∥ψ ε ∥L2 (Uε ) ≤ c , with the constant c independent of ε and for i = 1, . . . , n,

 ∂ Nα1 ∂ T0 ∂ Nα1 α2 ∂ 2 T0 ∂ Mα1 ∂ T0 ∂ψ ε 1 ∂ Nα1 ∂ T0 ∂ 2 T0 = mε + + Nα1 + + ∂ξi ε ∂ηi ∂ξα1 ∂ξi ∂ξα1 ∂ξα1 ∂ξi ∂ηi ∂ξα1 ∂ξα2 ∂ηi ∂ξα1   2 3 2 ∂ Nα1 α2 ∂ T0 ∂ Mα1 ∂ T0 ∂ T0 ∂ T0 +ε + Nα1 α1 + + Mα1 ∂ξi ∂ξα1 ∂ξα2 ∂ξα1 ∂ξα2 ∂ξi ∂ξi ∂ξα1 ∂ξα1 ∂ξi    2 ∂ mε ∂ T0 ∂ T0 ∂ T0 + Nα1 + ε Nα1 α2 + Mα1 . ∂ξi ∂ξα1 ∂ξα1 ∂ξα2 ∂ξα1 According to Lemma 1

∥T0 ∥H 1 (Uε ) ≤ c ε 1/2 ,

(25)

(26)

(27)

Q. Ma et al. / Journal of Computational and Applied Mathematics 306 (2016) 87–115

95

then, it follows that

∥ψ ε ∥H 1 (Uε ) ≤ c ε −1/2 . Therefore

∥Gε ∥H 1/2 (Γ¯1 ) = ∥ψ ε ∥H 1/2 (Γ¯1 ) ≤ c ∥ψ ε ∥H 1 (Θ ) = c ∥ψ ε ∥H 1 (Uε ) ≤ c ε −1/2 .

(28)

To estimate the term α ik , from the definition (26) and the cell problem (13), it can be easily proved that

 Θ

∂α ik = 0, ∂ηi

α ik dξ = 0,

(29)

and then from the Lemma 2.1 in Oleinik et al. [2],

 Sti

α ik dη = 0,

(30)

where Sti denotes the set Sti = {ξ|ξj = t , 0 ≤ ξi ≤ 1, i ̸= j}. On the outer boundary, ∂ Θ consists of two-dimensional faces of ε(Υ + z ) for some z ∈ I. lj

Denote the two-dimensional faces by σj1 , . . . , σj such that σjk is parallel to the hyperplane ξj = 0 and lies on Γ¯ 2 . We thus have

Γ¯ 2 =

lj n  

σjs .

j =1 s =1

qsj

Denote by



the cell ε(Υ + z ) whose surface contains the set σjs . For any σjs , due to the condition (30), we have

α nξi dη = ik

σjs

 σjs

α jk nξj dη = 0,

(31)

since nξi = 0, for i ̸= j. Note that there is no summation over j in above equation. Set mv (ξ) =

1

|qsj |

 qsj

v dξ ,

ξ ∈ σjs , j = 1, . . . , n, s = 1, . . . , lj .

Therefore, mv (ξ ) is constant on each surface σjs . Taking into account (31) and the fact that no more that 2d qsj have a non-empty intersection, we obtain from Lemma 2

 2  2     ik ik     α n v d η = α n (v − m ) d η ξi ξi v  ¯   ¯  Γ2

Γ2

 ≤c Γ¯ 2

|v − mv |2 dη = c

lj  n   j=1 s=1

σjs

|v − mv |2 dη ≤ c ε

 Θ

|∇v|2 dξ .

(32)

From the equation for T − T ε,2 in Eq. (24), set the function T¯ = T − T ε,2 + ε G¯ ε ,

¯ε

(33)

¯ε

ε

where G ∈ H (Θ ) such that G = G on Γ¯ 1 , and T¯ satisfy the following problem 1

     ∂ T¯ ε ε 1 ∂ ∂ G¯ ε 1 ∂    in Θ , − J ∂ξ Cij a ∂ξ J = J F − ε J ∂ξ Cij a ∂ξ J i j i j  ∂ T¯ ∂ T0 ∂ G¯ ε   ECij aε nξi = E α ij nξi − ε E H ε + ε ECij aε nξ on Γ¯ 2 . T¯ = 0 on Γ¯ 1 , ∂ξj ∂ξj ∂ξj i

(34)

The weak form of Eq. (34) is

 Θ

Cij a

∂ T¯ ∂ Φ Jdξ = ∂ξj ∂ξi

 Θ

ε F ε Φ dξ − ε

 Θ

Cij a

∂ G¯ ε ∂ Φ Jdξ + ∂ξj ∂ξi

 Γ¯ 2

α ij

∂ T0 nξ J Φ d ξ − ∂ξj i

 Γ¯ 2

ε H ε Φ Jdξ ,

(35)

96

Q. Ma et al. / Journal of Computational and Applied Mathematics 306 (2016) 87–115

where Φ ∈ H 1 (Θ ) is the arbitrary test function satisfying Φ = 0 on Γ¯ 1 . Taking Φ = T¯ into Eq. (35), and by Schwarz inequality and the regularity of the coefficients aε and Cij , one gets



            ∂ G¯ ε ∂ Φ   2 ik ∂ T0 ε ε    ¯ Jdξ  +  nξi J Φ dξ  +  |∇ T | dξ ≤ c ε  F Φ dξ | + ε| Cij a α ε H Φ Jdξ  . ∂ξ ∂ξ ∂ξ j i k Θ Θ Θ Γ¯ 2 Γ¯ 2 ∂T

Taking v = ∂ξ0 J T¯ into Eq. (32) gives k

    ik ∂ T0  ≤ c ε1/2 ∥T¯ ∥H 1 (Θ ) .  ¯ J T d ξ n α ξi   ¯ ∂ξk Γ2

(36)

Now, applying again the Schwarz and Poincáre inequalities gives the estimation for T¯ as

∥T¯ ∥H 1 (Θ ) ≤ c ε 1/2 . Taking into account (33) and the residual estimation (28) and (36), we finally get

∥T − T ε,2 ∥H 1 (Θ ) ≤ ∥T¯ ∥H 1 (Θ ) + ∥ − ε G¯ ε ∥H 1 (Θ ) ≤ c ε 1/2 . From the expression above and the cell problems of Eqs. (13), (19) and (20), it is a quasi-periodic problem and the macroscopic variable dependence of the cell functions makes it complicated to find the final approximate solution numerically. Among various coordinate transformations, the simplest is of course the Identity with ξi = xi , and the two domains Ω and Θ are coincident. It can be easily verified that the conventional SOTS expansion of the temperature in Cartesian coordinates [3,18] will be obtained. In the following sections, some other typical coordinate transformations will be discussed to show the adaptability of the asymptotic expansions. 3.3. The general one-dimension problem If the composite is a layered material, for instance, the coefficient depends only on one coordinate ξ1 . And there are C11 = C11 (ξ1 ),

D1 = D1 (ξ1 ),

E = E (ξ1 ),

J = J (ξ1 ).

The source term f = f (ξ1 ), then the heat conduction problem Eq. (5) can be reduced to a one-dimensional problem, from which the explicit expressions of the SOTS expansion can be obtained. Without loss of generality, set the one-dimensional cell domain Υ = [0, 1] and the macroscopic domain

¯ = ε[z1 , z2 ] × Θ ¯, Θ = [a0 , b0 ] × Θ

z1 , z2 ∈ Z, 0 < z1 < z2 ,

¯ denotes the domain in other direction ξi , i = 2, . . . , n. Consider the following general one dimensional problem in which Θ    1 d dT ε (ξ1 )  ε  − C (ξ ) a (ξ ) J (ξ ) = f (ξ1 ) in [a0 , b0 ],  11 1 1 1 J (ξ1 ) dξ1 dξ1 dT ε   T ε (a0 ) = 0; E (ξ1 )C11 (ξ1 )aε (ξ1 ) (b0 ) = q(ξ1 ). dξ1

(37)

From Eq. (21), the SOTS expansion for temperature T ε (ξ1 ) is ε

T (ξ1 ) = T0 (ξ1 ) + ε N1 (ξ1 , η1 )

dT0 (ξ1 ) dξ1

+ε

2



N11 (ξ1 , η1 )

d2 T0 (ξ1 ) dξ12

+ M1 (ξ1 , η1 )

dT0 (ξ1 ) dξ1



+ O(ε3 ),

(38)

where N1 (ξ1 , η1 ) satisfies the following cell problem

  ∂ ∂ N1 da a = C11 ∂η1 ∂η1 dη1  N1 (ξ1 , 0) = N1 (ξ1 , 1) = 0.   −C

11

in Υ = [0, 1],

(39)

It is noted that the cell function satisfies the zero boundary condition instead of the zero integral average, because for the scalar function, the only condition of periodicity is the same value in the two endpoints of the interval Υ . By directly setting the 0 boundary, it is convenient to deal with the Dirichlet boundary condition for T ε,2 . The cell solutions N11 , M1 are also set to have the zero boundary. Remark 2. For the two or three-dimensional cell domain, if the cell domain Υ has the symmetry in the center plane in all axes, the zero boundary condition can also be applied and it can be proven that the cell functions have the H 1 continuity cross the boundary between the neighboring cells. We refer the corresponding work to Cao and Cui [18].

Q. Ma et al. / Journal of Computational and Applied Mathematics 306 (2016) 87–115

97

T0 (ξ1 ) satisfies the homogenized problem

   d ∂ T0 (ξ1 ) 1  0  C (ξ ) a (ξ ) J (ξ ) = f (ξ1 ) in [a0 , b0 ],  11 1 1 1 J (ξ1 ) dξ1 ∂ξ1 dT0   T0 (a0 ) = 0, E (ξ1 )C11 (ξ1 )a0 (ξ1 ) (b0 ) = q, dξ1

(40)

with the coefficient a0 defined by a0 (ξ1 ) = a011 (ξ1 ) =

  Υ

a(η1 ) + a

 ∂ N1 (ξ1 , η1 ) dη, ∂η1

(41)

N11 (ξ1 , η1 ) and M1 (ξ1 , η1 ) satisfy the following cell problems, respectively

    ∂ ∂ N1 ∂ N11 ∂ 0 a = C11 a + (aN1 ) + a − a in Υ = [0, 1], 11 ∂η1 ∂η1 ∂η1 ∂η1  N11 (ξ1 , 0) = N11 (ξ1 , 1) = 0,          ∂ M1 ∂ N1 ∂ N1 ∂ 2 N1 da0 ∂  −C ∂ 0 a = D a a + a + a − a − + C 11 1 11 ∂η1 ∂η1 ∂η1 ∂η1 ∂ξ1 ∂ξ1 ∂η1 dξ1  M1 (ξ1 , 0) = M1 (ξ1 , 1) = 0.   −C

(42)

in Υ = [0, 1],

(43)

From the expressions and equations above, we have the following proposition: Proposition 1. If N1 , N11 and M1 satisfy (39), (42) and (43), respectively, the homogenized coefficient a0 is defined in Eq. (41), then the following statements hold

 1 ∂ N1 (η1 ) 1 (a) a (ξ1 , η1 ) = a (η1 ) = a(η1 ) + a(η1 ) = 1/ dt . ∂η1 a ( t) 0  η1 1 (b) N1 (ξ1 , η1 ) = N1 (η1 ) = −η1 + a0 dt . a ( t) 0    η1  1 t  η1  t 1 1 1 1 1 (c) N11 (ξ1 , η1 ) = N11 (η1 ) = η12 − a0 dsdt − a0 − a0 dsdt dt . 2 2 a(t ) 0 0 a( s ) 0 0 0 a(s) 0

0

(d) M1 = 0.

(44) (45)

(46) (47)

Proof. (a) From the boundedness of Cij , C11 ̸= 0, the equation of N1 in Eq. (39) further reduces to

−

∂ ∂η1

 a

∂ N1 ∂η1



da

=

dη1

.

(48)

Then, N1 (ξ1 , η1 ) = N1 (η1 ). By integration of Eq. (48), we have dN1 dη1

= −1 +

c1 a

,

(49)

where c1 is a constant to be determined. From the periodicity of N1 , integrating Eq. (49) in Υ , thus 1



c1 = 1/

0

1 a(t )

dt .

Taking Eq. (49) into Eq. (41), it follows that a0 =

 



a + a −1 +

c1 

Υ

a

dη = c1 ,

(50)

and (44) is obtained by substituting Eq. (50) into Eq. (49). (b) From Eq. (49), by integrating from 0 to η1 , one gets N1 (η1 ) = −η1 + a0

η1

 0

1 a(t )

dt + c2 ,

where the constant c2 is to be determined. From the boundary condition N1 (0) = N1 (1) = 0, then c2 = 0. (c) From Eq. (44), the equation for N11 in Eq. (42) can be reduced to

−

∂ ∂η1

 a

∂ N11 ∂η1

 =

∂ (aN1 ). ∂η1

(51)

98

Q. Ma et al. / Journal of Computational and Applied Mathematics 306 (2016) 87–115

Therefore, N11 (ξ1 , η1 ) = N11 (η1 ). By integration of Eq. (51), it follows dN11 dη1

= −N1 +

c3 a

,

(52)

with the undetermined constant c3 . By the same argument, integrating Eq. (52) in Υ and taking use of the periodicity of N11 gives c3 = a0

1



N1 (t )dt .

(53)

0

Integrating again Eq. (52) from 0 to η1 , then N11 (η1 ) = −

η1



N1 (t )dt + c3

η1



dt + c4

a(t )

0

0

1

(54)

and taking into account N11 (0) = N11 (1) = 0, we have c4 = 0. Eq. (46) follows by taking the expression (45) of N1 (η1 ) into Eq. (54). (d) From Eq. (44), the right hand side of the equation for M1 in Eq. (43) vanishes, then M1 = 0. Finally, the SOTS expansion for one-dimensional problem Eq. (37) is simplified to T ε (ξ1 ) = T0 (ξ1 ) + ε N1 (η1 )

dT0 (ξ1 ) dξ1

+ ε 2 N11 (η1 )

d2 T0 (ξ1 ) dξ12

+ O(ε3 ),

(55)

and all the cell functions N1 (η1 ), N11 (η1 ) and homogenized coefficients a0 have the explicit expressions by Proposition 1. Moreover, we have the following error estimation: Theorem 3. In the curvilinear coordinate ξ1 , T ε (ξ1 ) satisfies the heat conduction problem (37) and has the SOTS expansion (55), T0 (ξ1 ) ∈ H 3 ([a0 , b0 ]), N1 (η1 ) and N11 (η1 ) ∈ H 1 ([0, 1]), C11 (ξ1 ) ∈ C ([a0 , b0 ]), then the following inequality holds

∥T ε (ξ1 ) − T ε,2 (ξ1 , η1 )∥H 1 ([a0 ,b0 ]) ≤ c ε,

(56)

where c is a constant independent of ε . Proof. By the same argument in the proof of Theorem 2, taking T ε − T ε,2 into the original heat equation and considering the zero boundary condition for the cell functions, we have the following boundary value problem

   d(T ε − T ε,2 ) 1 d   C11 aε J = ε F˜ ε in [a0 , b0 ], − J dξ1 dξ1  d  ε ε, 2 (T − T )(a0 ) = 0; (T ε − T ε,2 )(b0 ) = −ε EC11 G˜ ε EC11 aε dξ1

(57)

where

F˜ ε =



H d

C11

J dξ1

G˜ ε = H

d2 T0 dξ12

d2 T0 dξ12

+ ε aN11

 J

d3 T0 dξ13

+ε ,

1 d J dξ1

 C11 aN11

H = a0

d3 T0 dξ13

 J

,

1



N1 (t )dt = a 0



 ∂ N11 + N1 . ∂η1

Multiply the equation in Eq. (57) by the test function φ ∈ H 1 ([a0 , b0 ]), which satisfies φ(a0 ) = 0, and integrate in the whole macroscopic domain Θ . Note that the variable should be changed to ξ1 and all the integrands do not depend on ξ2 and ξ3 , except the Jacobian J. Finally, we have



b0

C11 a

ε d(T

a0

ε

− T ε,2 ) dφ G(ξ1 )dξ1 = dξ1 dξ1



b0

ε F˜ ε Gφ dξ1 − (εEC11 Gε G¯ φ)|ξ1 =b0 ,

a0

where G(ξ1 ) =

 ¯ Θ

J (ξ1 , . . . , ξn )dξ2 · · · dξn ,

¯ (ξ1 ) = G

 Γ¯ 2

J (ξ1 , . . . , ξn )dξ2 · · · dξn .

Let φ = T ε − T ε,2 and make use of the Schwartz and Poincáre inequalities, we obtain (56). The general one-dimensional problem is discussed in this section, which does not refer to any specific transformation. The simplest is of course, the identity transformation in one dimension. In the following sections, two typical orthogonal coordinate systems are considered, and it can be shown that in the plane axisymmetric and spherical symmetric transformations, the one-dimensional problems obtained are similar to the one discussed herein.

Q. Ma et al. / Journal of Computational and Applied Mathematics 306 (2016) 87–115

99

3.4. Axisymmetric problem The cylindrical coordinate system is shown in Fig. 2(a). Let (ξ1 , ξ2 , ξ3 ) be the cylindrical coordinate instead of the conventional (r , θ , z ) and consider the following transformation x1 = ξ1 cos ξ2 ,

x2 = ξ1 sin ξ2 ,

x3 = ξ3 .

At this time C11 = 1,

C22 =

J = ξ1 ,

D1 =

1

ξ1

1

ξ12

,

,

C33 = 1, D2 = 0,

Cij = 0,

i ̸= j,

D3 = 0,

and we obtain the following heat equation

−

1 ∂

ξ1 ∂ξ1



ξ1 aε (ξ)

∂ T ε (ξ) ∂ξ1

 −

∂ ξ12 ∂ξ2 1



aε (ξ)

∂ T ε (ξ) ∂ξ2

 −

∂ ∂ξ3



aε (ξ)

∂ T ε (ξ) ∂ξ3



= f (ξ),

which implies that the coefficients aε are periodic in all the radial, circumferential and axial directions. The structure is often used in the real applications, such as the hollow cylinder composed of composite materials, shown in Fig. 2(c) with the representative cell depicted in Fig. 2(b). Due to the expression of C22 , we obtain a real quasi-periodic problem. To further simplify the SOTS expansion, consider the axisymmetric case. The typical structure is illustrated in Fig. 2(e) with periodic axisymmetric plane shown in Fig. 2(d). The normalized cell domain Υ is chosen as in Fig. 1(c). Assume that the coefficient aε , the source term f and boundary flux q are all independent of ξ2 . Then, T ε does not depend on the polar angle. For convenience of our SOTS deduction, we consider the following axisymmetric heat conduction problem

     1 ∂ ∂ T ε (ξ) ∂ ∂ T ε (ξ)  ε ε  − ξ a (ξ) − a (ξ) = f (ξ) in Θ , 1  ξ1 ∂ξ1 ∂ξ1 ∂ξ2 ∂ξ2 (58)  ∂ T0 (ξ)  T0 (ξ) = 0 on Γ¯ 1 , aε (ξ) nξi = q(ξ) on Γ¯ 2 ∂ξi where Θ represents the two-dimensional axisymmetric plane. It should be emphasized that ξ = (ξ1 , ξ2 ) with ξ1 and ξ2 being the radial and axial coordinates, respectively. Γ¯ 1 and Γ¯ 2 are the boundary lines in the axial direction. According to Eq. (21), by some manipulations, the SOTS expansion is given by T ε (ξ) = T0 (ξ) + ε Nα1 (η)

  1 ∂ T0 (ξ) ∂ T0 (ξ) ∂ 2 T0 (ξ) + ε 2 Nα1 α2 (η) + Mα1 (η) + O(ε3 ), ∂ξα2 ∂ξα1 ∂ξα2 ξ1 ∂ξα2

(59)

where Nα1 (η) satisfies the cell problem

   ∂ Nα1 (η) ∂ a(η) ∂   a (η) = in Υ , −  ∂ηi ∂ηi ∂ηα1     Nα1 (η) = 0 Υ

and the homogenized problem for T0 (ξ) becomes

   1 ∂ 0 ∂ T0 (ξ)   − ξ a = f (ξ) in Θ , 1 ij  ξ ∂ξ ∂ξj 1 i ∂ T0 (ξ)   T0 (ξ) = 0 on Γ¯ 1 , a0ij nξi = q(ξ) on Γ¯ 2 . ∂ξj The expressions of the homogenized coefficients a0ij are Eq. (16) formally, but does not depend on the macroscopic variable ξ . Nα1 α2 (η) and Mα1 (η) satisfy the following cell problems, respectively.

   ∂ Nα1 α2 (η) ∂ Nα1 (η) ∂(a(η)Nα1 (η)) ∂   a(η) = a(η) + + a(η)δα1 α2 − a0α1 α2 in Υ , − ∂ηi ∂ηi ∂ηα2 ∂ηα2     Nα1 α2 (η) = 0. Υ    ∂ Mα1 (η) ∂ Nα1 (η) ∂   a(η) = a(η) + a(η)δ1α1 − a01α1 in Υ , − ∂ηi ∂ηi ∂η1     Mα1 (η) = 0. Υ

100

Q. Ma et al. / Journal of Computational and Applied Mathematics 306 (2016) 87–115

(a) Cylindrical coordinate system.

(b) Three-dimensional cell.

(c) Composite hollow cylinder.

(d) Composite axisymmetric plane.

(e) Axisymmetric hollow cylinder.

Fig. 2. Cylindrical coordinate system and the quasi-periodic composite structures.

Q. Ma et al. / Journal of Computational and Applied Mathematics 306 (2016) 87–115

101

Now the SOTS finite element computation can be carried out directly on the axisymmetric plane Θ . If the cells are arranged in irregular manner in the Θ , for example, the quasi-periodically distributed composites laid on the outer shell of the space vehicle because of its curved profile. At this time, the SOTS method cannot be applied directly. We can instead make use of another transformation based on the axisymmetric problem (58), and transform the axisymmetric plane to a new domain in the body-fitted coordinate system, for instance, in which the cells are periodically distributed. Then, the new SOTS expansion is obtained. If the temperature only depends on ξ1 , the heat equation further reduces to

  dT ε ε ξ1 a (ξ1 ) = f (ξ1 ), − ξ1 dξ1 dξ1 1 d

which is the one-dimensional equation described in Section 3.3. Thus, it immediately follows that Theorem 4. Let T ε (ξ1 ) satisfy the following plane axisymmetric heat conduction problem

   dT ε (ξ1 ) 1 d  ε  ξ a (ξ ) = f (ξ1 ) in [a0 , b0 ], −  1 1 ξ1 dξ1 dξ1 dT ε  T ε (a ) = 0;  aε (ξ1 ) (b0 ) = q(ξ1 ). 0 dξ1 Then, T ε (ξ1 ) has the following SOTS expansion T ε (ξ1 ) = T0 (ξ1 ) + ε N1 (η1 )

dT0 (ξ1 ) dξ1

+ ε 2 N11 (η1 )

d2 T0 (ξ1 ) d2 ξ 1

+ O(ε 3 ).

The expressions of the homogenized coefficients a0 , and the cell functions N1 (η1 ) and N11 (η1 ) are Eqs. (44)–(46), respectively, and T0 satisfies the homogenized problem

   1 d  0 dT0 (ξ1 )  ξ a = f (ξ1 ) in [a0 , b0 ], −  1 ξ1 dξ1 dξ1 dT0   T0 (a0 ) = 0; a0 (b0 ) = q(ξ1 ). dξ1 If the regularity assumptions in Theorem 3 are satisfied, then the following error estimation holds

∥T ε (ξ1 ) − T ε,2 (ξ1 , η1 )∥H 1 ([a0 ,b0 ]) ≤ c ε where c is a constant independent of ε . Remark 3. By the regularity supposition of Cij and J, here we do not consider the case ξ1 = 0, i.e. the transformation is non-degenerate. 3.5. Spherical symmetric problem The spherical coordinate system is shown in Fig. 3(a). Also let (ξ1 , ξ2 , ξ3 ) be the conventional coordinate (r , θ , φ). One representative macroscopic structure with periodicity in (ξ1 , ξ2 , ξ3 ) is a composite hollow sphere shown in Fig. 3(b) with the cell domain chosen in Fig. 2(b). Consider the following transformation x1 = ξ1 sin ξ2 cos ξ3 ,

x2 = ξ1 sin ξ2 sin ξ3 ,

x3 = ξ1 cos ξ2 .

We have C11 = 1,

C22 =

J = ξ12 sin ξ2 ,

1

ξ12

,

D1 =

2

ξ1

,

1

, ξ12 sin2 ξ2 cot ξ2 D2 = , ξ12

C33 =

Cij = 0,

i ̸= j,

D3 = 0

and the heat conduction equation becomes

−

∂ ξ12 ∂ξ1 1

      ∂ ∂ T ε (ξ) 1 ∂ T ε (ξ) sec2 ξ2 ∂ ∂ T ε (ξ) ε ξ12 aε (ξ) − 2 sin ξ2 aε (ξ) − a (ξ) = f (ξ). ∂ξ1 ∂ξ2 ∂ξ3 ξ1 sin ξ2 ∂ξ2 ξ12 ∂ξ3

It is also a quasi-periodic problem, and the cell problems and homogenized solutions have the complicated formulations due to the metric values C22 and C33 . Only by reducing to the one-dimensional spherical symmetric problem, the cell functions can be separated from the macroscopic variable, the corresponding equation then becomes

  dT ε (ξ1 ) 2 ε − 2 ξ a (ξ1 ) = f (ξ1 ), dξ1 ξ1 dξ1 1 1

d

which is also the same as the one-dimensional problem (37) formally, and immediately we have

102

Q. Ma et al. / Journal of Computational and Applied Mathematics 306 (2016) 87–115

(a) Spherical coordinate system.

(b) Hollow sphere with quasi-periodic structure.

Fig. 3. Spherical coordinate system and the quasi-periodic composite structure.

Theorem 5. Let T ε (ξ1 ) satisfy the following spherical symmetric heat conduction problem

   dT ε (ξ1 ) 1 d  2 ε  ξ a (ξ1 ) = f (ξ1 ) in [a0 , b0 ], − 2 dξ1 ξ1 dξ1 1  dT ε  T ε (a0 ) = 0; aε (ξ1 ) (b0 ) = q(ξ1 ). dξ1 Then, T ε (ξ1 ) has the following SOTS expansion T ε (ξ1 ) = T0 (ξ1 ) + ε N1 (η1 )

dT0 (ξ1 ) dξ1

+ ε 2 N11 (η1 )

d2 T0 (ξ1 ) d2 ξ1

+ O(ε3 ).

The expressions of the homogenized coefficients a0 , and the cell functions N1 (η1 ) and N11 (η1 ) are Eqs. (44)–(46), respectively, T0 satisfies the homogenized problem

   1 d dT0 (ξ1 )  2 0  ξ a = f (ξ1 ) in [a0 , b0 ], − 2 dξ1 ξ1 dξ1 1  dT0  T0 (a0 ) = 0; (b0 ) = q(ξ1 ). a0 dξ1 If the regularity conditions in Theorem 3 are satisfied, then the following error estimation holds

∥T ε (ξ1 ) − T ε,2 (ξ1 , η1 )∥H 1 ([a0 ,b0 ]) ≤ c ε, where c is a constant independent of ε . Remark 4. We did not consider the degenerate case such as ξ1 = 0, ξ2 = 0 or ξ2 = π , but in one-dimensional problem, the general SOTS expression (55) does not contain 1/ξ1 as the two-dimensional axisymmetric case Eq. (59) as well as the term 1/ sin ξ2 , both of which will cause the singularity. Therefore there is not any problem in the practical finite element computation and assembly of the SOTS approximation solutions for the spherical symmetric problem. 4. SOTS finite element algorithm Based on the second-order asymptotic expansion (21), we present the corresponding finite element algorithm. It is noted that the cell functions and homogenized coefficients are defined point-by-point in the macroscopic domain, so the macroscopic and microscopic finite element meshes are coupled when computing these terms and assembling the final SOTS approximate solutions. Here, we give the computing procedures as follows, 1. Set the material conductivity aε , configure the periodic cell domain Υ , the macroscopic domains Θ and Ω , and make the v space discretization for the three domains. Set the number of nodes and elements for Υ to be NΥv and NΥe , also NΘ and e NΘ for Θ , respectively. Choose a proper coordinate transformation and compute Cij , J and Dj .

Q. Ma et al. / Journal of Computational and Applied Mathematics 306 (2016) 87–115

103

2. For the sake of simplicity, in each degree of freedom in Θ , let the nodal solutions of the cell functions be also expressed by Nα1 (ξ, η), Nα1 α2 (ξ, η) and Mα1 (ξ, η), and they are presented in terms of node values by Nα1 (ξ, η) = HΥe (η)Nαe 1 (ξ, η), Nα1 α2 (ξ, η) = HΥe (η)Nαe 1 α2 (ξ, η), Mα1 (ξ, η) = HΥe (η)Mαe 1 (ξ, η), respectively in each element, where HΥe (η) denotes the shape function of the cell domain, and Nαe 1 (ξ, η), Nαe 1 α2 (ξ, η) and Mαe 1 (ξ, η) are the element cell solutions. By finite element assembling process, the global space discretized systems are given by

 Ne Υ     KΥe (ξ, η)Nαe 1 (ξ, η) = Rαe 1 (ξ, η),      e=e1  NΥ  KΥe (ξ, η)Nαe 1 (ξ, η) = Rαe 1 α2 (ξ, η),   e=1   e  NΥ      KΥe (ξ, η)Mαe 1 (ξ, η) = R¯ αe 1 (ξ, η), 

(60)

e=1

in which

∂ HΥe (η) ∂ HΥe (η) dη, ∂ηi ∂ηj Υe  ∂ HΥe (η) Rαe 1 (ξ, η) = − Ckα1 (ξ)a(η) dη, ∂ηk Υe    ∂ Nαe 1 (ξ, η) + Cα1 α2 (ξ)a(η) − Cα2 k (ξ)a0kα1 (ξ) HΥe (η)dη Rαe 1 α2 (ξ, η) = Cα2 k (ξ)a(η) ∂ηk Υ e ∂ HΥe (η) e dη, − Ckα2 (ξ)a(η)Nα1 (ξ, η) ∂ηk Υe     ∂ a0jα1 (ξ) e ∂ Nαe 1 (ξ, η) 0 e ¯Rαe (ξ, η) = + δ HΥ (η)dη D (ξ) a (η) a (η) − a (ξ) H (η) d η − C (ξ) k α k ij k α Υ 1 1 1 ∂ξi Υe  ∂ηk  Υe  e 2 e ∂ Nα1 (ξ, η) ∂ HΥe (η) ∂ Nα1 (ξ, η) e − HΥ (η) dη − Cij (ξ)a(η) , ∂ξj ∂ηi ∂ξi ∂ηj Υe KΥe (ξ, η) =



Cij (ξ)a(η)

where Υ e is the element of the cell domain, and a0ij (ξ) is computed by Eq. (16). Applying periodic boundary condition, the cell functions are obtained. It is observed that the discretized systems (60) share the same global stiff matrix, which is very convenient and time-saving for obtaining all the cell functions. 3. Compute the homogenized heat conduction problem for T0 (ξ) in Θ . Analogously, define the nodal solution of the homogenized temperature also by T0 (ξ), which is expressed in each element by e T0 (ξ) = HΘ (ξ)T0e (ξ), e where HΘ (ξ) denotes the shape function of the homogenized domain and T0e (ξ) is the element homogenized solution. The global discretized system is given by e

NΘ 

e KΘe (ξ)T0e (ξ) = RΘ (ξ),

e=1

where KΘ (ξ) = e

e RΘ (ξ) =

 Θ

e

Γ¯ 2e

∂ HΘe (ξ) ∂ HΘe (ξ) J (ξ)dξ , ∂ξk ∂ξi  1 e e q(ξ)HΘ (ξ)J (ξ)ds + f (ξ)HΘ (ξ)J (ξ)dξ . E (ξ) Θe

Cij (ξ)a0jk (ξ)

Θ ε and Γ¯ 2e are the elements of the homogeneous domain Θ and boundary Γ¯ 2 , respectively. Applying the boundary condition on Γ¯ 1 , the homogenized solution can be obtained. 4. Assemble the first-order and second-order solutions T ε,1 (ξ, η) and T ε,2 (ξ, η) by the expressions in Eq. (22).

104

Q. Ma et al. / Journal of Computational and Applied Mathematics 306 (2016) 87–115

(b) Cross section domain of the layered material (ε = 1/4).

(a) Layered cell domain.

Fig. 4. Multi-layered thick wall cylinder and the cell domain.

5. Compute the first-order and second-order thermal gradient solutions by

  ∂ Nα1 ∂ T0 ∂ Nα1 ∂ T0 ∂ T ε,1 ∂ T0 ∂ 2 T0 = + +ε + Nα1 , ∂ξk ∂ξk ∂ηk ∂ξα1 ∂ξk ∂ξα1 ∂ξα1 ∂ξk   ∂ Nα1 α2 ∂ 2 T0 ∂ Mα1 ∂ T0 ∂ T ε,1 ∂ T ε,2 = +ε + ∂ξk ∂ξk ∂ηk ∂ξα1 ∂ξα2 ∂ηk ∂ξα1   2 ∂ N ∂ Mα1 ∂ T0 ∂ T ∂ 2 T0 ∂ 3 T0 α1 α2 0 + ε2 + Mα1 . + Nα1 α2 + ∂ξk ∂ξα1 ∂ξα2 ∂ξα1 ∂ξα2 ∂ξk ∂ξk ∂ξα1 ∂ξk ∂ξα1 6. Obtain the approximations T ε,1 (ξ, η) and T ε,2 (ξ, η) in the real domain Ω , and compute the real first-order and secondorder thermal gradient solutions by

∂ T ε,1 (x) ∂ξk ∂ T ε,1 (ξ, η) = , ∂ xi ∂ xi ∂ xi

∂ T ε,2 (x) ∂ξk ∂ T ε,2 (ξ, η) = . ∂ xi ∂ xi ∂ xi

Here, the mapping from Ω and Θ can be implemented by the mesh generation method that is often used in the finite difference computation. The transformation from the finite difference mesh to finite element mesh is also easy, and for instance, we can use the bilinear quadrilateral element to perform the computation on the transformed finite element mesh, but the quality of each element should be carefully investigated. Then, both the Cij and J should be computed point-by-point. Using this technology, we can expand the present asymptotic analysis finite element algorithm to more complex situations, i.e. the multi-scale simulation of the heat transfer problem for the large-area TPS of the spacecraft, which is composed of periodic composite materials with complicated structures and curved profiles. It is also seen that because of the coupling v of these two scales, when ε is smaller, the macroscopic mesh size also becomes smaller with much bigger NΘ , and we have to resolve with the aid of parallel computing to obtain all of the cell functions. Of course, if the two scales are decoupled, like the axisymmetric and spherical symmetric problems, only one group of representative cell functions and the constant homogenized coefficients are sufficient. 5. Numerical examples and discussions In this section, three examples are presented to demonstrate the validity of the proposed asymptotic analysis finiteelement algorithm. Throughout this section, the piecewise linear one-dimensional line or two-dimensional triangle element is used, and all of the material coefficients and the length of the structure are non-dimensional. 5.1. Plane axisymmetric heat transfer problem The first example is the heat conduction problem for the infinite thick wall cylinder with inner and outer radius as 1 and 2, respectively. The two materials are arranged periodically in radial direction r and homogeneous in axial direction z. Fig. 4(a) and (b) show the cell domain and the cross section with ε = 1/4, respectively. Let the inner temperature be

Q. Ma et al. / Journal of Computational and Applied Mathematics 306 (2016) 87–115

(a) N1 .

105

(b) N11 . Fig. 5. Cell functions in the domain Υ . Table 1 Non-dimensional conduction coefficients for each material. Heat conductivity (a) Material 1 (a1 ) Material 2 (a2 )

2 0.05

zero, and the uniform heat gradient 1 is imposed on the outer surface. Now the temperature T ε is the function of only r and satisfies the following boundary value problem

   1 d dT ε   raε = 1, − r dr

dr

dT ε

  T ε (1) = 0,

dr

1 < r < 2, (61)

(2) = 1.

The non-dimensional conductivities for the two materials are listed in Table 1. From the explicit expressions of the cell functions and homogenized coefficients in Proposition 1, we have a0 =

2a1 a2 a1 + a2

= 0.09756.

(62)

The cell solutions N1 and N11 are presented as

 a −a 1 2  η1 , 0 ≤ η1 ≤ 0.5, − a1 + a2 N1 =   a1 − a2 (η1 − 1), 0.5 < η1 ≤ 1. a1 + a2    a1 − a2 1 2 a2   − η1 + η1 , 0 ≤ η1 ≤ 0.5, − a1 + a2 2 2(a1 + a2 )   N11 =  a − a2 1 a1   1 − (η1 − 1)2 − (η1 − 1) , 0.5 < η1 ≤ 1. a1 + a2 2 2(a1 + a2 ) respectively, which are shown in Fig. 5. The homogenized temperature T0 satisfies

   1 d dT0   ra0 = 1, − r dr

  T0 (1) = 0,

dr

dT0 dr

(2) =

1 < x1 < 2, a˜ ε a0

,

(63)

(64)

106

Q. Ma et al. / Journal of Computational and Applied Mathematics 306 (2016) 87–115

(a) ε = 1.

(b) ε = 1/2.

(c) ε = 1/4.

(d) ε = 1/8. Fig. 6. Comparisons of the computed temperature approximations with different ε .

where a˜ ε = 0.05 is the heat conductivity of the material lying on the outer surface. The explicit expression and the first and second derivatives of the temperature are obtained as

  1  2  −r + 8(˜aε + 1) ln r + 1 , T0 =   0  4a      dT0 4(˜aε + 1) 1 = 0 −r + , dr 2a r        d2 T 1 4(˜aε + 1)   0 = − 1 − . 2 0 2 dr

2a

(65)

r

According to the expansion, the asymptotic approximate solutions are assembled by Eqs. (63)–(65). Fig. 6 shows the temperature distributions of T0 , T ε,1 and T ε,2 by different lines types in various ε . Since the exact solutions of Eq. (61) are not easily obtained, we take T ε as the finite-element solution that is solved by 1024 elements in the fine domain shown in Fig. 6 by the circle symbols, and compare with the asymptotic solutions. Because of the discontinuous coefficients in Table 1, the temperature is continuous but not smooth at the interface of the two materials. As ε approaches to 0, it is seen Figs. 6(a) to 6(d) that the homogenized solution T0 and the corrections solutions T ε,1 and T ε,2 reach the exact solution T ε asymptotically, but the local oscillation within the periodicity always exists. The homogenized solution T0 is smooth enough everywhere and by adding the correction terms, T ε,1 and T ε,2 give the better approximation. It is also seen from Fig. 6 that on the inner boundary r = 1, we obtain better results than that on the outer flux boundary r = 2, due to the residual −εEC11 G˜ ε obtained in Eq. (57), and this term will also approach to 0.

Q. Ma et al. / Journal of Computational and Applied Mathematics 306 (2016) 87–115

107

Fig. 7. Comparison of the convergence curves for the first and second order solutions. Table 2 Comparison of the computed relative errors with different ε .

ε

∥e0T ∥L2 ∥T ε ∥L2

∥e1T ∥L2 ∥T ε ∥L2

∥e2T ∥L2 ∥T ε ∥L2

∥e0T ∥H 1 ∥T ε ∥H 1

∥e1T ∥H 1 ∥T ε ∥H 1

∥e2T ∥H 1 ∥T ε ∥H 1

1 1/2 1/4 1/8 1/16 1/32 1/64

2.0564E−0 5.5305E−1 2.2741E−1 1.0402E−1 4.9815E−2 2.4383E−2 1.2063E−2

1.2665E−0 3.3280E−1 1.3586E−1 6.2023E−2 2.9688E−2 1.4530E−2 7.1881E−3

9.0000E−1 2.8532E−1 1.2614E−1 5.9798E−2 2.9154E−2 1.4399E−2 7.1558E−3

1.5200E−1 9.2340E−1 7.8497E−1 7.3287E−1 7.1009E−1 6.9942E−1 6.9426E−1

8.7095E−1 2.8875E−1 1.2751E−1 6.0477E−2 2.9501E−2 1.4575E−2 7.2448E−3

7.8378E−1 2.5176E−1 1.0904E−1 5.1168E−2 2.4822E−2 1.2227E−2 6.0669E−3

To compare these solutions quantitatively, define the error functions e0T = T ε − T0 ,

e1T = T ε − T ε,1 ,

e2T = T ε − T ε,2 .

2

Table 2 lists the relative errors in the L and H 1 (energy) norms for all the three approximate solutions with different periodicities. For bigger ε , the asymptotic solutions cannot give good approximations as the existence of residual terms F˜ ε and EC11 G˜ ε in Eq. (57). However, as ε are becoming smaller, all of the errors are decreasing. In the same ε , T ε,2 is better than the other two solutions. To check the convergence of the present asymptotic analysis finite-element algorithm, Fig. 7 shows the convergence curves of the relative errors for the first and second order solutions, from which we can see that the convergence order for both the solutions is ε . Here we choose the semi-norm as the H 1 norm, and the relative errors are better than the errors in L2 norm due to the simplicity of the problem. The thermal gradient in radial direction is shown in Fig. 8 with ε = 1/8 and 1/16, and there exists a clear jump on the interface of the two materials, which will never disappear as the ε becomes smaller. As is illustrated above, the present second-order asymptotic expansion can also reproduce the oscillating behavior of the solution and it is necessary to add the second-order corrector to simulate the thermal gradient more precisely. 5.2. Two-dimensional shear transformation In this section, we consider a non-orthogonal transformation. The periodic cell domain is shown in Fig. 9(a), the original macroscopic domain is depicted in Fig. 9(b) with ε = 1/16, which is obtained by the following shear transformation x1 = ξ1 + ωξ2 ,

x2 = ξ2 ,

based on the regular domain Θ shown in Fig. 9(c). It follows that C11 = 1 + ω2 ,

C12 = C21 = −ω,

C22 = 1,

J = 1,

D1 = D2 = 0

and the heat conduction problem in Ω considered is formulated in the domain Θ as

         ε ∂ ∂T ε ∂ ∂T ε ∂ ∂T ε −(1 + ω2 ) ∂ aε ∂ T +ω aε + aε − aε =f ∂ξ1 ∂ξ1 ∂ξ1 ∂ξ2 ∂ξ2 ∂ξ1 ∂ξ2 ∂ξ2  ε T = 0 in ∂ Θ .

in Θ ,

Since Cij is constant, the cell functions Nα1 , Nα1 α2 and Mα1 are all decoupled from macroscopic variable ξ . Set ω = 0.25 and the heat source term f = 10, and then the heat conductivities a1 = 1 and a2 = 0.001.

108

Q. Ma et al. / Journal of Computational and Applied Mathematics 306 (2016) 87–115

(a) ε = 1/8.

(b) ε = 1/16. Fig. 8. Comparison of the computed thermal gradients in radial direction r with different ε .

(a) Cell domain Υ .

(b) Macroscopic domain Ω (ε = 1/16).

(c) Macroscopic domain Θ (ε = 1/16).

Fig. 9. Cell domain and two macroscopic domains.

The computational meshes are listed for the cell and homogenized domains as well as the original fine mesh with different

ε in Table 3. As observed that it will take so much time for the FE fine computation when the mesh size is becoming smaller with decreasing ε , while for the SOTS method, the computation time is fixed as its independence of ε . Therefore, it is preferable to use the asymptotic computation for the composite materials involving smaller periodicity.

Q. Ma et al. / Journal of Computational and Applied Mathematics 306 (2016) 87–115

109

Table 3 The information of the computational meshes for the shear transformation. Original fine mesh

Nodes Elements Mesh size

ε = 1/4

ε = 1/8

ε = 1/16

ε = 1/32

2577 4960 3.3146E−2

10 113 19 840 1.5673E−2

40 065 79 360 8.2864E−3

159 489 317 440 4.1432E−3

Cell mesh

Homogenized mesh

1203 2276 5.3961E−2

1089 2048 2.2097E−2

(a) FE fine solution T ε .

(b) Homogenized solution T0 .

(c) First-order approximation T ε,1 .

(d) Second-order approximation T ε,2 .

Fig. 10. Comparisons of the FE solution on fine mesh and the asymptotic solutions.

It is noted that the heat conductivity becomes the second-order tensor a0ij , which is computed by a0ij

0.57427 = 0.044433



0.0428539 . 0.581825



Fig. 10 shows the finite element solution on the refined mesh and the asymptotic solutions with the same range of scales. It is clearly seen that, compared with the fine solution T ε in Fig. 10(a), the first-order approximation T ε,1 shown in Fig. 10(c) is not enough to give the accurate correction of the homogenized solution T0 in Fig. 10(b), whereas there is not any difference between the second-order finite-element solution T ε,2 in Fig. 10(d) and the FE fine solution. Table 4 shows the computed relative errors for the asymptotic solutions with different periodicities and Fig. 11 gives the convergence curves of the asymptotic approximate solutions with respect to ε 1/2 . It can be seen that there is an

110

Q. Ma et al. / Journal of Computational and Applied Mathematics 306 (2016) 87–115

Table 4 Comparisons of computed errors with different ε .

ε

∥e0T ∥L2 ∥T ε ∥L2

∥e1T ∥L2 ∥T ε ∥L2

∥e2T ∥L2 ∥T ε ∥L2

∥e0T ∥H 1 ∥T ε ∥H 1

∥e1T ∥H 1 ∥T ε ∥H 1

∥e2T ∥H 1 ∥T ε ∥H 1

1/4 1/8 1/16 1/32

9.5807E−1 6.5234E−1 2.4342E−1 6.6925E−2

9.5910E−1 6.5270E−1 2.4304E−1 6.6627E−2

1.6607E−2 1.3145E−2 5.9679E−3 2.5106E−3

9.6712E−1 9.0658E−1 9.0605E−1 9.0398E−1

9.0643E−1 9.0626E−1 9.0555E−1 9.0278E−1

7.9774E−2 7.6182E−2 6.4112E−2 5.8982E−2

Fig. 11. The convergence curves of the asymptotic solutions with respect to ε 1/2 .

approximately linear relationship between the H 1 errors and ε 1/2 for T ε,1 and T ε,2 . We could hardly get any accurate information from T ε,1 due to the large H 1 relative errors, while T ε,2 gives out the very fine approximation to T ε quite well. Fig. 12 shows the computed thermal gradients in x1 direction of all the temperature fields. It is noted that different range of scales is utilized for the homogenized solution ∂ T0 /∂ x1 in Fig. 12(b) and the first-order approximate solution ∂ T ε,1 /∂ x1 in Fig. 12(c) compared with the other solutions because of the big difference between ∂ T ε,1 /∂ x1 and ∂ T ε,2 /∂ x1 , which further confirms the first-order solution is incapable of capturing the locally oscillating behavior of the thermal gradient, and the second-order correctors play an important role in reflecting the physical reality. The real thermal gradient within one periodical cell and on the interface of the two materials is correctly obtained by adding Nα1 α2 and Mα1 . 5.3. Polar coordinate transformation In this section we simulate a real quasi-periodic heat conduction problem. Still use cell domain in Fig. 9(a) and consider the cylinder in Section 5.1, but there is periodicities in both radial direction x1 and circumferential direction x2 . Fig. 13(a) and (b) show the cross sections with two different ε , and it is seen that all of the cells in the macroscopic domain are different from the reference cell Υ . Set the temperature T ε = 0 on both the inner and outer surfaces of the cylinder and the heat source term f = 10. The heat conductivities are the same as in Section 5.2. We then consider the temperature distribution in this cylinder and obtain the asymptotic solutions. By the symmetry of the structure, we choose the computational domains shown in Fig. 13(c) and (d). And by the following polar coordinate transformation x1 = ξ1 cos

π

 ξ2 ,

x2 = ξ1 sin

π

 ξ2 ,

2 2 we change the fan-shaped domains to the rectangle domains shown in Fig. 13(e) and (f) with 1 ≤ ξ1 ≤ 2, 0 ≤ ξ2 ≤ ε , and C11 = 1,

C22 =

4

π 2 ξ12

,

C12 = C21 = 0,

J = ξ1 ,

D1 =

1

ξ1

The heat conduction problem becomes now

  1 ∂  ε ε 4 ∂ ε ∂T ε ∂T  − ξ a − a =f  1   ξ1 ∂ξ1 ∂ξ1 ∂ξ2 π 2 ξ12 ∂ξ2  T ε = 0 on ξ1 = 1, 2,   ε ε  aε ∂ T n + 4 aε ∂ T n = 0 on ξ = 0, ε.  ξ1 ξ 2 ∂ξ1 π 2 ξ12 ∂ξ2 2

in [1, 2] × [0, ε],

,

D2 = 0.

Q. Ma et al. / Journal of Computational and Applied Mathematics 306 (2016) 87–115

(a) FE fine solution ∂ T ε /∂ x1 .

(b) Homogenized solution ∂ T0 /∂ x1 .

(c) First-order approximation ∂ T ε,1 /∂ x1 .

(d) Second-order approximation ∂ T ε,2 /∂ x1 .

111

Fig. 12. Comparison of the computed thermal gradients in x1 direction with ε = 1/16.

As is explained in Section 3.4, due to the existence of C22 , all of the cell functions and the homogenized coefficients are dependent on ξ1 , continuously and non-linearly. For instance, the cell functions Nα1 satisfy the following equations,

   ∂ ∂ N1 4   − a − 2 2  ∂η ∂η π ξ1 1 1    ∂ ∂ N 4 2  − a − 2 2 ∂η1 ∂η1 π ξ1

  ∂ ∂ N1 ∂a a = , ∂η2 ∂η2 ∂η1   ∂ ∂ N2 4 ∂a a = 2 2 , ∂η2 ∂η2 π ξ1 ∂η2

and the homogenized problem for T0 (ξ) becomes

      1 ∂ 4 ∂  0 ∂ T0 0 ∂ T0 0 ∂ T0 0 ∂ T0  − ξ1 a11 + a12 − 2 2 a + a22 = f in [1, 2] × [0, ε],    ∂ξ1 ∂ξ2 ∂ξ2 π ξ1 ∂ξ2 21 ∂ξ1  ξ1 ∂ξ1 T0 = 0 on ξ1 = 1, 2,        0 ∂ T0 4  0 ∂ T0 0 ∂ T0 0 ∂ T0  + a12 nξ1 + a + a22 nξ2 = 0 on ξ2 = 0, ε.  a11 ∂ξ1 ∂ξ2 ∂ξ2 π 2 ξ12 21 ∂ξ1 The meshes are also listed for different ε in Table 5 with same cell mesh in Table 3. Note that the homogenized mesh is also changing with ε . Since the cell functions are dependent only on ξ1 , by a regular triangulation on the homogenized domain [1, 2] × [0, ε], it is sufficient to solve only ( 4ε + 1) × 8 cell functions for each ε . Compared with the decoupled cases in Sections 5.1 and 5.2, although much time is spent on computing the cell functions, it is worthwhile as the asymptotic

112

Q. Ma et al. / Journal of Computational and Applied Mathematics 306 (2016) 87–115

(a) Macroscopic domain Ω (ε = 1/5).

(b) Macroscopic domain Ω (ε = 1/10).

(c) Computational domain (ε = 1/5).

(d) Computational domain (ε = 1/10).

(e) Regular computational domain (ε = 1/5).

(f) Regular computational domain (ε = 1/10).

Fig. 13. Two-dimensional domains and computational domains.

Table 5 The information of the computational meshes for the polar coordinate transformation. Original fine mesh

Nodes Elements Mesh size

Homogenized mesh

ε = 1/5

ε = 1/10

ε = 1/20

ε = 1/40

2301 4360 1.5621E−2

4581 8720 7.8103E−3

9141 17 440 3.9051E−3

18 261 34 880 1.9526E−3

( 4ε + 1) × 5 4 ×4×2 ε 5ε × 7.0711E − 2

behavior of the temperature can be obtained and it will be seen that by the correctors, the original solution with oscillating fluctuation can be also correctly approximated. According to the finite element procedure proposed in Section 4, we take the computed results into the fan-shaped domain and extend it to 1/4 circle by periodicity, illustrated in Fig. 14 with ε = 1/10. The finite element solution on the fine mesh is still used to compare with the asymptotic solutions. In contrast to the homogenization T0 in Fig. 14(b) and the first-order solution T ε,1 in Fig. 14(c), there is not any difference between the T ε in Fig. 14(a) and T ε,2 in Fig. 14(d).

Q. Ma et al. / Journal of Computational and Applied Mathematics 306 (2016) 87–115

(a) FE fine solution T ε .

(b) Homogenized solution T0 .

(c) First-order approximation T ε,1 .

(d) Second-order approximation T ε,2 .

113

Fig. 14. The finite element solutions on fine mesh and the asymptotic solutions with ε = 1/10.

Table 6 Comparison of computed errors with different ε .

ε

∥e0T ∥L2 ∥T ε ∥L2

∥e1T ∥L2 ∥T ε ∥L2

∥e2T ∥L2 ∥T ε ∥L2

∥e0T ∥H 1 ∥T ε ∥H 1

∥e1T ∥H 1 ∥T ε ∥H 1

∥e2T ∥H 1 ∥T ε ∥H 1

1/5 1/10 1/20 1/40

7.8706E−1 4.0471E−1 3.9423E−2 1.3740E−2

7.8748E−1 4.0451E−1 3.8870E−2 1.2876E−2

5.8762E−2 3.4741E−2 1.7906E−2 1.2822E−2

9.9720E−1 9.9439E−1 9.8353E−1 9.4514E−1

9.9366E−1 9.8970E−1 9.7437E−1 9.1961E−1

2.3017E−1 2.2843E−1 2.2519E−1 2.1455E−1

Table 6 shows the evaluated relative errors of the asymptotic solutions in L2 and H 1 norms, from which we see that T ε,1 is not satisfactory to give the detailed information. Fig. 15 depicts the convergence curves of T ε,1 and T ε,2 , showing that there is still linear relationship between the errors and ε 1/2 . The errors converge quite slowly in H 1 norm and the error for T ε,2 is about 20%, with much improvement from about 90% for T ε,1 . The errors in computing the second-order correctors may originate from evaluating the partial directives of ∂ Nα1 /∂ηj with respect to the macroscopic coordinate ξi in solving the cell solution Mα1 as it involves calculating the second-order directive of the Nα1 , while it is not very accurate to obtain the term by just using the piecewise linear element in the cell domain. To compare with the FE fine solution T ε , all the asymptotic solutions are mapped into the fine mesh, the errors may also come partly from the interpolation from the homogenized and cell domains into the fine mesh. By the refinement on the cell mesh and homogenized mesh, the errors can be further reduced at the cost of more computational time.

114

Q. Ma et al. / Journal of Computational and Applied Mathematics 306 (2016) 87–115

Fig. 15. The convergence curves of the asymptotic solutions with respect to ε 1/2 .

(a) FE fine solution ∂ T ε /∂ x1 .

(b) Homogenized solution ∂ T0 /∂ x1 .

(c) First-order approximation ∂ T ε,1 /∂ x1 .

(d) Second-order approximation ∂ T ε,2 /∂ x1 .

Fig. 16. Comparison of the thermal gradient in x1 direction with ε = 1/10.

The approximate thermal gradients in x1 direction are also illustrated in Fig. 16 with ε = 1/10 only on the fan-shaped domain. It is noted that the different range of scales is also used as the big difference again exists between ∂ T ε,1 /∂ x1 in Fig. 16(c) and ∂ T ε,2 /∂ x1 in Fig. 16(d). The gradient jump occurs on the interface of the two materials. The local information is shown clearly in the second-order approximate solution ∂ T ε,2 /∂ x1 in Fig. 16(d), while the oscillating behavior is not correctly captured in ∂ T ε,1 /∂ x1 in Fig. 16(c) by comparing the FE fine solution ∂ T ε /∂ x1 in Fig. 16(a). 6. Conclusions and future work In this work, the asymptotic analysis finite-element algorithm has been presented and applied to analyze the heat conduction problem with small periodic configurations in any proper coordinate transformation. The quasi-periodic problem is obtained due to the adjoint metric tensor Cij existing in the heat equation. The finite element algorithms are carried out for three typical coordinate transformations and the convergence order of numerical results agree well with the theoretical prediction.

Q. Ma et al. / Journal of Computational and Applied Mathematics 306 (2016) 87–115

115

It is obviously concluded that if the transformation is linear, all the cell functions are separated from the macroscopic variable, for Cij are all constants. In other cases, if the cell functions and homogenized coefficient are only dependent on a single variable, for instance ξ1 in Section 5.3, then we can compute a series of typical cell functions and homogenized coefficients in advance, with the other ones obtained by interpolation to reduce the computation. Otherwise, there is not a generally effective method to save the computational cost. For the three dimensional problem, it is worthwhile to use parallel computing to obtain cell functions as all of them are independent of each other. Based on the idea of coordinate transformation, the present asymptotic analysis expansion and corresponding finite element computation can also be extended to the linear elastic problems, especially for the axisymmetric and spherical symmetric structures composed of composite materials, and moreover the static and dynamic thermo-elastic problems can also be considered in our future work. Acknowledgments The research is supported by National Key Basic Research and Development Program (2014CB744100), National Natural Science Foundation of China (11325212, 11501140, 11426074, and 91530319), and also supported by China Post-doctoral Science Foundation (2014M562616). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13]

[14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29]

A. Bensoussan, J.L. Lions, G. Papanicolaou, Asymptotic Analysis for Periodic Structures, North-Holland, Amsterdam, 1978. O.A. Oleinik, A.S. Shamaev, G.A. Yosifian, Mathematical Problems in Elasticity and Homogenization, North-Holland, Amsterdam, 1992. D. Cioranescu, P. Donato, An Introduction to Homogenization, Oxford University Press, New York, 1999. J.L. Lions, Some Methods in the Mathematical Analysis of Systems and their Control, Science Press, Beijing, China, 1981, Gordon and Breach, New York. V.V. Jikov, S.M. Kozlov, O.A. Oleinik, Homogenization of Differential Operators and Integral Functionals, Springer Verlag, Berlin-Heidelberg, 1994. D. Cioranescu, J. Saint Jean Paulin, Homogenization in open sets with holes, SIAM J. Numer. Anal. 24 (1997) 1077–1094. E. Sanchez-Palencia, Non-Homogenous Media and Vibration Theory, in: Lecture Notes in Physics, vol. 127, Springer, Berlin, 1980. G. Allaire, Homogenization and two-scale convergence, SIAM J. Math. Anal. 23 (6) (1992) 1482–1518. G. Allaire, K.E. Ganaoui, Homogenization of a conductive and radiative heat transfer problem, Multiscale Model. Simul. 7 (3) (2009) 1148–1170. 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. A. Abdulle, W. E, B. Engquist, E. Vanden-Eijnden, The heterogeneous multiscale method, Acta Numer. 21 (2012) 1–87. 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. J.Z. Cui, Multi-scale Computational Method For Unified Design of Structure, Components and their Materials. Invited Presentation on ‘‘Chinese Conference of Computational Mechanics, CCCM-2001’’, Dec. 5–8, 2001, Proc. On ‘‘Computational Mechanics in Science and Engineering’’, Peking University Press, 2001. J.Z. Cui, L.Q. Cao, Finite element method based on two-scale asymptotic analysis, Math. Numer. Sin. 1 (1998) 89–102. L.Q. Cao, J.Z. Cui, D.C. Zhu, Multiscale asymptotic analysis and numerical simulation for the second order Helmholtz equations with rapidly oscillating coefficients over general domains, SIAM J. Numer. Anal. 40 (2) (2002) 543–577. L.Q. Cao, J.Z. Cui, D.C. Zhu, J.L. Luo, Spectral analysis and numerical simulation for second order elliptic operator with highly oscillating coefficients in perforated domains with a periodic structure, Sci. China Ser. A 45 (12) (2002) 1588–1602. L.Q. Cao, J.Z. Cui, Asymptotic expansions and numerical algorithms of eigenvalues and eigenfunctions of the dirichlet problems for second order elliptic equations in perforated domains, Numer. Math. 96 (2004) 528–581. 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. 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. F. Su, J.Z. Cui, Z. Xu, Q.L. Dong, Multi-scale method for the quasi-periodic structures of composite materials, Appl. Math. Comput. 217 (2011) 5847–5852. 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, Adv. Mater. Res. 629 (2013) 160–164. 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. Q. Ma, J.Z. Cui, Second-order two-scale analysis method for the heat conduction problem with radiation boundary condition in periodical porous domain, Commun. Comput. Phys. 14 (2013) 1027–1057. Z.Q. Yang, J.Z. Cui, Qiang Ma, The second-order two-scale computation for integrated heat transfer problem with conduction, convection and radiation in periodic porous materials, Discrete Contin. Dyn. Syst. Ser. B 19 (3) (2014) 827–848. Z.Q. Yang, J.Z. Cui, Y.F. Nie, Q. Ma, The second-order two-scale method for heat transfer performances of periodic porous materials with interior surface radiation, CMES Comput. Model. Eng. 88 (5) (2012) 419–442. Q.F. Zhang, J.Z. Cui, Existence theory for Rosseland equation and its homogenized equation, Appl. Math. Mech. 33 (12) (2012) 1595–1612. 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, Internat. J. Numer. Methods Engrg. 60 (2004) 241–269. J.J. Wan, Multi-scale analysis method for dynamic coupled thermoelasticity of composite structures (Ph.D. thesis), The Academy of Mathematics and Systems Science, CAS, China, 2007. Z.H. Yang, Y. Chen, Z.Q. Yang, Qiang Ma, Dynamic thermo-mechanical coupled response of random particulate composites: A statistical two-scale method, Chin. Phys. B 23 (7) (2014) 605–616.