Meccanica DOI 10.1007/s11012-014-0031-x
MULTI-SCALE AND MULTI-PHYSICS MODELLING FOR COMPLEX MATERIALS
Particulate random composites homogenized as micropolar materials Patrizia Trovalusci • Maria Laura De Bellis Martin Ostoja-Starzewski • Agnese Murrali
•
Received: 26 May 2014 / Accepted: 17 July 2014 Ó Springer Science+Business Media Dordrecht 2014
Abstract Many composite materials, widely used in different engineering fields, are characterized by random distributions of the constituents. Examples range from polycrystals to concrete and masonry-like materials. In this work we propose a statistically-based scale-dependent multiscale procedure aimed at the simulation of the mechanical behavior of a two-phase particle random medium and at the estimation of the elastic moduli of the energy-equivalent homogeneous micropolar continuum. The key idea of the procedure is to approach the so-called Representative Volume Element (RVE) using finite-size scaling of Statistical Volume Elements (SVEs). To this end properly defined Dirichlet, Neumann, and periodic-type non-classical boundary value problems are numerically solved on the SVEs defining hierarchies of constitutive bounds. The results of the performed numerical simulations point out the importance of accounting for spatial randomness as well as the additional degrees of freedom of the continuum with rigid local structure.
P. Trovalusci (&) M. L. De Bellis A. Murrali Department of Structural and Geotechnical Engineering, Sapienza University of Rome, via A. Gramsci 53, 00197 Rome, Italy e-mail:
[email protected] M. Ostoja-Starzewski Department of Mechanical Science and Engineering, Institute for Condensed Matter Theory and Beckman Institute, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA
Keywords Random composites Cosserat continua Scale-dependent statistical homogenization Representative volume element
1 Introduction Various composite materials, encountered in nature and those extensively adopted in many engineering fields, are characterized by particulate random microstructures. Examples are polymer, ceramic, metal matrix composites or also concrete, granular materials and porous rocks, as well as masonry made of crushed stones with different dimensions, casually arranged in the mortar (Fig. 1). A key issue in mechanics of materials characterized by microstructural randomness is that the classical concept of the Representative Volume Element (RVE), well established in periodicity based homogenization techniques for many years [12, 21], loses its validity [15]. In the last few years, various procedures based on the solution of specific Boundary Value Problems (BVPs) have been proposed to perform classical homogenization for nonperiodic assemblies [4, 6, 19, 20, 22]. In order to account for the effects of the microstructural size, heterogeneous non-periodic materials have been also studied by extending the homogenization schemes to gradient-enhanced continua, although applied to a single fixed mesoscale [8, 9]. Stochastic approaches grounded on finite-size scaling homogenization have
123
Meccanica Fig. 1 Examples of particulate composites: ceramic/metal matrix composites (a, b); Roman concrete (c), tuffaceous rock (d)
proved to be among the most effective for identifying the RVE size and the overall constitutive moduli for linear elastic and thermoelastic, as well as non-linear and non-elastic, material systems [7, 14, 16]. In this paper we adopt the statistically-based scaledependent homogenization procedure developed in [26], in order to reproduce the actual microstructure of two-phase random media and to estimate the constitutive moduli of energy-equivalent continuum with rigid local structure (micropolar) [1, 13]. This procedure applies to random composites perceived micropolar continua both at the micro and macro level. The use of a Cosserat continuum at the local level is appropriate when the size of a heterogeneity is comparable to the characteristic length of its inner microstructure (lc ). In this way it is possible to account for local bending deformation mechanisms which can be prominent when particles sizes and/or inter-particle spacings are comparable to lc [2, 3, 17]. It is worth noting that when lc approaches zero (in the case of material isotropy, or at least orthotropy [23]) the microlevel is represented as a classical (Cauchy) continuum. A finite-size scaling of so-called statistical volume elements (SVEs) is used which approaches the RVE through two hierarchies of constitutive bounds, respectively stemming from the numerical solution of Dirichlet
123
and Neumann non-classical BVPs, set up on mesoscale material cells. For defining the classical and micropolar BVPs we use a generalized macro-homogeneity (Hill– Mandel type [5]) condition, which accounts for nonsymmetric stress and strain as well as couple-stress and curvature tensors. In particular, for a two-dimensional elastic medium made of a base matrix and a random distribution of disk-shaped inclusions of given density, two hierarchies of constitutive bounds are obtained by considering mesoscale test-windows of different sizes supposed to be placed anywhere in a random material domain. Under the hypotheses of statistical homogeneity and mean-ergodicity of the medium: the convergence trend of the bounds is detected as function of the SVE size; the RVE size is attained on the basis of a statistical criterion; the average homogenized, classical and micropolar, elastic moduli are estimated. The results of the numerical simulations, performed adopting Dirichlet, Neumann, and generalized periodic boundary conditions, point out the importance of taking into account the spatial randomness of the medium, and in particular the presence of inclusions that intersect the edges of the test windows. The value of accounting for the additional stress and strain measures of the micropolar continuum is also discussed.
Meccanica
2 Micropolar continuum
YK YC rij ¼ AYY ijhk ehk þ Aijhk ahk þ Aijhk khk ;
A material point of a micropolar continuum is characterized by its position and an orientation; the kinematical descriptors are displacements and rotations, represented by the vectors ðui Þ and (ui ), respectively. Within the framework of a linearized theory, in which displacements and rotations stand for velocities and angular velocities, respectively, and works for rate of works, the kinematics of the continuum is governed by the following relations: cij ¼ ui;j þ ekij uk ; jij ¼ ui;j :
ð1Þ
where ðcij Þ and ðjij Þ are the generally asymmetric strain and curvature tensors, respectively, and where eijk is the Levi–Civita tensor, with i; j; k ¼ 1; 3. The balance equations in the absence of body forces and couples are: sij;j ¼ 0; lkj;j þ ekji sij ¼ 0 ;
mi ¼ lij nj :
ð3Þ
In order to separately investigate the classical and micropolar components we divide the strain and stress tensors in their symmetric and skew-symmetric part. That is: cij ¼ eij þ aij ; sij ¼ rij þ bij ;
lij ¼
ACY ijhk ehk
þ
ACK ijhk ahk
þ
ACC ijhk khk
ð5Þ
;
with i; j; h; k ¼ 1; 3 and where (AYY ijhk ) is the classical KY KK YC constitutive tensor, while (AYK ), (A ijhk ijhk ), (Aijhk ), (Aijhk ), CY CK CC (AKC ijhk ), (Aijhk ), (Aijhk ) and (Aijhk ) are the micropolar constitutive tensors (i; j; h; k = 1; 3). For hyperelastic YY KK KK CC materials: ðAYY ijhk Þ = ðAhkij Þ, ðAijhk Þ = ðAhkij Þ, ðAijhk Þ = YK KY YC CY KC ðACC hkij Þ, ðAijhk Þ = ðAhkij Þ, ðAijhk Þ = ðAhkij Þ and ðAijhk Þ = CK ðAhkij Þ. In the specific case in which an isotropic behavior is taken into account, the Eq. (5) are specialized to: rij ¼ kehh dij þ 2leij ; bij ¼ 2lc aij ; 1 1 khh dij þ l2t 2l2c kji þ 2l2c kij ; lij ¼ l l2t W ð6Þ
ð2Þ
where ðsij Þ and ðlkj Þ are respectively the generally asymmetric stress and couple stress tensors. Denoting with ðti Þ and ðmi Þ the tractions and surface couples on the boundary of a control volume of outward normal ðni Þ, always with i; j ¼ 1; 3, we also have: ti ¼ sij nj ;
KK KC bij ¼ AKY ijhk ehk þ Aijhk ahk þ Aijhk khk ;
ð4Þ
where ðeij Þ and ðrij Þ are the classical symmetric strain and stress tensors, while ðaij Þ and ðbij Þ are the skewsymmetric strain and stress tensors characterizing, together with the curvature and the couple stress tensor ðjij Þ and ðlij Þ, a micropolar medium. The quantity aij ¼ 12 ðui;j uj;i Þ ekij uk is the relative rotation between the macrorotation and microrotation. The constitutive relations in the linear elastic anisotropic case are:
ði; j; h ¼ 1; 3Þ, where k and l are the Lame´ constants, lc is the micropolar shear modulus, lc is the characteristic length for bending, lt is the characteristic length for torsion and W is the dimensionless polar ratio, in agreement with [11].
3 Computational homogenization for random composites We study the scale-dependent effective response of heterogeneous random materials described as twodimensional and two-phase composites, under the assumption that the medium is characterized by statistical homogeneity and mean-ergodicity. We consider a simplified model made from a base matrix with randomly distributed disk-shaped inclusions of fixed radius d. Two antipodal cases, represented in Fig. 2, are considered: (a) a material with stiff inclusions in a soft matrix and (b) a material with soft inclusions in a stiff matrix. Thus, case (a) is representative of a so-called high contrast material, significant examples of which are metal ceramic composites (MCC) or concrete masonry-like materials; case (b) is representative of a low contrast material, as, for example, porous ceramic matrix composites (CMC) or porous filled rubble masonry or
123
Meccanica
Fig. 2 Scheme of simplified random models analyzed: a high contrast medium; b low contrast medium
magmatic rock structures. In the heterogeneous medium, isotropic linear elastic constitutive relationships are assumed for both material phases: rij ¼ kehh dij þ 2leij ; bij ¼ 2lc aij ; l3j ¼
2l l2c k3j
ð7Þ ;
(i; j = 1, 2) with evident meaning of the symbols. In order to determine the constitutive tensors components of the homogenized continuum and to detect the RVE size, we follow the statistical procedure presented in [26] and briefly recalled in the following steps: 1.
2.
3.
Set a nominal area fraction q, defined as the ratio between the total area of the inclusions and the area of a test window. Select a window size L and define the dimensionless scale factor d ¼ Ld. For each d determine the number of disks and simulate (uniform) random dispositions of disks’ centers using a hard-core Poisson point field (preventing the overlapping between disks). This corresponds to simulate (independent) realizations Bd ðxÞ of portions of the random medium sampled in a Monte-Carlo sense, x being an elementary event. For each Bd , solve both Dirichlet and the Neumann BVPs, consistently with a generalized macrohomogeneity condition, accounting for the presence of infinitesimal deformation gradients and curvatures: Z 1 ðrij eij þ bij aij þ lij jij ÞdV V Bd ðxÞ ð8Þ ¼ rij eij þ bij aij þ lij jij ;
123
4.
5.
where V is the volume of Bd and the overlined symbols define homogenized macroscopic measures, with the meaning of volume average quantities. The condition (8), in which the contributions of the classical and micropolar variables are considered separately, establishes a correspondence between the average internal work over a mesoscale window and the mechanical internal work density of the macromodel, expressed in terms of homogenized stress and strain measures. Repeat steps 2–3 until the confidence interval of the average homogenized constitutive parameters set at 95 %, evaluated over a normal standard distribution, is less than a small desired value. If the number of realizations necessary for ensuring the requirement in step 4 is small enough, this means that the values of the homogenized constitutive coefficients are distributed around their average with a vanishing variation coefficient. At this point the RVE is achieved. Otherwise choose an increased value of d and go to step 2.
The solution of the BVPs is required for any window randomly moving within the domain of the wall, which corresponds to a given realization of the heterogeneous medium Bd ðxÞ. In the following, the explicit dependence on x will be dropped for simplicity of notation. Let consider a square-shaped Euclidean region occupied by Bd (window), and let the origin of the reference frame coincide with the axes of symmetry of this region. The Dirichlet boundary conditions (D-BCs), consistent with the condition (8), can be written as: ui ¼ eij xj ;
1 u3 ¼ eij3 aij þ j3j xj 2
on oBd : ð9Þ
(i; j= 1, 2). The Neumann boundary conditions (NBCs) are: ti ¼ ðrij þ bij Þnj ; where mo3 ¼
m3 ¼ mo3 þ l3j nj
on oBd ; ð10Þ
Z
e3ij xi bjl nl is the moment imposed to oB
ensure the moment balance in the presence of skewsymmetric shear (i; j; l = 1,2). In addition to the BCs presented in Eqs. (9) and (10) a further case is herein considered: the adoption of
Meccanica
generalized periodic boundary conditions (P-BCs). Different authors [22, 25] have highlighted that, adopting P–BCs, it is possible to have a better estimate of the effective properties of the homogenized medium (for a single mesostructure size), even if the microstructure is not characterized by a periodic distribution of the inclusions, at least limited to the case of classical continua. In this work an extension of the standard P-BCs conceived for classical continua is proposed to handle with the case of micropolar continua. Generalized conditions are therefore deduced and briefly presented below. Let xþ i and xi (i ¼ 1 ; 2) be the coordinates of two corresponding points (pair) located at opposite edges of Bd . The P-BCs can be written as follows. First, considering to assign the classical components of the homogenized macroscopic strain eij , we have (i, j = 1, 2): –
conditions imposed on the vertices of oBd ui ¼ eij xj ;
–
ð11Þ
conditions imposed on the pairs of points on the edges of oBd þ uþ i ui ¼ eij ðxj xj Þ ;
ð12Þ
Considering the skew-symmetric shear strain aij , we have (i, j = 1, 2): –
conditions imposed on the vertices of oBd 1 u3 ¼ eij3 aij 2
–
ð13Þ
conditions imposed on the pairs of points on the edges of oBd : uþ 3
u 3
¼ 0:
ð14Þ
Finally, associated to the imposition of the curvature j3i there are (i = 1, 2): –
–
ð15Þ
conditions imposed on the pairs of points on the edges of oBd parallel to the ith axis uþ 3 u3 ¼ 0 ;
conditions imposed on the pairs of points on the edges of oBd orthogonal to the ith axis uþ 3 þ u3 ¼ 0 ;
ð17Þ
The same statistical procedure detailed before has then been adopted by using these PBCs and solving the corresponding boundary value problems for any mesoscale window until attaining the convergence to the RVE detected as described in step 5.
4 Numerical results The BVPs described in Sect. 3 have been numerically solved by using COMSOL Multiphysicsany window randomly moving within the domainr , a finite element code in which the related field equations can be directly implemented. Unstructured meshes of quadratic Lagrangian triangular finite elements have been adopted. This study is focused on the contrast between material properties of the two samples considered, which may be representative of different kind of materials. Table 1 reports the adopted material parameters in a dimensionless form, expressing the ratio between corresponding quantities of inclusions and matrix. The homogenized model is generally anisotropic and the constitutive coefficients are those reported in Eq. (5), specialized for the two-dimensional medium (i; j; h; k ¼ 1; 2). In the following we focus our attention on the most significant components for the particular heterogeneous medium selected, correYY
KK
CC
sponding to (Aijhk ), (Aijhk ) and (Aijhk ) and consider these elastic coefficients: A KK
KK A1212
YY
CC
= ðA1111 þ A2222 Þ=2, CC
classical; A = and A = trA , micropolar. The mesoscale window Bd ideally corresponds to a portion of the actual random medium in which Table 1 Ratios between material parameters of inclusions and matrix
conditions imposed on the vertices of oBd u3 ¼ j3i xi
–
ð16Þ
Material
Parameters ki =km
li =lm
lci =lcm
lci =lcm
a
46
4.93
4.93
101
b
0.021
0.202
0.202
101
123
Meccanica
inclusions are not prevented from intersecting the window edges. Thus, the numerical simulations are performed by taking into account non-homogeneous boundaries (crossing inclusions). We also consider the less realistic case of homogeneous boundaries (noncrossing inclusions). The comparison between the homogenized responses obtained by performing numerical simulations for the two cases, either applying Dirichlet and Neumann boundary conditions, allows us to emphasize the influence of positions of the inclusions with respect to the windows’ boundary. The solutions under D-BCs and N-BCs are also compared to the solutions obtained by applying periodic boundary conditions. These P-BCs are herein adopted limited to the case of homogeneous boundaries. In the case of non-homogeneous boundaries, in fact, the generalized P-BCs cannot be applied in the large, since the crossing inclusions are randomly distributed along the window edges and an priori oneto-one correspondence of displacements and rotations between pair of points cannot be established. Figure 3 reports the average of the classical YY
coefficient A versus the scale parameter d, for both materials (a) and (b). This value is normalized to the average value of the convergence coefficient evaluYY
ated in the case of crossing inclusions, hARVE i, obtained imposing D-BCs and N-BCs in the case of non-crossing inclusions. The convergence trend to RVE depends on whether inclusions cross or do not cross the windows’ boundaries. In particular, for the material (a) the RVE size, dRVE , is equal to 20 in the case of crossing inclusions, while it is dRVE ¼ 25 in the case of non-crossing inclusions. The BVPs results obtained under P-BCs show a faster convergence, with respect to the lower/upper bounds approach, to a value of the elastic modulus between the Dirichlet and Neumann BVPs solutions (closer to the Dirichlet solution) obtained in the case of non crossing inclusions. The material (b) shows a slower convergence trend. Accordingly, the RVE is attained for dRVE ¼ 25 in the case of crossing inclusions, while in the case of non-crossing inclusions dRVE [ 25. Analogously to the case of material (a), the results under P-BCs converge, as d increases, to a value near the convergence value obtained adopting D-BCs and N-BCs in the presence of non-crossing inclusions. Figure 4 reports the micropolar results in terms of the average of the coefficient A
123
KK
(normalized to the
KK
average convergence value hARVE i obtained for crossing inclusions) versus the scale parameter d, always for both materials (a) and (b). The two material coefficients exhibit convergence trends similar to the YY
KK
coefficients A . Note that the terms A do not vanish as d increases, indicating that the micropolar modulus KK
ARVE is significant also in the presence of inclusions of small size. Also in this case, the P-BCs provide results that are between those obtained with D-BCs and N-BCs for homogeneous boundaries. Thus, a faster convergence is guaranteed. Figure 5 shows the micropolar results in terms of the average of the homogenized characteristic length qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi CC KK parameter lc ¼ A =A (normalized to the average convergence value hlcRVE i) evaluated for inclusions crossing the windows boundary) versus the scale parameter d. Both materials exhibit differences between the curves obtained in the case of crossing and non-crossing inclusions greater than in the classical case. For the material (a) the RVE is attained at dRVE ¼ 15 in the case of crossing inclusions and at dRVE ¼ 20 in the case of non-crossing inclusions. For the material (b) the RVE is attained at dRVE ¼ 20 in the case of crossing inclusions and at dRVE ¼ 25 in the case of non-crossing inclusions. The results under the generalized P-BCs seems to confirm what highlighted in the case of classical continuous, i.e. a faster convergence to a value of the micropolar elastic modulus between the Dirichlet and Neumann BVPs solutions (closer to the Dirichlet solution) obtained in the case of homogeneous boundaries. This value significantly differs from hlcRVE i obtained in the case of non-homogeneous boundaries. Moreover, it is observed that, for both materials (a) and (b), when d increases, hlc i tends to the value of the characteristic length lc of the matrix (material (a): hlcRVE i ¼ 0:1; material (b): hlcRVE i ¼ 1Þ. This shows micropolar bending effects to be weaker in the medium (a) than in the medium (b). These findings are in agreement with some experimental results [10, 11]. Overall all these results show that with the proposed scale-dependent procedure, based on the estimation of hierarchies of upper and lower bounds, the RVE size and the corresponding effective, classical and micropolar, constitutive moduli can be statistically detected for materials with different contrast between elastic
Meccanica
5 Final remarks The scale-dependent homogenization of composite materials made of random distribution of inclusions in a matrix described as micropolar continua has been here reported. The proposed procedure exploits the solution of hierarchies of bounds stemming from the
1.3
/RV E
1.2 1.1 1 0.9 0.8 D−BC no−cr D−BC cr N−BC no−cr N−BC cr P−BC no−cr
0.7 0.6 0.5
5
10
15
20
25
δ = L/d Medium (b) with soft inclusions 1.3
YY
>/RV E
1.2
RV E
KK
>/RV E
1.6
1.1 1 0.9 0.8 D−BC no−cr D−BC cr N−BC no−cr N−BC cr P−BC no−cr
0.7 0.6 0.5
5
10
15
20
1.4
1.2
1
0.8
0.6
25
5
10
δ = L/d
< lc > / < lc >RV E
>RV E KK
>/
