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.
(a)
(b)
(c)
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 [7–9] 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 [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 [15–20]. 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 [23–28]. 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 can be applied; see [30–35] 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 [38–41], 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 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 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 such that with 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 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, where is the stress tensor at the macroscopic scale, is a suitably chosen weight function such that , 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: 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 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: Here, 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 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: 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 By comparing (3.9) to the general elastic constitutive relation, the macroscopic modulus tensor is given by 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:
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 : 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.
(a)
(b)
(c)
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.
(a) 5% porosity
(b) 20% porosity
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: 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: 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 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.
(a) Change of the norm |
(b) Change of the energy
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 GPa, , GPa, and . 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.
(a)
(b)
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.