Abstract

A simple multiscale analysis framework for heterogeneous solids based on a computational homogenization technique is presented. The macroscopic strain is linked kinematically to the boundary displacement of a circular or spherical representative volume which contains the microscopic information of the material. The macroscopic stress is obtained from the energy principle between the macroscopic scale and the microscopic scale. This new method is applied to several standard examples to show its accuracy and consistency of the method proposed.

1. Introduction

In the computer modeling of heterogeneous materials, the length scale of interest of the system should be taken into consideration. Though, taking the microstructure into account will often introduce excessive computation cost. To avoid this, through homogenization techniques, the multiscale models were invented to model the microscopic scale behavior with reasonable accuracy and cost.

Multiscale methods can be categorized into hierarchical, semiconcurrent and concurrent methods [1], Figure 1. In hierarchical multiscale methods, information is passed from the fine scale to the coarse scale but not vice versa. Computational homogenization [2] is a classical upscaling technique. Hierarchical multiscale approaches are very efficient.

The basic idea of semiconcurrent multiscale methods is illustrated in Figure 1(c). In semiconcurrent multiscale methods, information is passed from the fine scale to the coarse scale and vice versa. Semiconcurrent methods are of the same computational efficiency as concurrent multiscale methods [3, 4]. Numerous concurrent multiscale methods [5, 6] have been developed that can be classified into “Interface” coupling methods and “Handshake” coupling methods. Interface coupling methods seem to be less effective for dynamic applications as avoiding spurious wave reflections at the “artificial” interface seem to be more problematic.

Early studies [79] in this direction were focused on finding the effective material properties and their bounds. Nevertheless, these models are restricted to limited class of constitutive models and shapes of the constituents. Guedes and Kikuchi [10] were among the first researchers to use the computational homogenization techniques which exploited the flexibility of numerical methods. Application and the solution techniques to inelastic regime and highly nonlinear cases were studied in [11].

A very interesting class of the semiconcurrent multiscale methods is referred to as FE2 [13, 14]. In this multilevel approach, a boundary value problem is solved at the microlevel and this provides the response at the higher level. The low-level boundary value problem is called a representative volume element (RVE) [2]; see Figure 2. Using this approach, one can predict the mechanical behavior at the macrolevel without needing to analytically identify the constitutive response a priori. Therefore, the overall material behavior including nonlinearities arises from the RVE.

This method has proven to be effective in several areas such as structural stress and damage analysis [1520]. For nonlinear cases of RVE, in which the computational cost is significant, Yuan and Fish [21] and Somer et al. [22] offered some solutions to improve the computational efficiency. Some practical examples of exploiting this framework can be found in [2328]. The discretization of the representative volume can be made by the well-known meshfree or finite element methods [29].

Recently, one of the hot areas of research is multiscale modeling of fracture where the FE2 can be applied; see [3035] for example. However it appears that there are some problems with using the conventional rectangular cell as the RVE. The first problem is the questionable ability of the rectangular RVE to model the edge effect as noted in [14]. This means the method will not be accurate for cracks. Moreover, as it is extensively discussed in [30], to treat the response of cracks and shear bands, neither the unit cell or the geometry can be assumed to be periodic. In other words, when a crack touches a boundary, the displacement jump violates the periodic boundary conditions (PBC) on that boundary. Furthermore, applying PBC especially in three dimensions is too cumbersome. This is because when the nodes of the opposite sides of the RVE are not aligned (and this is often the case in unstructured meshes), one has to use projection methods to enforce the boundary conditions.

Therefore, in [30, 36, 37] a virtual unit cell was introduced to provide the energetic consistency of the effective discontinuity. Therefore, in this paper, we propose to use a circular (or spherical) RVE where the macroscopic scale strain is constrained to the displacement of the circular boundary of the representative volume kinematically. This approach is similar to the kinematically constrained microplane model. In that model, the strain tensor is decomposed to the corresponding strain vector to obtain the corresponding stress vector [3841], and so forth. The constitutive relation is given between the strain vector and the stress vector. The new RVE can be later used in crack multiscale propagation problems. However, in this paper we only focus on the formulation of the new method and its applicability to various problems other than the one-to-one comparison to the conventional rectangular RVE. The latter will then be the subject of another paper.

The paper is arranged as follows: the governing equations are given in Section 2. The scale transition scheme of the method proposed in this paper is given in Section 3. Some numerical examples are discussed in Section 4. Finally, we draw conclusions of this paper in Section 5.

2. Governing Equations

The class of multiscale solid mechanics problems of concern in this paper is characterized by the well-known equilibrium boundary value problem at the macroscopic scale where the constitutive relation at each point is defined by the response of a representative volume as shown in Figure 2. The representative volume itself is modeled as a conventional continuum with a simple constitutive law which can be defined by the response of the lower scales if more scales are needed for the accuracy of the solution. When the microscopic structural length scale is comparable with the macroscopic length scale, second-order homogenization may be used, for example, [42].

Consider a bounded domain Ω𝑅𝑑, where 𝑑 is the dimension of the physical space and here in this paper 𝑑=2 unless it is stated otherwise. The extension of the proposed method to three dimensions is straightforward. We assume that the material is heterogeneous in the macroscopic scale and its microstructure includes inclusions and voids (Figure 2). The structure is subjected to a set of displacement and traction boundary conditions on the disjoint complementary parts of the boundary, that is, Γ𝑢Γ𝑡=. Here we are looking for a solution to the boundary value problem at the macroscopic and microscopic scales given by𝜎𝑖𝑗,𝑗+𝑏𝑖𝑢=0,𝐱Ω,𝑖=𝑢𝑖,𝐱Γ𝑢,𝜎𝑖𝑗,𝑛𝑗=𝑡𝑖,𝐱Γ𝑡,(2.1) where 𝜎𝑖𝑗, 𝑏𝑖, 𝑢𝑖, 𝑛𝑗, and 𝑡𝑖 are the Cauchy stress, the body force, the displacement, the unit normal, and the traction measured at point 𝐱, respectively, and the superimposed bar denotes the imposed boundary conditions.

Equivalent to the strong form of the governing equations, the equilibrium condition can be written as the following variational form: Find 𝑢𝑖𝑖, for all 𝛿𝑢𝑖0 such thatΩ𝜎𝑖𝑗𝛿𝜀𝑖𝑗𝑑ΩΩ𝑏𝑖𝛿𝑢𝑖𝑑ΩΓ𝑡𝑡𝑖𝛿𝑢𝑖𝑑Γ=0,(2.2) with𝑖=𝑢𝑖𝑢𝑖1(Ω),𝑢𝑖=𝑢𝑖onΓ𝑢,𝛿𝑖=𝛿𝑢𝑖𝛿𝑢𝑖1(Ω),𝛿𝑢𝑖=0onΓ𝑢,(2.3) in which 𝛿𝜀𝑖𝑗 is the variation of the strain tensor associated with the test function 𝛿𝑢𝑖.

3. Scale Transition Using Circular Substructure

3.1. The Calculation Procedure for the Macroscopic Stress from the Macroscopic Strain

In the context of multilevel finite element analysis, we extract the macroscopic stresses directly by solving a local finite element problem for the circular substructure Ω𝑠 shown in Figure 2. In the following, the script 𝑠 denotes the corresponding components in the microscopic level. The strain tensor measured at the macroscopic scale is transferred to the displacement of the boundary Γ𝑠 of the circular substructure Ω𝑠. The displacement 𝑢𝑠𝑖 at a point on Γ𝑠 is given by𝑢𝑠𝑖(𝜃)=𝑟0𝜀𝑖𝑗𝑛𝑗𝑑𝑟𝑟𝜀𝑖𝑗𝑛𝑗,(3.1) where 𝑟 is the radius of the substructure and 𝜀𝑖𝑗 is the strain tensor at the macroscopic scale as shown in Figure 2. Please note that 𝑛𝑗 in (3.1) is a function of 𝜃. The macroscopic scale and the subscale are then kinematically constrained to each other.

The stress tensor at the macroscopic scale is constructed from the response of the subscale. The simplest choice would be a weighted average of the stress at the subscale, that is,𝜎𝑖𝑗=Ω𝑠𝑤𝜎𝑠𝑖𝑗𝑑Ω,(3.2) where 𝜎𝑖𝑗 is the stress tensor at the macroscopic scale, 𝑤 is a suitably chosen weight function such that Γ𝑠𝑤dΓ=1, and 𝜎𝑠𝑖𝑗 is the subscale stress tensor.

Here, we propose another method based on the energy principle. According to the principle of virtual work, the virtual work done by the stress tensor in the representative volume should be equal to the sum of the virtual work done by the traction on the boundary of the substructure structure:𝑉𝜎𝑖𝑗𝛿𝜀𝑖𝑗=Γ𝑠𝑡𝑠𝑖𝛿𝑢𝑠𝑖𝑑Γ.(3.3) Here, 𝑉 is the volume of the substructure, 𝑡𝑠𝑖 is the traction on Γ𝑠 and 𝑢𝑖 is the displacement on Γ𝑠. Since the displacement is given by (3.1), that is, 𝛿𝑢𝑠𝑖=𝑟𝑛𝑗𝛿𝜀𝑖𝑗, the macroscopic stress tensor can be expressed as𝜎𝑖𝑗=𝑟𝑉Γ𝑠𝑡𝑠𝑖𝑛𝑠𝑗𝑑Γ.(3.4) The boundary of the substructure may be approximated by using a shape function. If the shape function satisfies partition of unity, the integrand of (3.4) can be rewritten as follows:𝜎𝑖𝑗=𝑟𝑉Γ𝑠𝐼Ψ𝐼𝑡𝑠𝑖𝑛𝑠𝑗=𝑟𝑑Γ,(3.5)𝑉𝐼Γ𝑠𝐼Ψ𝐼𝑡𝑠𝑖𝑛𝑠𝑗𝑑Γ.(3.6) Here, 𝐼Ψ𝐼=1 and Γ𝑠𝐼 is a part of the boundary of the substructure associated with node 𝐼. Γ𝑠Ψ𝐼𝑡𝑠𝑖 is equal to the nodal force 𝑓𝑠𝐼𝑖 at node 𝐼. If the normal 𝑛𝑠𝑗 in (3.6) is further approximated by the normal measured at nodes, the macroscopic stress tensor is given by𝜎𝑖𝑗=𝑟𝑉𝐼𝑓𝑠𝐼𝑖𝑛𝑠𝐼𝑗.(3.7) The normal, the nodal force, and the nodal displacement of the substructure are shown in Figure 3.

In the case of a linear elasticity, the nodal force 𝑓𝑠𝐼𝑖 is given in terms of the stiffness matrix and the nodal displacements:𝑓𝐼𝑖=𝑃𝐾𝐼𝑖𝑃𝑝𝑢𝑃𝑝=𝑃𝐾𝐼𝑖𝑃𝑝𝑟𝜀𝑝𝑞𝑛𝑃𝑞,(3.8) in which 𝐾𝐼𝑖𝑃𝑝 is the stiffness matrix of the boundary of the substructure, 𝑢𝑃𝑝 is the nodal displacement at node 𝑃, and 𝑛𝑃𝑞 is the normal at node 𝑃. It is trivial to mention that stiffness matrix 𝐾𝐼𝑖𝑃𝑝 in (3.8) can easily be obtained by using the standard static condensation. However, there is still a computational cost regarding the inversion of the reduced stiffness matrix during the static condensation process.

From (3.7) and (3.8), the macroscopic stress tensor is simply given by𝜎𝑖𝑗=𝑟2𝑉𝐼𝑃𝐾𝑠𝐼𝑖𝑃𝑝𝑛𝑠𝐼𝑗𝑛𝑠𝑃𝑞𝜀𝑝𝑞.(3.9) By comparing (3.9) to the general elastic constitutive relation, the macroscopic modulus tensor is given by𝐶𝑖𝑗𝑝𝑞=𝑟2𝑉𝐼𝑃𝐾𝑠𝐼𝑖𝑃𝑝𝑛𝑠𝐼𝑗𝑛𝑠𝑃𝑞.(3.10) Equation (3.9) can be extended to the case of nonlinear elastic materials. To do so, 𝐾𝑠𝐼𝑖𝑃𝑝 needs to be interpreted as the tangential (or algorithmic) stiffness matrix:Δ𝜎𝑖𝑗=𝑟2𝑉𝐼𝑃𝐾𝑠𝐼𝑖𝑃𝑝𝑛𝑠𝐼𝑗𝑛𝑠𝑃𝑞Δ𝜀𝑝𝑞.(3.11)

3.2. The Relation of this Method with the Microplane Model

It is instructive to discuss the similarity of this method to the microplane constitutive models developed for concretes and rocks [38, 39], and so forth. Figure 4 shows a sketch for the concept of the microplane model. The interaction between aggregates in a representative volume of concrete (Figure 4(a)) is modeled separately on a plane as shown in Figure 4(b). Strain tensor 𝜀𝑖𝑗 is decomposed to the normal and tangential components on a plane defined by normal 𝑛𝑖:𝜀𝑁=𝜀𝑖𝑗𝑛𝑗𝑛𝑖,𝜀𝑇=𝜀𝑖𝑗𝑛𝑗𝑡𝑖,(3.12) where 𝜀𝑁 and 𝜀𝑇 are the normal and tangential components and 𝑡𝑖 is the unit normal in a tangential direction on the plane. Such a projection is made on all the planes for the Gauss points on a unit sphere as shown in Figure 4(c). The stress vectors calculated on the individual planes are gathered to obtain the corresponding macroscopic stress tensor according to the principle of virtual work. Here a numerical integration on the sphere is necessary in the process.

Although it is beneficial to work with strain and stress vectors to develop the constitutive relation, taking into account the characteristics of the material, the constitutive relation becomes semiphenomenological because it is not easy to consider interactions between different microplanes in many cases. Also, the accuracy of the constitutive behavior depends on the order of the Gauss quadrature. So, the behavior of a microplane model is sometimes hinged on the numerical quadrature taken for the model although it is converging as the order of the quadrature increases [43].

The model proposed in this paper adopts the idea of the microplane model. However, we do not need to develop a constitutive model on each microplane. Instead, the behavior of the circular substructure to a given boundary displacement is directly calculated. The macroscopic stress tensor is constructed without a numerical integration by using the partition of unity of the finite element mesh as in (3.9) and (3.11).

4. Numerical Results

4.1. Plate Made of a Porous Material

Here we considered an elastic medium with many circular voids subjected to uniaxial tension as shown in Figure 5. It is assumed that the constitutive relation of each material point is completely described by Hooke’s law.

We considered two different length scales. The first length scale was the dimension of the specimen shown in Figure 5. The width and the height were 100 mm and 100 mm, respectively.

The stress at this scale was calculated from the second length scale. The second length scale was the dimension of the circular substructure radius which is 25 mm maximum. Note that the radius of the substructure is dependent on the porosity percentage and number of pores. The linear elastic material was used at this scale. The elastic modulus and Poisson’s ratio were 210 GPa and 0.3, respectively.

While the same elastic properties were used, various levels of the porosity from 5% to 20% were considered. The distribution of the linear elastic stress in the circular substructure is shown in Figure 6.

The overall properties may be obtained from the classical homogenization based on the rule of mixture. We solved this problem with the classical homogenization and the multiscale approach proposed in this paper. In computational homogenization, the macro-to-microtransition is achieved by enforcing the following condition:1𝐹=𝑉Ω𝑠Ω𝑠𝐹𝑑𝑉.(4.1) Here 𝐹 can be the any homogenized macroscopic system variable such as stress 𝜎𝑖𝑗 and 𝐹 is the corresponding microscopic variable. In other words, (4.1) essentially imposes a volumetric averaging of the system variables.

In order to analyze the convergence of the problem the normalized error in displacement and energy are computed by:𝑒𝐿2=𝑢𝑒𝑥(𝑥)𝑢(𝑥)𝐿2𝑢𝑒𝑥(𝑥)𝐿2=Ω𝑢𝑒𝑥(𝑥)𝑢(𝑥)2𝑑ΩΩ(𝑢𝑒𝑥(𝑥))2,𝑑Ω𝑒𝑒𝑛=𝑢𝑒𝑥(𝑥)𝑢(𝑥)𝑒𝑛𝑢𝑒𝑥(𝑥)𝑒𝑛=Ω𝐸𝜖𝑒𝑥(𝑥)𝜖(𝑥)2𝑑ΩΩ𝐸(𝜖𝑒𝑥(𝑥))2.𝑑Ω(4.2) Figure 7 shows the normalized error of the response of the system with respect to the reference model which is the detailed microstructure model. For a given porosity, as the number of pores was increasing the results became almost independent of the number of pores. Other factors such as size and geometry of the pores should not influence the result of the numerical homogenization scheme proposed here. Therefore we also investigated the effects of the number and size of the pores as well as the geometry of them and obtained only marginal changes in the results.

4.2. Cantilever Beam

The next example is an elastic cantilever beam as shown in Figure 8. The length of the beam is 48 m and the height is 12 m. The problem is in the plane stress condition. The beam is loaded with a parabolic traction at the end as shown in the figure. It was assumed that the constitutive relation of the problem should be obtained from a multiscale approach. The representative volume element of this example is assumed as a circular cell with a hole in the center was considered as shown in Figure 8. The substructure is not shown in Figure 8. The radius of the cell was 1 mm. The cell included a 0.5 mm radius of a hole in its center. The finite element mesh of the circular cell is shown in Figure 8. The elastic modulus and Poisson’s ratio of the material at the microscopic scale were 210 GPa and 0.3, respectively.

Here we would like to check the convergence rate of the homogenization model. This is necessary to certify that our model does not demand extremely fine meshes to give accurate solution with reasonable computation cost. Therefore with keeping the mesh of the circular cell unchanged, we increase the refinement of the mesh at the macroscopic scale. The error in the 𝐿2 norm and energy norm decreased as the mesh refinement increased; see Figure 9. The proposed multiscale method worked more successfully than the classical homogenization method.

4.3. Cracked Plate Made of a Bimaterial

The last problem is a mode 𝐼 crack problem shown in Figure 10(a). Due to symmetry, only a quarter of the problem was modeled and the rectangle is 40 m by 40 m and the length of the crack is 10 mm. It was considered that the material consisted of two different materials (labeled as 1 and 2, resp. in Figure 10(b)) and the microscopic structure could be modeled as a cell with an inclusion. The elastic properties for the microscopic cell are 𝐸1=60 GPa, 𝜈1=0.2, 𝐸2=26 GPa, and 𝜈2=0.2. Figure 11 shows that the energy norm and the stress intensity factor were quickly converging to the reference value as the mesh refinement was increasing. Here, the mesh of the cell at the microscopic scale was kept unchanged. Figure 12 shows the deformed mesh of this problem.

5. Conclusion

A new method for multiscale analysis was presented. This new method uses a circular cell in two dimensions (or a sphere in three dimensions) as a substructure at lower scales. The macroscopic strain tensor is projected to the boundary displacement of the cell by a kinematical constraint. This method is a generalization of the microplane model developed for the constitutive models of concrete. The stress on the macroscopic scale is obtained from the principal of virtual work applied between the traction and the stress of the circular substructure. Unlike the microplane model, the numerical integration on the boundary of the substructure is not needed. This new method is applied to several problems to show its performance. Because of the simple nature of this method, this method can be applied to many multiscale problems.

Acknowledgments

This project is also supported by a Korea University Grant. The financial support provided by German research foundation (DFG), under Grants CMS-0310596, 0303902, 0408359, Rolls-Royce Contract 0518502, Automotive Composite Consortium Contract 606-03-063L, and AFRL/MNAC MNK-BAA-04-0001 Contract, is gratefully acknowledged, too.