A p-adaptive multi-node extended multiscale finite element method for 2D elastostatic analysis of heterogeneous materials

Computational Materials Science 73 (2013) 79–92

Contents lists available at SciVerse ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

A p-adaptive multi-node extended multiscale finite element method for 2D elastostatic analysis of heterogeneous materials H. Liu, H.W. Zhang ⇑ Department of Engineering Mechanics, Faculty of Vehicle Engineering and Mechanics, State Key Laboratory of Structural Analysis for Industrial Equipment, Dalian University of Technology, Dalian 116024, People’s Republic of China

a r t i c l e

i n f o

Article history: Received 30 November 2012 Accepted 19 February 2013

Keywords: Multi-node extended multiscale finite element method Multiscale computation Adaptive algorithm Heterogeneous materials

a b s t r a c t A multi-node extended multiscale finite element method is developed for 2D elastostatic analysis of the computational models with various microstructures. In addition, an adaptive algorithm is proposed for the coarse-scale mesh based on the developed multiscale method. Then, the basic principles of the multi-node extended multiscale finite element method are introduced in detail. To verify the effectiveness of the adaptive algorithm, some representative numerical experiments are carried out. By comparing with the reference solutions, which are obtained by the standard finite element method on the fine-scale mesh, it can be seen that multiscale solutions with high accuracy will be calculated by combining the developed multiscale method and the proposed adaptive algorithm. Finally, a nearly optimal distribution of macroscopic nodes on the fixed coarse-scale mesh can be found by using the proposed adaptive algorithm. Thus, the contradiction between the amount of calculation and computational accuracy, to some extent, can be balanced by using the adaptive multi-node extended multiscale finite element method. Ó 2013 Elsevier B.V. All rights reserved.

1. Introduction Composite materials are commonly applied into the engineering practice due to their various advantages, such as high stiffness, high damping and corrosion resistant. In addition, most composite materials are multiscale in nature [1]. Generally speaking, the scale of the constituents (microscopic scale) is far less than the scale of the structure (macroscopic scale). Therefore, the grids must be refined to make their sizes less than the smallest scale when simulating the mechanical behavior of complex microstructure composite materials by using the direct methods, such as the finite element method and the finite difference method. This will need significant resources in terms of computer memory and CPU time. In this context, various multiscale computational methods which cannot only save the computational resources but also ensure sufficient accuracy are developed in recent decades. Many of the developed multiscale methods are based on the homogenization theory, such as the asymptotic homogenization method [2–9] and the representative volume element (RVE) method [10–15]. These methods are usually limited to the problems with periodic microstructures. E et al. proposed the heterogeneous multiscale method (HMM) [16–18] and provided a general framework for the design of the multiscale methods. For the elliptic problems, this method can be carried out at the macroscopic scale

⇑ Corresponding author. Tel./fax: +86 411 84706249. E-mail address: [email protected] (H.W. Zhang). 0927-0256/$ - see front matter Ó 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.commatsci.2013.02.025

with some undetermined parameters which can be obtained by solving the elliptic problems on the sub domains with special boundary conditions, then the microscopic heterogeneous properties are reflected to the macroscopic scale and the original multiscale problems can be solved at the macroscopic scale. The multiscale eigenelement method was developed by Xing et al. [19,20], which brings the microscopic characteristics within the macroscopic element to the boundary nodal DOF of the macroscopic element by static condensation. Thus this method is efficient and easy to implement and has been used for the static and dynamic analyses [20]. The multiscale finite element method (MsFEM) was developed by Hou et al. [21,22] and has been applied to solve the second order elliptic boundary value problems with high oscillating coefficients. To deal with the multiscale problems in the solid mechanics, Zhang et al. [23] developed a coupling multiscale finite element method (CMsFEM) for consolidation analysis of heterogeneous saturated porous media. Furthermore, an extended multiscale finite element method (EMsFEM) was also proposed by Zhang et al. for the elastic and elasto-plastic analysis of the periodic lattice truss materials [24] and the continuum heterogeneous materials [25]. In EMsFEM, the multiscale base functions are numerically constructed by solving the static equilibrium equations on the unit cell and these functions can capture the heterogeneous properties at the microscopic scale effectively. In addition, the additional coupling terms of the multiscale functions among different directions in multi-dimensional problems are taken into consideration, which can improve the accuracy significantly.

80

H. Liu, H.W. Zhang / Computational Materials Science 73 (2013) 79–92

Γσ

Γ

Homogeneous microstructure

f

Periodic heterogeneous microstructure

Ω Γu Non-periodic heterogeneous microstructure

Fig. 1. Illustration of the structure with heterogeneous microstructures.

On the other hand, in order to reduce the computation resources and obtain the approximate optimal solutions, many adaptive multiscale methods are developed for the heterogeneous media in resent years. Zheng et al. [26] studied the stationary incompressible Navier–Stokes equations by the variational multiscale (VMS) method with h-adaptive technique. Vernerey and Kabiri [27] proposed a concurrent adaptive multiscale method for the elasticity problems with periodic microstructure. Temizer and Wriggers [28] investigated the finite deformation analysis of microscopic heterogeneous structures by an adaptive multiscale resolution strategy based on the homogenization formulations. Lee et al. [29] described a sequential implicit multiscale finite-volume framework for coupled flow and transport with general prolongation and restriction operations for both pressure and saturation where three adaptive prolongation operators for the saturation are used. Larson et al. [30] developed an adaptive variational multiscale methods for the elliptic problems by using the variational multiscale framework together with a systematic technique for approximation of the fine scale part of the solution. He and Ren [31] presented an adaptive multiscale finite element method for solving the unsaturated water flow problems in heterogeneous porous media with many scales. Hajibeygi and Jenny [32] introduced a space–time adaptive iterative multiscale finite volume method for the multiphase flow problems. In this paper, a multi-node extended multiscale (MEM) finite element method is developed based on our previous work for the mechanical analysis of the heterogeneous materials. In addition, a p-version adaptive elastostatic analysis of the heterogeneous material is carried out by the proposed multiscale computational method. In the previous methods (MsFEM and EMsFEM), for 2D elliptic problems, the number of the nodes of each coarse element is only four. Thus one must refine the coarse-scale mesh to deal with the problems with local large strain gradient or stress concentration. This will reduce the efficiency of the multiscale algorithm greatly. To overcome this deficiency, a kind of coarse element with more than four nodes is presented in MEM. In this context, the nodes can be distributed freely according to the practical needs. For example, we can place more nodes in the domain with large strain gradient and fewer nodes in the domain where the strain gradient is small without refining the coarse-scale mesh. Finally, we can obtain sufficiently accurate solutions with approximation optimal distribution of the nodes on the coarse-scale mesh. The heterogeneous structure and its boundary conditions are shown in Fig. 1, where three kinds of microstructures are taken into con-

sideration in this paper, i.e., homogeneous, periodic heterogeneous and non-periodic heterogeneous microstructures. This paper is organized as follows. In the next section (Section 2), the governing equations of the elliptic problems with oscillating coefficients are introduced briefly. In Section 3, the multinode extended multiscale (MEM) finite element method is presented in detail, which includes the construction of the multiscale base functions, calculation of the equivalent macroscopic matrix, downscaling technique, etc. Then the relative error estimations of the multiscale solution between two adjacent adaptive steps are defined in Section 4 and the macroscopic node adaptive strategy of the coarse-scale mesh is presented in Section 5. Some typical numerical experiments are carried out by the adaptive multi-node extended multiscale finite element method in Section 6, and compared with the solutions which are obtained by using the standard finite element method (FEM) on the fine-scale mesh. Finally, some conclusions are given in Section 7. 2. Governing equations The heterogeneous structure occupying a region X and having a boundary C is shown in Fig. 1. The static equilibrium equations of the heterogeneous materials can be described as

@ reij @xej

eeij ¼

þ fi ¼ 0 in X e 1 @uei @uj þ e e 2 @xj @xi

reij ¼ Deijkl eekl in X

ð1Þ ! in X

ð2Þ ð3Þ

where rij and eij are the components of the stress and strain tensors on the microscopic scale, respectively. fi and ui are the vector components of the volume loads and displacements, respectively. The superscript e is the small parameter which signifies the multiscale nature of the problem. The Dirichlet’s and Neumann’s boundary conditions are defined as

 i on Cu ui ¼ u e rij nj ¼ si on Cr

ð4Þ ð5Þ

 i and s i are the values of with Cu [ Cr = C and Cu \ Cr = ;, where u displacements and forces which are imposed on the boundary, respectively. nj are the components of an outward vector normal to the boundary Cr.

81

H. Liu, H.W. Zhang / Computational Materials Science 73 (2013) 79–92

The equations described above are referring to the fine-scale model. For the multiscale problems, the size of the macroscopic scale is usually much larger than the size of the microscopic scale. A meaningful solution can be obtained by using the standard FEM on the fine-scale mesh where the size of the element hf must be smaller than the finest scale e, i.e., hf  e. Even this is not feasible due to that it needs a large number of the computing resources. While for the proposed multiscale method in the following sections, the heterogeneous properties of the microscopic scale can be captured effectively by the multiscale numerical base functions. Thus the size of the mesh for proposed method can be larger than the size of the finest scale, i.e., hc  e. The original problems can be calculated only at the macroscopic scale and this will save the computing resources greatly. 3. Multi-node extended multiscale finite element method The main purpose of EMsFEM is to construct the multiscale base functions numerically by solving the static equilibrium equations on the sub grids in different directions separately. These functions can capture the heterogeneous properties at the microscopic scale effectively and efficiently. Then the macroscopic equivalent matrices of the coarse element can be obtained by virtue of these functions. Finally, the original problems are only solved at the macroscopic scale, which can reduce the computational cost significantly. Consider a region X with the boundary C as shown in Fig. 2, where there are two set of meshes, i.e., coarse-scale mesh (thick blue lines) and fine-scale mesh (thin black lines). The coarse-scale mesh is used for the multiscale solution at the macroscopic scale by the proposed method and the fine-scale mesh is used for the reference solution at the microscopic scale by the standard FEM. In Fig. 2, the sub domains Xo and Xc denote the oversampling element and the coarse element, respectively. As mentioned above, one must refine the coarse-scale mesh for the mechanical analysis of the heterogeneous structure with large microstructure and large strain gradient by using the original fournode coarse element. This will reduce the computational efficiency greatly. In order to describe the deformation more accurately, a kind of coarse element with more than four nodes is proposed as shown in Fig. 3. The construction of multiscale base functions of the multi-node coarse element is described in the following subsection. 3.1. Construction of the multiscale base functions of the multi-node coarse element From the previous work [23–25], one can conclude that the multiscale base functions with oscillating boundary conditions can provide higher accurate solutions. The oversampling technique

3

4

4

10

8

9

7 11 6 1

2 Four-node coarse element

Sub grids

1

2

5 Multi-node coarse element

Fig. 3. Sub grids and the coarse elements.

is a frequently-used approach to obtain the oscillating boundary conditions to construct the multiscale base functions. This technique is effective for the four-node coarse element as shown in the original EMsFEM, but for the multi-node coarse element, it will be invalid. Therefore a new approach is introduced here to obtain the piecewise oscillating boundary conditions to construct the multiscale base functions of the multi-node coarse element. From Fig. 3, One can see that the multi-node coarse element has four boundaries and there are more than or equal to two macroscopic nodes distributed on each boundary. The construction of the multiscale base functions of the multi-node coarse element can be divided into three steps. First of all, the oscillating boundary values can be obtained by applying the oversampling technique to the corresponding four-node coarse element. Then the piecewise oscillating boundary functions of the multi-node coarse element can be evaluated by the linear combination of the oscillating boundary values obtained in the previous step. Finally, the multiscale base functions of the multi-node coarse element can be obtained by solving the static equilibrium equations on the sub grids with the special boundary constraint conditions which are composed of the piecewise oscillating boundary functions selectively. In order to consider the heterogeneities of the materials within the coarse element and to deal with the vector field problems in the solid mechanics, the multiscale base functions are constructed separately in different coordinate directions. For the 2D problems, two forms of the multiscale base functions are constructed for the interpolation of the displacement field. One is used for the displacement interpolation in the x direction and another is used in the y direction. Here, the construction of the multiscale base functions of the multi-node coarse element will be described step by step as follows.

3' 4' 4

3

e

Sub grids

Γ

Ωc Ω

y

Ωc

Ωο

1 Coarse element

x

Unit displacement

2 Ωo

1' Fig. 2. A region with coarse-scale mesh and fine-scale mesh.

3

Fig. 4. Illustration of the oversampling technique.

2'

H. Liu, H.W. Zhang / Computational Materials Science 73 (2013) 79–92

" w¼

#

w10 xx

0

w20 xx

0

w30 xx

0

w40 xx

0

0

w10 yy

0

w20 yy

0

w30 yy

0

w40 yy

ð6Þ

where the dimensions of wj0 xx and wj0 yy are all no  1, where no is the number of the microscopic nodes on the sub grids within the oversampling element. Then the temporary base functions ui of the node i of the target four-node coarse element can be obtained by the linear combination of wj0 , i.e.

uixx ¼

4 X cxij wj0 xx

ð7Þ

j¼1

uiyy ¼

4 X cyij wj0 yy

ð8Þ

j¼1

where cxij and cyij are the coefficients of linear combination determined by the Kronecker delta properties uixx(xj) = dij and uiyy (xj) = dij, respectively, where xj represents the coordinate of the node j (j = 1, 2, 3, 4) of the target coarse element Xc. It can be proved that the temporary base functions ui can satisfy the normalization conP P ditions, i.e., 4i¼1 uixx ¼ 1 and 4i¼1 uiyy ¼ 1. Without considering the additional coupling terms, the temporary base function matrix u of the target four-node coarse element Xc can be given as

"

u¼

u1xx

0

u2xx

0

u3xx

0

u4xx

0

0

u1yy

0

u2yy

0

u3yy

0

u4yy

# ð9Þ

where the dimensions of uixx and uixx are all nc  1 where nc is the number of the microscopic nodes on the sub grids within the target four-node coarse element. After obtaining the temporary base functions of the four-node coarse element, the oscillating boundary values of the four-node coarse element can be given easily. To make the matter simple, we take si(i = 1, 2, 3, 4) to denote the four boundaries of the multi-node coarse element (see Fig. 5). Thus the values of the temporary base functions ui on the boundary sj of the coarse element can be expressed as uixx jsj in the x direction and uiyy jsj in the y direction, respectively, where uixx jsj and uiyy jsj are the oscillating boundary values of the target coarse element. It should be noted that the oversampling technique is only used once for each multi-node coarse element. The size of the oversampling element is 3  3 microstructures in this paper and this has been stud-

s3 m=1

3

8

6 m=1

m=2

7

m=1

m=2

s2

4

m=2

3.1.1. The oscillating boundary values of the target coarse element Consider a larger domain Xo that covers the coarse element Xc as illustrated in Figs. 2 and 4, in which D1234 is the target coarse element and D10 20 30 40 is the oversampling element. Firstly, the 0 temporary base functions wj0 ðj ¼ 10 ; 20 ; 30 ; 40 Þ of the oversampling element are constructed with linear boundary conditions as shown   0 wj0 xx where wj0 xx (or wj0 yy Þ is the tempoin Fig. 4, and wj0 ¼ 0 0 wj yy rary base function of node j0 of the oversampling element in the x (or y) direction. Taking w10 xx as an example (see Fig. 4), unit displacement is applied on node 10 in the positive x direction, and the values vary linearly along the boundaries 10 20 and 10 40 , just as those in the standard bilinear base functions of the rectangular element in the standard FEM. In order to avoid rigid displacements, the displacements of node 30 are set to zero in the x and y directions. Using the boundary conditions mentioned above, the internal displacement field of the oversampling element can be obtained by solving the static equilibrium equations in the oversampling domain Xo by using the standard FEM. Then, the temporary base function w10 xx can be obtained. Similarly, the rest temporary base functions of the other nodes of the oversampling element can also be constructed with the same approach. Finally, the temporary base function matrix w of the oversampling element can be given as

s4

82

1

m=1

5

m=2

2

s1 Fig. 5. An 8-node coarse element.

ied in our previous work [23–25]. For the structure with periodic heterogeneous microstructure, the coarse element is placed in the center of oversampling element. While for the structure with non-periodic microstructure, the location of coarse element in oversampling element is the same with that in Ref. [23]. 3.1.2. The piecewise oscillating boundary values of the multi-node coarse element For the multi-node coarse element, there are more than or equal to two macroscopic nodes on each boundary, and the number of the macroscopic nodes on each boundary can be different (see Fig. 3). Taking mi (i = 1, 2, 3, 4) to denote the number of the macroscopic nodes on the boundary si, then we have 2 6 mi 6 msi , where msi means the number of the microscopic nodes on sub grids located on the boundary si of the multi-node coarse element. For mi = 2, the multi-node coarse element degenerates into the fournode coarse element. To illustrate the construction process of the piecewise oscillating boundary values of the multi-node coarse element, an 8-node coarse element (mi = 3) is taken as an example (see Fig. 5). Generally speaking, the construction processes of piecewise oscillating boundary values are the same in different directions and similar for each boundary of the multi-node coarse element, so only the process of s1 in the x direction is considered here. The oscillating boundary conditions of s1 can be expressed as

‘11x ¼ u1xx js1 ; ‘12x ¼ u2xx js1 ; ‘13x ¼ u3xx js1 ; ‘14x ¼ u4xx js1

ð10Þ

‘jix

where ¼ uixx jsj represents the value of the temporary base function uixx (i = 1, 2, 3, 4) on the boundary sj, where uixx is the temporary base function for the macroscopic node i of the coarse element in the x direction. By the normalization condition, i.e., P4 i¼1 uixx ¼ 1, we have

‘11x þ ‘12x þ ‘13x þ ‘14x ¼ 1

ð11Þ

To make the multiscale base functions of the multi-node coarse element satisfy the normalization conditions, two oscillating functions L11x and L21x are defined as follows

L11x þ L21x ¼ 1 L11x

ð12Þ

where ¼ þ (or ¼ þ Then another function L21x 2 can be obtained by Eq. (12), i.e., L1x ¼ 1  L11x . Similarly, the oscillating functions L1ix and L2ix of the boundary si can be given as ‘11x

‘14x

L11x

‘12x

‘13x ).

L1ix ¼ ‘iix þ ‘ikx ¼ uixx jsi þ ukxx jsi ; ðk ¼ i  1Þ

ð13Þ

L2ix

ð14Þ

¼1

L1ix

83

H. Liu, H.W. Zhang / Computational Materials Science 73 (2013) 79–92

4

3

0.1×10

7

6

y 1

x

5

2

N 4x xx

N 4x yx

0.1×10

N1xxx

N 7x yx

N 7x xx Fig. 6. Illustration of the numerical multiscale base functions.

where L1ix is the summation of the value on the boundary si of the functions uixx and ukxx, where k = i  1, if k is equal to zero, then let k = 4. Then, for the boundary si of the 8-node coarse element, the m piecewise oscillating boundary values hk;ix in the x direction can be obtained by the linear combination of Llix (l = 1, 2), i.e.

m

hk;ix ¼

2 X akl Llix

ð15Þ

l¼1

where m denotes the serial number of the segments which are divided by the macroscopic nodes on the boundary si, and m = 1  2 for the 8-node coarse element (see Fig. 5). The serial number of the segments is sorted in counterclockwise order, i.e., for the boundary s1, m = 1 stands for segment 15 (from the macroscopic node 1 to 5), and m = 2 for segment 52; m = 1 and m = 2 represent segment 48 and segment 81 for the boundary s4, respectively. k m (k = 1, 2) denotes the trend of the functions hk;ix on the segment m m, i.e., k = 1 means that hk;ix varies from 1 to 0 along the segment m m, while hk;ix varies from 0 to 1 along the segment m for k = 2. akl are the coefficients of the linear combination and are determined  m  by the conditions hk;ix xm ¼ dkl , where xm l l represent the endpoint coordinates of the segment m on the boundary si, in which xm 1 and xm 2 are the first and second endpoint coordinates of the segment m, respectively. The piecewise oscillating boundary values of the other boundaries (s2, s3, s4) can also be obtained by using the same approach. Then the multiscale base functions can be calculated by solving the static equilibrium equations on the sub grids with the specified displacement constraints which are composed of the piecewise oscillating boundary values selectively.

3.1.3. The multiscale base functions of the multi-node coarse element m Similarly, hk;iy is defined as the piecewise oscillating functions of the boundary si of the multi-node coarse element in the y direction, m where the denotations m, k and i have the same meaning with hk;ix described in the previous subsection. Taking the 8-node coarse element as an example (see Fig. 5), the construction process of the multiscale base functions will be given as follows. For the 2D vector field problems, the base functions Ni of the macroscopic node i of the multi-node coarse element are   composed of two parts, i.e., Ni ¼ Nxi Nyi , where Nxi denotes the base function matrix used for the displacement interpolation in the x direction and Nyi is the interpolation function matrix used in the y direction. The dimensions of the function matrices Nxi , Nyi and Ni are 2nc  1, 2nc  1 and 2nc  2, respectively, where nc is the number of the microscopic nodes on the sub grids within the coarse element. Since the additional coupling terms are considered here, the displacement function matrix Nxi can be expressed as  x Nixx Nxi ¼ , where Nxixx and Nxiyx are the values of Nxi in the x Nxiyx and y directions, respectively. Nxiyx represents the additional coupling term and means the displacement field in the y direction within the coarse element induced by imposing the specified displacement constraint conditions Ci on the boundaries of the multi-node coarse element. Such boundary conditions can make the multiscale base functions of the multi-node coarse element meet the rigid body motion conditions (translation and rotation). This can also be verified by numerical method. On the other hand, the boundary heterogeneities of the multi-node coarse element can be taken into considered in such boundary conditions due to the use of the piecewise oscillation boundary values. Furthermore, this can improve the calculation accuracy and it can also be seen in our

84

H. Liu, H.W. Zhang / Computational Materials Science 73 (2013) 79–92

Fig. 7. Illustration of p-extension of the coarse element.

previous work [23–25]. Taking the construction process of Nx1 as an example, it can be evaluated by solving the following static equilibrium equations with imposing the specified displacement constraint conditions C1xx in the x direction and C1yy in the y direction on the boundaries of the multi-node coarse element, i.e.

LNx1 ¼ 0;

in Xc

Nx1xx ¼ C1xx ; Nx1yx ¼ C1yy ¼ 0; on @ Xc

ð16Þ

where L is the elasticity operator and satisfies

 Lu ¼ div D : 12 ðru þ ðruÞT Þ , where D is a fourth order stiffness tensor representing material properties; C1xx and C1yy are the values of the specified displacement constraints C1 in the x and y  C1xx . The specified displacedirections, respectively, i.e., C1 ¼ C1yy ment constraints C1xx in the x direction are the combination of 1

2

h1;1x for the segment 15, h2;4x for the segment 81 and zeros for the other segments on the boundaries of the coarse element. On the other hand, the values of the specified displacement constraints C1yy in the y direction are zeros for all the segments on the boundaries of the coarse element, i.e., C1yy = 0. It should be

noted that only the specified displacement constraints are imposed on the boundaries of the multi-node coarse element. Thus only the boundary tractions are taken into account and the body forces are ignored when calculating the base functions. By solving Eq. (16), the multiscale base function matrix  x N1xx x N1 ¼ can be obtained numerically. Similarly, the rest base Nx1yx functions can also be calculated by solving the static equilibrium equations with the specified displacement constraints which are the combination of the piecewise oscillating boundary values selectively. Thus the numerical base functions for the multi-node coarse element can be expressed as N ¼ ½ N1

N2

    Nn  ¼ Nx1

Ny1

Nx2

Ny2

      Nxn

 Nyn ¼

"

Nx1xx

Ny1xy

Nx2xx

Ny2xy

 

Nxnxx

Nynxy

Nx1yx

Ny1yy

Nx2yx

Ny2yy

      Nxnyx

Nynyy

#

ð17Þ

where the dimensions of the function matrices N; Ni ; Nxi and Nxixx are 2nc  2n, 2nc  2, 2nc  1 and nc  1, respectively, where nc is the number of the microscopic nodes on the sub grids within the coarse element and n is the number of the macroscopic nodes of the multinode coarse element; Nxiyx and Nyixy are the additional coupling terms.

Fig. 8. Flow chart of the adaptive algorithm.

85

H. Liu, H.W. Zhang / Computational Materials Science 73 (2013) 79–92

Fig. 12. Comparisons of the von Mises stress on the left boundary.

It is obvious that the numerical base functions Nxixx (or Nyiyy ) constructed above can satisfy the Kronecker delta properties, i.e.

Nxixx ðXj Þ ¼ dij ¼

Fig. 9. Computational model and its boundary conditions.



1 j¼i 0 j–i

ði; j ¼ 1  nÞ

;

ð18Þ

0.2×5

where Xj denotes the coordinate of the macroscopic node j, i.e., Nxixx (or Nyiyy ) is equal to 1 at the macroscopic node i, while equal to 0 at the other macroscopic nodes. Furthermore, these functions can also Pn x and satisfy the normalization conditions, i.e., i¼1 Ni ¼ 1 Pn y N ¼ 1. The numerical multiscale base functions of the 7-node i¼1 i coarse element are shown in Fig. 6. After obtaining the numerical multiscale base functions of the multi-node coarse element, the relationship between the microscopic and macroscopic displacement fields can be expressed as

us ¼ NUc

ð19Þ

where us and Uc are the microscopic and macroscopic displacement fields of the coarse element.

0.2×5

Sub grids (unit cell)

3.2. Equivalent macroscopic matrices and vectors of the coarse element

Coarse-scale mesh

Once the numerical base functions are constructed, the equivalent macroscopic matrices and vectors of the coarse element can be

Fig. 10. Coarse-scale mesh and the sub grids.

1

1

1

1

1

5

5

5

1

1

1

1

1

1

1

3

3

4

4

1

1

1

1

1

1

3

4

3

3

3

1

1

1

1

1

3

3

3

3

3

1

1

1

1

1

3

3

3

3

3

Initial

Final

Fig. 11. Distributions of the macroscopic nodes and the orders of the coarse elements.

86

H. Liu, H.W. Zhang / Computational Materials Science 73 (2013) 79–92

Fig. 13. Comparisons of the von Mises stress: (a) MEM (Initial), (b) MEM (Final), (c) FEM.

Fig. 14. The L-shaped domain and its boundary conditions.

Fig. 16. The final p-distribution of the coarse-scale mesh.

Fig. 15. Coarse-scale mesh and sub grids of the L-shaped domain.

Fig. 17. Comparisons of von Mises stress on the line AB.

87

H. Liu, H.W. Zhang / Computational Materials Science 73 (2013) 79–92

Fig. 18. Comparisons of von Mises stress on the line CD.

derived easily. The static equilibrium equations of the microstructure (unit cell) can be given as

Ks us ¼ Fs

ð20Þ

where Ks, obtained by using the standard FEM, is the overall stiffness matrices of the microstructure; us and Fs are the microscopic nodal displacement vector and external force vector, respectively. Substituting Eq. (19) into Eq. (20) and then multiplying on the left by the vector NT on both sides of Eq. (20) lead to

Kc Uc ¼ Fc

ð21Þ

where Kc is the equivalent macroscopic stiffness matrix of the multi-node coarse element; Fc is the equivalent macroscopic external force vector of the coarse element. Kc and Fc can be expressed as

Kc ¼ NT Ks N T

Fc ¼ N Fs

ð22Þ ð23Þ

The dimensions of Ks, Kc, Fs and Fc are 2ns  2ns, 2nc  2nc, 2ns  1 and 2nc  1, respectively, where nc and ns are the numbers of the macroscopic nodes of the multi-node coarse element and the microscopic nodes on the sub grids within the coarse element, respectively. Generally speaking, it has ns  nc. Once the equivalent matrices and vectors are obtained, the process of the standard FEM can be used to assemble the global stiffness matrix and load vector on the macroscopic scale, respectively.

Fig. 20. A deep beam with coarse-scale mesh and its boundary conditions.

3.3. Macroscopic computation and downscaling technique In the previous sections, we have discussed the construction of the multiscale numerical base functions and the equivalent macro matrices in detail. Then the original problems can be solved on the macroscopic scale which will reduce the computational resources significantly. The global stiffness matrix of the coarse-scale mesh can be assembled as follows, i.e. i C KC ¼ AM i¼1 Kc

ð24Þ

C where AM i¼1 is the matrix assembly operator, where MC is the number of the coarse elements; Kic , obtained by Eq. (22), is the stiffness matrix of the coarse element i. On the other hand, the external force vector on the macroscopic scale can also be calculated by

i C FC ¼ AM i¼1 Fc

ð25Þ

where Fci is the load vector of the coarse element i. According to the content described above, one can analyze the multiscale mechanical behaviors of the heterogeneous structure on the macroscopic scale. The degrees of freedom are reduced greatly, and then the macroscopic displacement vector UC can be obtained by solving the static equilibrium equations on the macroscopic scale, i.e.

Fig. 19. Comparisons of von Mises stress: (a) MEM (Initial), (b) MEM (Final), (c) FEM.

88

H. Liu, H.W. Zhang / Computational Materials Science 73 (2013) 79–92

Fig. 21. Two kinds of microstructure: (a) homogeneous microstructure, (b) heterogeneous microstructure with a hole.

KC UC ¼ FC

ð26Þ

Meanwhile, the displacements on the microscopic scale can also be calculated conveniently by taking advantage of the mutiscale numerical base functions, i.e., Eq. (19). Thus, the microscopic information such as the microscopic strain and stress can be obtained naturally. 4. Error estimations In order to avoid confusion, we take ~ to denote the multiscale solution (obtained by the proposed method and downscaling technique). The strain energy over the single fine element e (see Fig. 2) can be determined by

1 Ue ¼ 2

Z

~Te

e  D  ~ee dX

Xe

ð27Þ

where ee is the strain of the fine element e. Then, the strain energies over the single coarse element and the whole structure can be expressed as N se X UC ¼ Ue MC X U iC i¼1

qC ¼

UC VC

ð30Þ

where VC represents the volume of the coarse element and qC is the average strain energy density of the coarse element. Then, the reference error of the strain energies for the whole structure between two adjacent adaptive steps can be defined as

U n  U nþ1 X X eU ðXÞ ¼ U nþ1 X

ð31Þ

where n denotes the adaptive step. Similarly, the reference error of average strain energy density of the coarse element can be given as

qn  qnþ1 C C eq ðCÞ ¼ qnþ1 C

ð32Þ

ð28Þ

On the other hand, the global displacement error in L2 norm between two adjacent adaptive steps can be defined as

ð29Þ

n

u

¼ ~ nþ1 ~X  u X

e¼1

UX ¼

where Nse and MC are the number of the fine elements within the coarse element and the number of the coarse elements, respectively. U iC is the strain energy over the coarse element i. Naturally, the average strain energy density of the coarse element can be written as

ffi rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi XNnod h 2  n 2 i ~ nþ1 ~ ni  u þ v~ i  v~ nþ1 u i i i¼1

Fig. 22. The final p-distributions of the coarse-scale mesh: (a) the model with homogeneous microstructure, (b) the model with heterogeneous microstructure.

ð33Þ

H. Liu, H.W. Zhang / Computational Materials Science 73 (2013) 79–92

Fig. 23. Comparisons of von Mises stress of the model with homogeneous microstructure: (a) MEM (Initial), (b) MEM (Final), (c) FEM.

Fig. 24. Comparisons of von Mises stress of the model with heterogeneous microstructure: (a) MEM (Initial), (b) MEM (Final), (c) FEM.

89

90

H. Liu, H.W. Zhang / Computational Materials Science 73 (2013) 79–92

coarse element Xm c and pj be the order of another coarse element that share Cjm . If pm is greater than or equal to pj, the macroscopic node should be increased by one on the edge Cjm . After increasing the macroscopic nodes on the edge Cjm , we need to re-distribute the macroscopic nodes on the edge Cjm uniformly. The p-extension of the coarse element is illustrated in Fig. 7 where CE denotes the coarse element number and p is the order of the coarse element. After increasing p-level of the coarse element Xm c , the multiscale numerical base functions of the coarse elements, whose numbers of the macroscopic nodes were changed, should be re-calculated. Then, another multiscale solution can be calculated by using the new p-distribution. This procedure will be continued until the predetermined convergence conditions (eU(X) < dU and eu(X) < du, where dU and du are the parameters to control the computational accuracy) are reached. The flow chart of the adaptive algorithm for the coarse-scale mesh is depicted in Fig. 8.

Fig. 25. Material distributions.

where Nnod is the number of the microscopic nodes in the fine-scale mesh; ui and vi are the displacements in the x and y directions of the microscopic node i, respectively. The global displacements in L2 norm of the multiscale solutions for the adaptive step n can be given as

n

u ~X ¼

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi XNnod h 2  2 i ~ ni þ v~ ni u i¼1

ð34Þ

respectively. Likewise, we define the reference error of the global displacement field between the adaptive step n and n + 1 as

eu ðXÞ ¼

n

u ~ nþ1 ~X  u

nþ1X

u ~

ð35Þ

X

5. Adaptive strategy for the macroscopic nodes As mentioned above, the number of the macroscopic nodes on the edge si of the multi-node coarse element is denoted by mi and let mk be the smallest among mi(i = 1, 2, 3, 4). The order p of the multi-node coarse element can be defined as p = mk  1. In this paper, the initial orders of each coarse element in the coarse-scale mesh are all p = 1 uniformly. Our p-extension strategy over the coarse-scale mesh is as follows. Firstly, we obtain an initial multiscale solution using low-p-level approximation. Next, we calculate the average strain energy density qC for each coarse element. Let Xm c be the coarse i element where qm C is the largest among qC , i = 1, 2, . . . , MC, where MC is the number of the coarse elements in the coarse-scale mesh. If this coarse element has reached the highest order, select the next coarse element with largest strain energy except this one. Then, the order of the coarse element Xm c will be elevated until the reference error of the average strain energy density eq(C) is lass than a(qC)  dC, where dC is a small error parameter and a(qC) is a correction parameter which is associated with the strain energy of the coarse element. a(qC) can be determined by

aðqC Þ ¼

minðqC Þ  qC þ1 maxðqC Þ

6. Numerical experiments In this section, some demonstrative examples are investigated to illustrate the performance of the adaptive algorithm. It should be mentioned that all of the parameters are dimensionless. The multiscale solution is calculated by using the proposed method and labeled by ‘MEM’. ‘FEM’ denotes the solution which is obtained by the standard FEM on the fine-scale mesh. Example 1. The computational model and its boundary conditions are shown in Fig. 9 where the uniform distributed load q is 105. The corresponding coarse-scale mesh and unit cell are shown in Fig. 10 where the size of the unit cell is 1  1 and each unit cell includes 5  5 fine elements. The material of the unit cell is homogeneous. The Young’s modulus and the Poisson’s ratio are 105 and 0.3, respectively. The tolerance parameters of the whole domain are dU = 5  105 and du = 5  105, respectively. In addition, the strain energy reference error parameter of the coarse element is dC = 1  103. A nearly optimal sequence of p-distribution is obtained and shown in Fig. 11, starting from an initial p = 1 uniform p-distribution. The comparisons of von Mises stress are depicted in Figs. 12 and 13, from which it can be found that the optimized multiscale solution is almost the same with the reference solution. After the

ð36Þ

where max (qC) and min (qC) represent the maximum and minimum strain energies among all of the coarse elements, respectively. The macroscopic node adaptive process is as follows. Let j Cjm ðj ¼ 1; 2; 3; 4Þ be the edges of the coarse element Xm c . If Cm is an edge on the boundary of the coarse-scale mesh, we increase the macroscopic node by one on the edge Cjm . If Cjm is an inter-element edge, we should compare the orders of both the coarse elements that share Cjm to determine whether need to increase the macroscopic nodes on the edge Cjm . Let pm be the order of the

Fig. 26. The final p-distribution of the coarse-scale mesh for the model with nonperiodic heterogeneous microstructures.

H. Liu, H.W. Zhang / Computational Materials Science 73 (2013) 79–92

91

Fig. 27. Displacements in the y direction: (a) MEM (Initial), (b) MEM (Final), (c) FEM.

optimization, the NDF (number of degrees of freedom) of the coarse-scale mesh is only 25.7% of that of the fine-scale mesh. This is to say that about 74.3% NDF is saved in this example. Example 2. A L-shaped domain with homogeneous microstructure is investigated in this example. The computational model and its boundary conditions are shown in Fig. 14 where the uniform distributed load q is 105/12. The coarse-scale mesh and sub grids (unit cell) are shown in Fig. 15 where the size of the unit cell is 1  1 and each unit cell includes 10  10 fine elements. The material properties are the same with those in Example 1. The parameters dU and du are all 1  106, respectively. In addition, the strain energy reference error parameter of the coarse element is dC = 1  103. The final p-distribution of the coarse-scale mesh is shown in Fig. 16 where the orders of the coarse elements in the domain with local stress concentration are significantly highest (p = 8). Thus an approximate optimal distribution of the macroscopic nodes is reached finally. In addition, the von Mises stresses on the lines AB and CD are compared in Figs. 17 and 18, from which one can see that the multiscale solution, after optimization, is in good agreement with the reference solution, especially in the area of stress concentration. The comparisons of the von Mises stress in the whole L-shaped domain are also shown in Fig. 19, where the same conclusion can be drew that the multiscale solution after optimization is very close to the reference solution. On the other hand, NDF of the coarse-scale mesh after optimization is only 7.4% of that of the fine-scale mesh, i.e., about 92.6% NDF is saved here.

Example 3. A cantilever with two kind of different microstructures is studied here. The computational model with coarse-scale mesh is shown in Fig. 20. The left boundary of the model is fixed and a uniform distributed load q with the magnitude of 104 is imposed on the right boundary. The homogeneous and periodic heterogeneous microstructures are taken into consideration as shown in Fig. 21. The Young’s modulus and the Poisson’s ratio are 105 and 0.3, respectively. The parameters dU and du are all 1  106 for the two models, respectively. In addition, the strain energy reference error parameters of the coarse element for these two models are also dC = 1  103. The final distributions of the orders of the coarse elements for the models with two different microstructures are shown in Fig. 22. Furthermore, the comparisons of von Mises stresses are depicted in Figs. 23 and 24, from which we can see that the multiscale solution which is optimized by the adaptive algorithm is very close to the reference solution. Similarly, the von Mises stresses of the unit cell A (see Fig. 20) are also compared with the reference solution in Figs. 23 and 24.

Example 4. The computational model with completely heterogeneous materials is considered in this numerical experiment. This model is same with that in Fig. 9 except the Young’s modulus E and the size L. Two kinds of the materials are considered in this example and their random distributions are shown in Fig. 25, where the Young’s moduli of Phase 1 and Phase 2 are 107 and

Fig. 28. Comparisons of von Mises stresses: (a) MEM (Initial), (b) MEM (Final), (c) FEM.

92

H. Liu, H.W. Zhang / Computational Materials Science 73 (2013) 79–92

106, respectively. In addition, the size of the model is 10  10, i.e., L = 10. There are 10  10 coarse elements in the coarse-scale mesh for the model in this example, and there are 10  10 fine elements which are included within each coarse element (see Fig. 25). The parameters are dU = 1  106 and du = 1  105, respectively. In addition, dC = 1  103 is used in this example. The final p-distribution of the coarse-scale mesh for the computational model with non-periodic heterogeneous microstructures is depicted in Fig. 26. In addition, the contours of the displacements in the y direction and the von Mises stresses are shown in Figs. 27 and 28, from which the same conclusion can be drew that the multiscale solutions after optimization are almost exactly the same with the reference solutions.

7. Conclusions The problem of finding a nearly optimal distribution of the macroscopic nodes on the fixed coarse-scale mesh is discussed by combining with the multi-node extended multiscale (MEM) finite element method and the adaptive algorithm. Firstly, the basic principles of MEM are introduced in this paper which includes the construction of the numerical base functions, the calculation of the equilibrium macroscopic matrices, downscaling technique, and so on. Then, some local and global error indictors are defined to estimate the errors of the multiscale solutions. Finally, the strategy of p-adaptive for the coarse-scale mesh is also formulated and some typical numerical experiments are carried out to verify the effectiveness of the adaptive algorithm. Three different kinds of the microstructures are taken into consideration in these examples, i.e., homogeneous, periodic heterogeneous and non-periodic heterogeneous microstructures. These numerical experiments show that the adaptive algorithm is effective and an approximate optimal distribution of the macroscopic nodes on the fixed coarse-scale mesh can be found by this algorithm. Thus, the contradiction between the amount of calculation and the computational accuracy can be balanced by the proposed adaptive multi-node extended multiscale finite element method effectively. Acknowledgements The supports of the National Natural Science Foundation (11072051, 10721062, 90715037, 10728205, 11232003, 51021140004), the 111 Project (B08014) and the National Key Basic Research Special Foundation of China (2010CB832704) are gratefully acknowledged.

References [1] P. Kanouté et al., Archives of Computational Methods in Engineering 16 (1) (2009) 31–75. [2] J. Pinho-da-Cruz, J.A. Oliveira, F. Teixeira-Dias, Computational Materials Science 45 (4) (2009) 1073–1080. [3] J.A. Oliveira, J. Pinho-da-Cruz, F. Teixeira-Dias, Computational Materials Science 45 (4) (2009) 1081–1096. [4] J.A. Otero et al., Mechanics of Materials 37 (1) (2005) 33–44. [5] K.S. Challagulla, A.V. Georgiades, A.L. Kalamkarov, Composite Structures 79 (3) (2007) 432–444. [6] J. Clayton, P. Chung, Journal of the Mechanics and Physics of Solids 54 (8) (2006) 1604–1639. [7] H. Berger et al., Computational Mechanics 33 (1) (2003) 61–67. [8] H.W. Zhang et al., International Journal for Numerical Methods in Engineering 69 (1) (2007) 87–113. [9] C.M. Chen, N. Kikuchi, F. Rostam-Abadi, Computers & Structures 82 (4–5) (2004) 373–382. [10] Z. Xia, Y. Zhang, F. Ellyin, International Journal of Solids and Structures 40 (8) (2003) 1907–1921. [11] T. Kanit et al., International Journal of Solids and Structures 40 (13–14) (2003) 3647–3679. [12] Z.H. Shan, A.M. Gokhale, Computational Materials Science 24 (2002) 361–379. [13] V. Pensée, Q.C. He, International Journal of Solids and Structures 44 (7-8) (2007) 2225–2243. [14] M. Galli, J. Cugnoni, J. Botsis, European Journal of Mechanics – A/Solids 33 (2012) 31–38. [15] M. Dondero et al., International Journal of Heat and Mass Transfer 54 (17–18) (2011) 3874–3881. [16] W. E, Communications in Mathematical Sciences 1(3) (2003) 423–436. [17] W. E, P.B. Ming, P.W. Zhang, Analysis of the heterogeneous multiscale method for elliptic homogenization problems. Journal of the American Mathematical Society, 2004, 18(1): 121-156. [18] W. E, B. Engquist, Z. Huang, Heterogeneous multiscale method: A general methodology for multiscale modeling. Physical Review B, 2003, 67(9), http:// dx.doi.org/10.1103/PhysRevB.67.092101. [19] Y.F. Xing, Y. Yang, X.M. Wang, Composite Structures 92 (9) (2010) 2265–2275. [20] Y.F. Xing, Y. Yang, Composite Structures 93 (2) (2011) 502–512. [21] T.Y. Hou, X.H. Wu, Z.Q. Cai, Mathematics of Computation 68 (227) (1999) 913– 943. [22] T.Y. Hou, X.H. Wu, Journal of Computational Physics 134 (169–189) (1997). [23] H.W. Zhang, Z.D. Fu, J.K. Wu, Advances in Water Resources 32 (2) (2009) 268– 279. [24] H.W. Zhang, J.K. Wu, Z.D. Fu, Computational Mechanics 45 (6) (2010) 623–635. [25] H.W. Zhang, J.K. Wu, J. Lv, Computational Mechanics 49 (2) (2011) 149–169. [26] H. Zheng, Y. Hou, F. Shi, Journal of Computational Physics 229 (19) (2010) 7030–7041. [27] F.J. Vernerey, M. Kabiri, Computer Methods in Applied Mechanics and Engineering (2012), http://dx.doi.org/10.1016/j.cma.2012.04.021. [28] I_ Temizer, P. Wriggers, Computer Methods in Applied Mechanics and Engineering 200 (37–40) (2011) 2639–2661. [29] S.H. Lee, H. Zhou, H.A. Tchelepi, Journal of Computational Physics 228 (24) (2009) 9036–9058. [30] M.G. Larson, A. Målqvist, Computer Methods in Applied Mechanics and Engineering 196 (21–24) (2007) 2313–2324. [31] X. He, L. Ren, Journal of Hydrology 374 (1–2) (2009) 56–70. [32] H. Hajibeygi, P. Jenny, Adaptive iterative multiscale finite volume method, Journal of Computational Physics 230 (3) (2011) 628–643.