- Email: [email protected]
RVE for media with randomly distributed inclusions
CHAPTER 9
RVE for media with randomly distributed inclusions 9.1 Introduction Modern materials often exhibit different internal architectures at different length scales or in different perspectives. Some of those architectures might be of completely irregular physical and geometric characteristics, such as that in particulate reinforced composites, in which reinforcing particulates are usually dispersed at random, or the random fibre distribution in the transverse cross-section of UD composites. Micromechanical modelling of such problems requires the use of representative volume elements (RVEs). The definition of an RVE and its associated considerations have been presented in Chapter 4, followed by a critical review of commonly observed misperceptions and misuse of RVEs in the literature in Chapter 5. Having reached the level of understanding thus far, the need is apparent for RVE formulation and analysis techniques that are logically derived based on the correct concepts and deliver results that are mathematically accurate. The task of this chapter is to offer such an approach to address the problem and elaborate on the steps taken in order to deliver the subject without ambiguity. Exact solution would be obtained if the problem could be analysed over an infinite domain, or the correct boundary conditions had been prescribed if a finite domain is employed as an RVE. However, none of these options is practical in general. The proposed approach is based on the principle of Saint-Venant. In the context of the present application, it states that the results obtained over the regions within the RVE sufficiently distant from the boundary, where inaccurate boundary conditions are prescribed, should be close enough to the exact solution, provided that the forces associated with such incorrect boundary conditions are statically equivalent to those appearing in the infinite domain or those in the exact boundary conditions. However, the principle of Saint-Venant itself does not provide any estimate of how large the distance from the boundary should be. In fact, it varies from problem to problem in general and the characteristic measure of this Representative Volume Elements and Unit Cells ISBN: 978-0-08-102638-0 https://doi.org/10.1016/B978-0-08-102638-0.00009-8
© 2020 Elsevier Ltd. All rights reserved.
319
j
320
Representative Volume Elements and Unit Cells
distance, or decay length, as defined in Chapter 4, will have to be determined in an ad hoc manner for a given type of problems.
9.2 Displacement boundary conditions and traction boundary conditions for an RVE The description of uniform displacement and uniform traction boundary conditions as referred to in Chapter 4 is in fact rather simplistic in a sense that it does not provide instructions of how such boundary conditions are to be implemented. The elaborated description of these boundary conditions will therefore be provided here to benefit potential users. Consider a rectangular area, 0 ya and 0 z b, as the volume element in the yz-plane in a 2D space for the ease of presentation, which can be easily extended to a cuboid in a 3D space. Assuming a uniaxial stress s0y is to be prescribed in the y-direction in order to determine the effective Young’s modulus and Poisson’s ratio in this direction, the traction boundary conditions can be given as follows. py y¼0 ¼ s0y and pz jy¼0 ¼ 0; and pz jy¼a ¼ 0; py y¼a ¼ s0y (9.1) and pz jz¼0 ¼ 0; py z¼0 ¼ 0 and pz jz¼b ¼ 0; py z¼b ¼ 0 where py and pz are the traction components. If one is interested in the determining the effective shear modulus, the traction boundary conditions for an RVE under a pure shear stress s0yz in the yz-plane are to be applied as and pz jy¼0 ¼ s0yz py y¼0 ¼ 0 py y¼a ¼ 0 and pz jy¼a ¼ s0yz (9.2) py z¼0 ¼ s0yz and pz jz¼0 ¼ 0 and pz jz¼b ¼ 0: py z¼b ¼ s0yz For FEM applications, these conditions will have to be topped up with appropriate constraints to rule out rigid body motions, in particular, rigid body rotations. The correct method of constraining rigid body rotations when the domain deviates from orthogonal regularity is elaborated in
321
RVE for media with randomly distributed inclusions
Section 5.5 of Chapter 5. Regarding the application of tractions, an important point is to be made about the vanishing traction components as involved in (9.1) and (9.2). Since traction boundary conditions are natural boundary conditions, they are prescribed when loads are applied. Loads of zero magnitude imply the absence of loads and therefore they do not need to be applied explicitly in an FE analysis. In other words, vanishing traction boundary conditions should be left unattended and this is strictly correct as far as FE modelling is concerned, whilst they will be satisfied approximately through energy minimisation in the same way as equilibrium conditions are satisfied, as well as the nontrivial boundary conditions. In (9.1) and (9.2), what needs to be applied is therefore reduced to the following. py y¼0 ¼ s0y and py y¼a ¼ s0y for uniaxial tension; and (9.3) pz jy¼0 ¼ s0yz py z¼0 ¼ s0yz
and and
pz jy¼a ¼ s0yz py z¼b ¼ s0yz
for pure shear.
(9.4)
Note that if one wishes to apply uniaxial stress s0z in the z-direction, boundary conditions for such case will be similar to those given by (9.3) but will have to be applied on surfaces z¼0 and z¼b. The so-called uniform displacement boundary conditions are not usually prescribed purely in terms of displacements. The boundary conditions for the same loading case as (9.1) and (9.2) are actually defined in a mixed form as follows. vjy¼0 ¼ 0
and
pz jy¼0 ¼ 0
and pz jy¼a ¼ 0 vjx¼a ¼ aε0y py z¼0 ¼ 0 and pz jz¼0 ¼ 0 py z¼b ¼ 0 and pz jz¼b ¼ 0 py y¼0 ¼ 0 and wjy¼0 ¼ 0 py y¼a ¼ 0
and
vjz¼0 ¼ 0
and pz jz¼0 ¼ 0
vjz¼b ¼ bg0yz
and
wjy¼a ¼ 0 pz jz¼b ¼ 0
for uniaxial tension; and
for pure shear.
(9.5)
(9.6)
322
Representative Volume Elements and Unit Cells
These boundary conditions should be topped up with the constraints to eliminate rigid body motions remaining in the system, viz. translation in the z-direction for (9.5) and rigid body rotation for (9.6). Ignoring the vanishing natural boundary conditions, one is left with the following, respectively. vjy¼0 ¼ 0 (
and
vjy¼a ¼ aε0y
for uniaxial tension; and
wjy¼0 ¼ wjy¼a ¼ 0 vjz¼0 ¼ 0
and
vjz¼b ¼ bg0yz
for pure shear
(9.7) (9.8)
One can apply uniaxial strain ε0z in a similar manner as ε0y is prescribed in (9.7). Note that the vanishing traction boundary conditions in (9.5) and (9.6) should not be replaced by vanishing displacements. Otherwise, they will correspond to a completely different loading conditions. It should also be pointed out that although (9.1) and (9.2) or (9.3) and (9.4) are unique forms of prescribing the relevant traction boundary conditions, the way of constraining rigid body motions as a necessary part of the imposition of boundary conditions can take different forms, and rigid body motion constraints have not been included in (9.1) and (9.2) or (9.3) and (9.4). On the other hand, (9.5) and (9.6) or (9.7) and (9.8) are not unique forms of prescribing displacement boundary conditions for the relevant effective stress states. This is due to the fact that any combination of rigid body motions can be superimposed on top of them to present a different appearance without changing the stress state, as has been explained in Section 6.2 of Chapter 6. For instance, one can have 1 vjy¼0 ¼ aε0y 2
and
1 vjy¼a ¼ aε0y 2
for uniaxial tension; and (9.9)
wjy¼0 ¼ 0
and
1 wjy¼a ¼ ag0yz 2
for pure shear (9.10) 1 0 vjz¼0 ¼ 0 and vjz¼b ¼ bgyz 2 They will be to exactly the same effects as (9.7) and (9.8). Any lack of uniqueness is a source of potential confusions and it is the responsibility of the users to apply the displacement boundary conditions correctly.
RVE for media with randomly distributed inclusions
323
In much the same way as in UCs, when (9.7) and (9.8) are employed, average strains ε0y and g0yz serve as Kdofs at which loads can be applied. The concept of Kdofs and their applications have been introduced in Section 6.6 of Chapter 6. Whilst displacement boundary conditions as provided above could be perfectly correct for UCs, as obtained from the symmetry conditions due to the regularity of the architecture at the lower length scale, the traction boundary conditions are usually incorrect unless the material is homogeneous even at its lower length scale. In absence of regularity as is usually the case for RVEs, even the displacement boundary conditions cannot be prescribed correctly, in general, let alone the traction boundary conditions. Any application of such boundary conditions will have to be considered as some kind of approximation. The extent of such an approximation determines the validity of the RVE introduced, as will be the subject of this chapter. Hill (1963) suggested that the volume element was representative if one could obtain sufficiently close results under both the uniform traction and the uniform displacement boundary conditions. As commented in Chapter 4, this is unnecessary, as will be elaborated fully in the next a few sections.
9.3 Decay length for boundary effects Without perfect symmetry, neither uniform tractions nor uniform displacements will give exact boundary conditions for the volume element concerned. However, either of them will represent a statically equivalent uniaxial loading condition if homogeneity is present at the upper length scale. According to the principle of Saint-Venant, either of them should result in perfectly correct stress and strain distributions an appropriate distance into the volume element from its boundary. This distance varies from point to point along the boundary. The maximum value of this distance for all relevant loading conditions is the decay length, that has already been introduced in Chapter 4. The inner part of the volume element a decay length away from its boundary in all directions forms a subdomain of the volume element. As long as the original volume element is large enough to allow such a subdomain to be also representative in terms of volume fractions of constituents, the subdomain will be an RVE, a perfect one. The effective properties of the material should not be extracted directly from the original RVE itself, but from the subdomain where the stress and strain distributions are free from the effects of incorrect boundary conditions, in
324
Representative Volume Elements and Unit Cells
terms of either traction or displacement. The question that remains is what should be the decay length. The principle of Saint-Venant does not provide any clear guideline for it, and in fact, the decay length varies from material to material and from one effective property to another, hence it will have to be determined by the user on a case-by-case basis. In order to estimate the decay length for UD composites, one with regularly packed fibres in its transverse cross-section is considered as an example. In a composite so chosen, symmetry is present and hence correct boundary conditions can be obtained if the border of volume element is selected along symmetry lines as a stack of repeating cells. In order to examine the effect of inaccurate boundary conditions, one side of the stack of such cells is deliberately truncated away from the symmetry line as shown in Fig. 9.1, such that no correct boundary conditions can be prescribed on it. Two analyses were conducted with each stack of cells, one involving uniform traction boundary conditions, as defined by (9.3), whilst another was carried out with uniform displacement (9.7) being prescribed. Both types of boundary conditions are statically equivalent to uniaxial tension. Comparing the stress distributions from the two analyses in the truncated cells, one can observe how far the disturbance introduced by the incorrect boundary conditions propagates as a way of evaluating the decay length. The results shown in Figs. 9.2 and Fig. 9.3 are von Mises stress contours plotted over a squarely packed array and a hexagonally packed array, respectively, where Fig. 9.2(A) and Fig. 9.3(A) were obtained with uniform displacement being applied, whilst in results shown in Fig. 9.2(B) and Fig. 9.3(B) were obtained with uniform traction being prescribed. The von Mises stress is not usually relevant for composites in general. However, in this case, both the fibre and matrix are assumed to be isotropic and the von Mises stress should be meaningful in the context of the present discussion. Only complete cells have been shown in Figs. 9.2 and 9.3, since the stress distributions in the truncated ones are deemed to be wrong. In Figs. 9.2 and 9.3, the columns with title ‘Case-0’ show stress contour over the stacks that did not involve the truncated cells hence had correct boundary conditions applied to them. Such stacks were analysed in order to compare the results with those obtained for truncated stacks. The remaining columns were obtained under incorrect boundary conditions, as each of those stacks had on the top an incomplete cell truncated at different locations as shown in Fig. 9.1. The pattern of contour plots in the column of Case-0 in Fig. 9.2(A) is identical for every cell in the stack, because Case-0 represents a correct solution to the problem. Comparing the patterns in Case-1, Case-2 and Case-
RVE for media with randomly distributed inclusions
325
Fig. 9.1 Arrays of periodic cell of (A) square packing (B) hexagonal packing.
3 with that in Case-0, it can be observed that the further away the cells are from the top border where incorrect boundary conditions are prescribed, i.e. moving from cell No. 1 to cell No. 5, the closer the comparison of the stress contours with that in Case-0. In fact, there is no noticeable difference
326
Representative Volume Elements and Unit Cells
(a) L1 1→
(b) Case-0 No.1
Case-1
Casse-2
Case-3
Case-0 C
L1 →
Case-1
Case-2
Case-3
No.1
L2 2→
No.2 2
L2 →
No.2
L3 →
No.3 3
L3 →
No.3
L4 L →
No.4 4
L4 →
No.4
L5 5→
No.5 5
L5 → No.5
Fig. 9.2 Von Mises stress contours over complete periodic cells from the square packed arrays (A) under prescribed uniform displacement and (B) under prescribed uniform traction.
between stress contours from cell No.2 onwards. Comparing respective cases in Fig. 9.2(A) and (b), it becomes obvious that relatively, the boundary effects decay over a longer distance if a uniform traction is prescribed, although disturbances in the latter case settle down eventually by cell No. 4 or so. The same trends are repeated in Fig. 9.3 for hexagonal packing where the plots have been arranged in exactly the same order as those in Fig. 9.2. In other words, under uniform displacement the decay length is approximately equal to the length of one cell, whilst under uniform traction it can be several cells. Apparently, the boundary conditions corresponding to prescribed uniform displacement represent a better approximation to the exact boundary conditions than those of prescribed uniform traction. Therefore, it is advisable that the displacement boundary condition should be employed. In this case, the typical decay length will be a couple of characteristic lengths. More justifications can be found in (Wongsto and Li, 2005). For regular patterns, the characteristic length is the cell size. In the case of
327
RVE for media with randomly distributed inclusions
(a)
(b) Case-0
L1 o
Case-11
Casse-2
Case-3
L1 o
Case-0
Case-1
C ase-2
Case-3
No o.1
No.1
L2 o
L2 o
No o.2
No.2
L3 o
L3 o
No o.3
No.3
L4 o
L4 o
No o.4
No.4
L5 o
L5 o
No o.5
No.5
Fig. 9.3 Von Mises stress contours over complete periodic cells from the hexagonally packed arrays (A) under prescribed uniform displacement and (B) under prescribed uniform traction. Generating the microstructures of randomly distributed physical and geometric features.
irregular patterns, it can be expected to be the average distance between the centres of two adjacent fibres. This will be further verified in Section 9.5.
9.4 Generation of random patterns In order to facilitate theoretical construction of RVEs, a means of generating a randomly distributed features is desirable. A UD composite having random fibre distribution over its transverse cross-section will be taken as a special case in a 2D space. The process of reproducing the geometry of cross-section is equivalent to generating randomly spaced, nonoverlapping circular discs in a plane within a domain. There could be different approaches to achieve this goal. Digitising microscopic photographs such as that shown in Fig. 4.1 is an obvious one, but the procedure could be computationally costly and time-consuming. Employing a suitable numerical algorithm would be more effective in this respect. One possibility is to adopt a ‘coin-dropping’, or random insertion process, as was
328
Representative Volume Elements and Unit Cells
implemented by Bulsara et al. (1999). Whilst the random nature of the obtained distribution is obvious, it may be difficult to achieve a specified disc (constituent) volume fraction, especially at high volume fractions. Constituent volume fraction is often a dominant parameter in determining the performance of the composite and a reasonable control of it is essential. A different process was employed by Gusev et al. (2000) using a Monte Carlo method. An alternative scheme was proposed in Wongsto and Li (2005), which is capable of reproducing any practical volume fraction and yet maintaining the random nature in distribution. The outline of the procedure is as follows. A set of regularly packed discs of radius R in the yz-plane is generated. An example of such set with hexagonally packed discs is shown in Fig. 9.4. The spacing between discs, 2b, can be calculated making use of the constituent volume fraction, Vf, as 2 pR ffiffi . Next, two frames are introduced and fixed in the plane, Frame b2 ¼ 2p 3V f
1 and Frame 2, with the latter being placed inside the former. The distance from the border of Frame 2 to that of Frame 1 should be equal to a decay length, which, according to the conclusion reached in the previous section, is a few characteristic lengths, b. As in Fig. 9.4, it is approximately 4b. These frames do not have to be perfect squares, but they both have to contain
Frame-1
Frame-2
Fig. 9.4 Hexagonally packed discs at Vf ¼ 65% involving 105 discs.
RVE for media with randomly distributed inclusions
329
sufficient number of discs. They have been chosen to be square, i.e. b¼a, as shown in Fig. 9.4 for the sake of argument. As a next step, each disc is moved in a random manner as illustrated in Fig. 9.5. A random angle q between 0 and 360 is generated first to determine the direction for the disc to shift. Along this direction, a maximum distance of shift can be found, which is the smaller of the distance to the border of Frame 1 and that to the point when the disc comes into contact with another, marked as r in Fig. 9.5. Then a second random number, k, is generated between 0 and 1 and the actual distance of shift is given as kr. This process is performed for every disc to complete one iteration. One iteration is effectively one shake of the layout. A sufficiently large number of iterations will have to be made, e.g. 250, before the original regularity disappears. After that, a check on the volume fractions within both frames will be conducted after each iteration. The iterations terminate when the volume fractions in both frames are close enough to the original one, allowing a small discrepancy, e.g. 2% error, as the criterion for terminating iterations. The procedure can be easily programmed and, given parameters of R, Vf and the number of discs to be involved, a set of coordinates can be obtained as the centres of the randomly distributed discs. As examples of the applications of the scheme outlined above, three cases have been generated at a volume fraction of 65% as shown in Fig. 9.6(AeC). They are of the similar volume fraction to that shown in Figure 4.1 and share all the features present in Figure 4.1. A further case was generated at Vf ¼ 60% as shown in Fig. 9.6(D), for which some experimental data are available (Soden et al., 1998). Computationally, this procedure is very efficient, and any of the cases can be generated literally in a matter of seconds on an ordinary PC. Cases as shown in Fig. 9.6 were meshed and analysed with the results being presented in the next section.
U kU
T
z y
Fig. 9.5 Schematic diagram for perturbation process.
330
Representative Volume Elements and Unit Cells
(a) Vf1 = 65.69%, and Vf2 =64.17%
(b) Vf1 = 65.84% and Vf2 = 64.26%
(c) Vf1 = 66.38% and Vf2 = 64.59%
(d) Vf1 = 61.00% and Vf2 = 59.40%
Fig. 9.6 Random fibre distribution obtained for a case with original volume fraction Vf ¼ 65% after (A) 265 iterations; (B) 286 iterations and (C) 382 iterations. Fibre distribution for a case with Vf ¼ 60% after 262 iterations e plot (D). Volume fractions within Frame 1 and Frame 2 are denoted as V1f and V2f , respectively.
9.5 Strain and stress fields in the RVE and the subdomain The cases as generated in the previous section represent transverse cross-sections of UD composites when randomly distributed discs are viewed as fibres and the surrounding area as the matrix. To proceed with the subsequent micromechanical FE analysis, area delimited by Frame 1 will be treated as the RVE. For the purpose of extracting and processing
RVE for media with randomly distributed inclusions
331
the results afterwards, the areas inside and outside Frame 2 are meshed individually, whilst the compatibility of the mesh will be maintained across the border of Frame 2 so that the part of the mesh inside Frame 2 can be extracted as a subdomain. Examples of these meshes are shown in Figs. 9.7 and 9.8, where respective RVEs correspond to random distributions presented in Fig. 9.6(A) and (D). In Fig. 9.7, a further frame was introduced within Frame 2. This does not compromise the validity of the mesh except that it had an additional consideration during meshing. The purpose of it will be clear when the results are discussed. Any gap between two neighbouring fibres will be meshed into at least two elements across the gap for representativeness of the mesh. In the unlikely event of two touching fibres, compromises will have to be made using triangular elements of very sharp corners. The analysis is carried out based on the generalised plane strain idealisation as has been described in Chapter 6 in the context of UC applications. The quadratic triangular and quadrilateral generalised plane strain elements CPEG6 and CPEG8 of Abaqus (2016) were used for the longitudinal/transverse tension and transverse shear deformation whilst heat conduction elements DC2D6 and DC2D8 were adapted for the anticlastic problem of longitudinal shear based on an analogy as described in Section 6.4.1.1 of Chapter 6.
Fig. 9.7 The mesh for the RVE with fibres distributed at random for the case of Vf¼65%.
332
Representative Volume Elements and Unit Cells
Fig. 9.8 The meshes for one of the RVEs with fibres distributed at random at Vf¼60%.
It has been established in Section 9.3 that imposition of a uniform normal displacement is a more advantageous method of prescribing the loading to the RVE. Therefore, to achieve a macroscopically uniaxial stress state in the y-direction, displacements as defined in (9.7) were applied along the relevant sides of the boundary. This can be conveniently achieved by making use of the Kdof ε0y . As has been established in Chapter 7, a uniform displacement boundary condition can be prescribed by applying a concentrated ‘force’ As0y , where A¼a2 is the planar area of the RVE, to the Kdof instead of a prescribed nodal ‘displacement’ without changing the nature of boundary conditions, provided that displacements on the boundary have been associated with the Kdofs. This ensures a uniform displacement imposed on the boundary whilst leaving its magnitude for the subsequent analysis to determine as the ‘nodal displacement’ at the Kdof concerned. To facilitate the analysis, for all the cases with 65% fibre volume fraction, both the fibre and matrix are assumed to be isotropic and homogeneous with elastic properties Em ¼ 1 GPa, nm ¼ 0.3 and Ef ¼ 10 GPa and nf ¼ 0.2 for matrix and fibre, respectively. For the case with 60% fibre volume fraction, material properties of both fibre and matrix are listed in Table 9.1 (Soden et al., 1998).
333
RVE for media with randomly distributed inclusions
Table 9.1 Constituents properties of the Silenka Eglass 1200 tex and epoxy UD composite.
Matrix Young’s modulus Matrix Poisson’s ratio Fibre Young’s modulus Fibre Poisson’s ratio Fibre volume fraction
3.35 (GPa) 0.35 74 (GPa) 0.2 60%
The analyses were carried out on the domain Frame 1. To obtain the pattern of deformation and stress distribution in the RVE, a load as the ‘concentrated force’ equivalent to 1 MPa average uniaxial stress in the y-direction was applied to it. The results extracted from the subdomain Frame 2 were considered to represent the correct responses of the material, free from the effects of incorrectly prescribed boundary conditions along the border of Frame 1. The analysis and data extraction process have been applied to all the four cases shown in Fig. 9.6. The results are presented and discussed below. The von Mises stress contour plot over the subdomain Frame 2 corresponding to the case with Vf¼65% as shown in Fig. 9.6(A) is displayed in Fig. 9.9. The contour plot is shown in a deformed configuration with the magnitude of deformation being amplified by a factor of 300. The areas with the highest values of stress are always found at interfaces, especially
Fig. 9.9 Contours of von Mises stress in Frame 2 as a part of RVE for Case-1 under macroscopically uniaxial transverse tension.
334
Representative Volume Elements and Unit Cells
in areas where fibres are close to each other in the direction of loading (Zou and Li, 2002). The distortion of the shape of fibre cross-section is small relative to overall deformation because of fibre stiffness being higher than that of the matrix. Therefore, the deformation is mostly taken by the matrix although stresses in the fibres tend to show higher levels than those in the matrix. This explains why observed stress concentrations are so located. As a result of deformation, the edges of Frame 2 do not remain straight. This provides an indication how wrong it would be if Frame 2 was analysed alone and straight edges had been assumed when prescribing displacement boundary conditions. To further emphasise this, a case was generated where the domain Frame 2 was analysed directly, with uniform displacements prescribed along the straight edges of this domain. The stress contour for this case is shown in Fig. 9.10 with the same deformation magnification factor of 300. Comparing this case with results obtained from the more accurate analysis conducted on Frame 1 from which Frame 2 was extracted, as shown in Fig. 9.9, significant difference can be found. The maximum von Mises stress in the incorrect analysis were found to be 9.879 MPa in the fibres and 4.501 MPa in the matrix whilst their counterparts obtained in the more accurate analysis were 14.371 MPa and 10.997 MPa, respectively. A potential but significant risk of employing incorrect analysis is that if one uses the von Mises criterion with a given yield stress for the matrix, the load for the onset of plastic deformation in the matrix will be miscalculated.
Fig. 9.10 Von Mises stress contour over the Frame 2 alone analysed with inaccurate boundary conditions.
RVE for media with randomly distributed inclusions
335
In this particular case, for instance, the respective load was underestimated by approximately a factor of 2.5. Given the erroneous analysis as shown in Fig. 9.10, reasonable accuracy can still be achieved in the subdomain within Frame 2 in Fig. 9.7 as mentioned before. To demonstrate this, the von Mises stress contour plots for this same subdomain but from different analyses, one on Frame 1 and the other on Frame 2, are shown in Fig. 9.11. Both plots share many common features, although the distance to the borders of Frame 2 is not quite equal to a decay length and the innermost frame may not be large enough to be representative. Some of the differences in the patterns of the contour
Fig. 9.11 Frame 3 extracted from (a) Fig. 9.10 and (b) Fig. 9.11, both with a deformation magnification factor of 300.
336
Representative Volume Elements and Unit Cells
plots are due to the difference in the scales of the stress levels as shown in the legends. The consistency in their trends reinforces the justification made to the proposed methodology based on the fact that under approximate boundary conditions, correct solutions can be obtained from an inner subdomain at least one decay length away from the boundary where the approximate boundary conditions are prescribed. The underlying mechanical principle is that of Saint-Venant.
9.6 Post-processing for average stresses, strains and effective properties As was mentioned in Chapter 5, a proper description of the procedure for post-processing the results from the analyses of RVEs and UCs has been a topic often avoided the literature. In the case of RVEs, if stresses and strains are to be extracted from a subdomain, one cannot make use of Kdofs when post-processing the results as the Kdofs are for the complete RVE, not the subdomain. The stresses and strains will have to be averaged one way or another, which involves the integration of stresses and strains over the subdomain. Since implementing the integration procedure of the stresses and strains as a computer code is a demanding job for most users, the following simplifications, which keep the mathematical rigour, can be helpful. One will find Green’s formula Z I I ny vfy vfz ds (9.11) þ dA ¼ fz dy þ fy dz ¼ ½ fy fz vy vz nz A
vA
vA
and the Gauss theorem (divergence theorem) 8 9 > ZZZ = < nx > vfx vfy vfz þ þ dU ¼ % ½ fx fy fz ny dS > vx vy vz ; : > vU nz U
(9.12)
useful, where ½ fx fy and ½ fx fy fz are continuously differentiable vector fields in 2D and 3D, A and U are the domains of interest in 2D and 3D spaces, respectively, with vA and vU being their boundaries, ½ ny nz and ½ nx ny nz are the unit outward normals to vA and vU, s is the arc length of vA and S the surface area of vU. In particular, if one introduces ½ fx fy ¼ ½ gx hx gy hy for 2D cases and ½ fx fy fz ¼ ½ gx hx gy hy gz hz for 3D cases, respectively, where g and h are two separate continuously
337
RVE for media with randomly distributed inclusions
differentiable vector fields like f, Green’s formula and the Gauss theorem give the following integration by parts, respectively Z I Z ny vgy vhy vgz vhz hy þ hz dA ¼ ½ fy fz þ gz ds gy dA vy vz vy vz nz vA
A
A
(9.13) ZZZ U
vgy vgx vgz hx þ hy þ hz dU ¼ vx vy vz
¼ % ½ gx hx vU
8 9 nx > > > > > ZZZ = < > vhy vhx vhz gz hz ny dS þ gy þ gz gx dU > > vx vy vz > > > > : ; U nz (9.14)
gy hy
With these mathematical formulae, the average strains and stresses in 2D case can be derived. Specifically, making use of Green’s formula (9.13) and referring to Fig. 9.12 for symbols involved, average strains are determined as Z Z Z I 1 1 vv 1 vv v0 1 0 εy ¼ εy dA ¼ 0dy þ vdz dA ¼ dA ¼ A A vy A vy vz A A
¼ ¼
1 A 1 A
A
I vdz ¼ vA Zz2
1 A
A
ZC3 vdz þ C2
vjy¼y2 dz
Zz2
1 A
z1
1 A
ZC1 vdz ¼ C4
1 A
vA ZC4
ZC3 vdz C2
(9.15)
z1
z z2 z1
C4
C1
y1
C3
C2
y2
vdz C1
vjy¼y1 dz;
b
1 A
y
a
Fig. 9.12 An RVE and its subdomain.
338
Representative Volume Elements and Unit Cells
1 ε0z ¼ A
Zy2 y1
g0yz ¼
1 A
Z gyz dA ¼ A
(9.16)
1B @ A
Z
ZC2
1B @ A 0
1 A
I vw vv 1 vdy þ wdz ¼ þ dA ¼ vy vz A vA
A
ZC3 vdy þ
C1
¼
wjz¼z1 dy; y1
0 ¼
Zy2
1 wjz¼z2 dy A
wdz C2
Zy2
ZC1 vdy þ
C3
Zy2 vjz¼z2 dy
y1
ZC4
C wdzA ¼
C4
Zz2 vjz¼z1 dy þ
y1
1
Zz2 wjy¼y2 dz
z1
1 C wjy¼y1 dzA;
z1
(9.17) where notation A is used both as the subdomain for post-processing and the area of it for averaging without confusion. To apply (9.15)e(9.17), it is essential that the sides of the subdomain are straight and parallel to the coordinate axes. To evaluate the integrations along the sides of the subdomain in the expressions above consistently, appropriate numerical integration rules should be employed for all elements involved along the side. If the elements are linear, i.e. 3- or 4-noded, the trapezium rule should be applied. For quadratic elements, one should use the Simpson’s rule. These are to reflect the fact that the contributions from different nodes are of different weights depending on the location of the node. It is clear that the coordinates of the nodes on all sides must be available and used appropriately as a part of the post-processing. Similar arguments can be made for average strains in 3D cases, for which the Gauss theorem should be used instead, resulting in the following expressions for average strains 8 9 > ZZZ = < nx > 1 vu v0 v0 1 0 εx ¼ þ þ dU ¼ % ½ u 0 0 ny dS > U vx vy vz U vU ; : > nz U 1 0 ZZ ZZ 1B C (9.18) udS udS A; ¼ @ U Sx2
Sx1
339
RVE for media with randomly distributed inclusions
ε0y ¼
1 U
ZZZ U
v0 vv v0 1 þ þ dU ¼ % ½ 0 v vx vy vz U vU
1 0 ZZ ZZ 1B C vdS vdSA; ¼ @ U Sy2
ε0z ¼
1 U
(9.19)
8 9 nx > > > = < > w ny dS > > > ; : > nz
(9.20)
Sy1
ZZZ U
8 9 nx > > > = < > 0 ny dS > > > ; : > nz
v0 v0 vw 1 þ þ dU ¼ % ½ 0 vx vy vz U vU
0
1 0 ZZ ZZ 1B C wdS wdS A; ¼ @ U Sz2
g0yz
1 ¼ U
Sz1
ZZZ U
v0 vw vv 1 þ þ dU ¼ % ½ 0 vx vy vz U vU
w
8 9 nx > > > > > > > < > = v ny dS > > > > > > > : > ; nz 1
0 ZZ ZZ ZZ ZZ 1B C wdS wdS þ vdS vdSA; ¼ @ U Sy2
Sy1
Sz2
Sz1
(9.21)
g0xz
1 ¼ U
ZZZ U
vw v0 vu 1 þ þ dU ¼ % ½ w vx vy vz U vU
0
8 9 nx > > > > > > > < > = u ny dS > > > > > > > : > ; nz 1
0 ZZ ZZ ZZ ZZ 1B C wdS wdS þ udS udS A; ¼ @ U Sx2
Sx1
Sz2
Sz1
(9.22)
340
g0xy
Representative Volume Elements and Unit Cells
1 ¼ U
ZZZ U
vv vu v0 1 þ þ dU ¼ % ½ v vx vy vz U vU
u
8 9 nx > > > > > > > = < > 0 ny dS > > > > > > > ; : > nz 1
0 ZZ ZZ ZZ ZZ 1B C vdS vdS þ udS udSA; ¼ @ U Sx2
Sx1
Sy2
Sy1
(9.23) where U has been used both as the subdomain and its volume and Sx1 , Sx2 , Sy1 , Sy2 , Sz1 and Sz2 are the six faces of the subdomain, which are flat and parallel to the coordinate planes as indicated in the subscripts. The numerical integrations on the faces of the subdomain will have to be evaluated consistently using a 2D numerical integration according to the tessellations on the faces and the order of the elements involved. It should be noted that users tend to base the calculations of the average stresses and strains on their intuition to obtain some formulae seemingly similar but not exactly the same as rigorously derived above. The point to make through the above derivation is to demonstrate that mathematical rigour is often readily available without introducing additional numerical demand. To convert the volume integration to surface integration for average stresses, similar process can be followed, but the outcomes tend to be slightly more complicated than what one would expect intuitively. Using integration by parts again along with necessary mathematical manipulations, one has ZZZ ZZZ 1 1 vx vx vx 0 sx ¼ sx dU ¼ sx þ sxy þ sxz dU U U vx vy vz U 8 U9 > ZZZ = < nx > 1 1 vsx vsxy vsxz ¼ % x½ sx sxy sxz ny dS þ þ xdU > U vU U vx vy vz ; : > nz U 0 0 1 ZZ ZZ ZZ 1B C 1B ¼ @x2 sx dS x1 sx dS A þ @ xsxy dS U U Sx2
Sx1
Sy2
0 1 ZZ ZZ ZZ C 1B C xsxy dS A þ @ xsxz dS xsxz dS A U Sy1
1
Sz2
Sz1
(9.24)
341
RVE for media with randomly distributed inclusions
vs
vsxz xy x where equilibrium equation vs vx þ vy þ vz ¼ 0 has been made use of. Similarly, 1 0 ZZZ ZZ ZZ 1 1B C s0y ¼ sy dU ¼ @ ysxy dS ysxy dS A U U U
Sx2
0 þ
1B @y2 U
ZZ
ZZ sy dS y1
Sy2
Sy1
1
Sx1
1 0 ZZ ZZ C 1B C sy dS A þ @ ysxz dS ysxz dS A U Sz2
Sz1
(9.25) 1 0 ZZ ZZ 1 1B C s0z ¼ sz dU ¼ @ zsxz dS zsxz dS A U U Sx2 Sx1 1 1 0U 0 ZZ ZZ ZZ ZZ 1B C 1B C þ @ zsyz dS zsyz dS A þ @z2 sz dS z1 sz dS A U U ZZZ
Sy2
Sy1
Sz2
Sz1
(9.26) s0yz
1 ¼ U
ZZZ U
1 syz dU ¼ U
1 ¼ % y½ sxz U vU
syz
ZZZ U
vy vy vy sxz þ syz þ sz dU vx vy vz
8 9 > > > nx > > > = > < sz ny dS > > > > > > > : nz ;
1 U
ZZZ U
vsxz vsyz vsz þ þ ydU vx vy vz
1 1 0 0 ZZ ZZ ZZ ZZ 1B C 1B C ¼ @ ysxz dS ysxz dS A þ @y2 syz dS y1 syz dS A U U Sx2
Sx1
Sy2
Sy1
0 1 ZZ ZZ 1B C þ @ ysz dS ysz dS A U Sz2
Sz1
(9.27)
342
Representative Volume Elements and Unit Cells
s0xz ¼
1 U
1 0 ZZ ZZ 1B C sxz dU ¼ @ zsx dS zsx dS A U
ZZZ U
Sx2
Sx1
1 1 0 0 ZZ ZZ ZZ ZZ 1B C 1B C þ @ zsxy dS zsxy dS A þ @z2 sxz dS z1 sxz dS A U U Sy2
Sy1
Sz2
Sz1
(9.28)
s0xy ¼
1 U
0
ZZZ sxy dU ¼ U
1B @x2 U
ZZ
ZZ sxy dS x1
Sx2
1 C sxy dS A
Sx1
1 1 0 0 ZZ ZZ ZZ ZZ 1B C 1B C xsy dS xsy dS A þ @ xsyz dS xsyz dS A þ @ U U Sy2
Sy1
Sz2
Sz1
(9.29) Apparently, the shear components of the average stress tensor can be derived differently resulting in different expressions as follows that should nevertheless produce the same numerical outcomes. 1 1 0 0 ZZ ZZ ZZ ZZ 1B C 1B C s0yz ¼ @ zsxy dS zsxy dS A þ @ zsy dS zsy dS A U U 0 1B þ @z2 U
Sx2
Szx
ZZ
ZZ syz dS z1
Sz2
1
Sy2
Sy1
C syz dS A
Sz1
(9.27a)
343
RVE for media with randomly distributed inclusions
0 s0xz ¼
1B 2 @x U
ZZ
ZZ sxz dS x1
Sx2
Sx1
1
1 0 ZZ ZZ C 1B C sxz dS A þ @ xsyz dS xsyz dS A U Sy2
1 0 ZZ ZZ 1B C þ @ xsz dS xsz dS A U Sz2
Sy1
Sz1
(9.28a) 1 1 0 0 ZZ ZZ ZZ ZZ 1B C 1B C s0xy ¼ @ ysx dS ysx dS A þ @y2 sxy dS y1 sxy dS A U U Sx2
Sx1
Sy2
1 0 ZZ ZZ 1B C þ @ ysxz dS ysxz dS A U Sz2
Sy1
Sz1
(9.29a) One can easily reduce the above from their 3D form to their 2D counterparts. 0 1 ZZ Zy2 Zy2 1 1B C s0x ¼ sx dA ¼ @x2 sx jx¼x2 dy x1 sx jx¼x1 dyA A U y1 y1 A (9.30) 0 x 1 Z2 Zx2 1 þ @ xsxy jy¼y2 dx ysxy jy¼y1 dxA U x1
s0y ¼
1 A
x1
0
ZZ sy dA ¼ A
0 þ
1@ y2 U
1B @ U
Zy2
Zy2 ysxy jx¼x2 dy
y1
Zx2
y1
Zx2 sy jy¼y2 dx y1
x1
x1
1
sy jy¼y1 dxA
1 C ysxy jx¼x1 dyA (9.31)
344
Representative Volume Elements and Unit Cells
s0xy ¼
1 A
0
ZZ sxy dA ¼
þ
1@ y2 A
Zy2 ysx jx¼x2 dy
y1
A
0
1B @ A
Zy2
y1
Zx2
Zx2 sxy jy¼y2 dx y1
x1
1 C ysx jx¼x1 dyA
1
(9.32)
sxy jy¼y1 dxA
x1
Equally, one can also have 1 0 Zy2 Zy2 1B C s0xy ¼ @x2 sxy jx¼x2 dy x1 sxy jx¼x1 dyA A 0 þ
1@ A
y1
y1
Zx2
Zx2 xsy jy¼y2 dx
x1
1
(9.32a)
xsy jy¼y1 dxA
x1
Relatively, the expressions of average stresses as surface integrations derived above are less likely to be obtained intuitively. For instance, intuition easily leads to the following two expressions as the average direct stress on each of the two surfaces perpendicular to the y-axis. ZZ ZZ y2 y1 y2 y1 0 0 sy ¼ sy dS or sy ¼ sy dS (9.33) U U Sy2
Sy1
If one has achieved complete convergence in terms of RVE size and absolute accuracy, there should not be any difference between either of (9.33) and mathematically derived expression (9.25) in this particular case. However, given the approximate nature of FEM, the accuracy of the derived expressions will be completely consistent with the FEM approximation whilst the intuitive ones lack such consistency. A point to be made here is that although the post-processing is not really exciting, it does not mean that it should be tampered casually, especially when consistent results can be obtained which are not too difficult to evaluate. FEM does not aim at exact solution. However, its accuracy erodes if it is not followed in a consistent manner. The consistency is underlined by the rationality, i.e. the strict definition of physical properties, mathematical rules and logic. For the integrity of science and engineering, it is worth one’s while to search for the rational solutions. Rational ones might demand a bit more skills and patience
RVE for media with randomly distributed inclusions
345
to obtain, but they are not necessarily more complicated to implement. The beauty of truth is that rational ones are often simpler. A good example on the same subject is the post-processing for UCs through the Kdofs as established in Chapter 7. Calling intuitive methods devised for obtaining the desired outcomes the ‘engineering approaches’ could sometimes be an insult to the name. Proper engineering approaches are taken by competent engineers to achieve approximate solutions in absence of rigorous theoretical solutions, whilst keeping necessary control of the errors introduced by the extra assumptions to deliver the approximations. Use of casual approaches that add unnecessary errors to the solutions where theoretically accurate solutions can be obtained can never be justified, especially if the casual approaches are not any simpler than their rigorous counterparts.
9.7 Conclusions It is the subject of this chapter how RVEs should be constructed and analysed and how the results should be extracted through appropriate postprocessing. Practices in this subject as found in the literature tend to be loosely presented without much mathematical rigour. However, through the elaborations and derivations, the coherence in the subject has been firmly established. Together with the conclusion from Chapter 4, an RVE can be constructed from the original configuration without falsifying periodicity. Uniform displacement boundary conditions are prescribed to generate a macroscopically uniaxial stress state or pure shear stress state. A subdomain of a number of characteristic lengths into the RVE should be extracted for post-processing in order to obtain average stresses and strains, from which effective material properties can be evaluated. The evaluation of average stresses and strains can be based on mathematically derived formulae in consistence with the FEM. Intuition is neither necessary nor helpful in this particular case.
References Abaqus Analysis User’s Guide, 2016. Abaqus 2016 HTML Documentation. Bulsara, V.N., Talreja, R., Qu, J., 1999. Damage initiation under transverse loading of unidirectional composites with arbitrarily distributed fibers. Composites Science and Technology 59, 673e682. Gusev, A.A., Hine, P.J., Ward, I.M., 2000. Fiber packing and elastic properties of a transversely random unidirectional glass/epoxy composite. Composites Science and Technology 60, 535e541. Hill, R., 1963. Elastic properties of reinforced solids: some theoretical principles. Journal of the Mechanics and Physics of Solids 11, 357e372.
346
Representative Volume Elements and Unit Cells
Soden, P.D., Hinton, M.J., Kaddour, A.S., 1998. Lamina properties, lay-up configurations and loading conditions for a range of fibre-reinforced composite laminates. Composites Science and Technology 58, 1011e1022. Wongsto, A., Li, S., 2005. Micromechanical FE analysis of UD fibre-reinforced composites with fibres distributed at random over the transverse cross-section. Composites Part A: Applied Science and Manufacturing 36, 1246e1266. Zou, Z., Li, S., 2002. Stresses in an infinite medium with two similar circular cylindrical inclusions. Acta Mechanica 156, 93e108.
RVE for media with randomly distributed inclusions 9.1 Introduction Modern materials often exhibit different internal architectures at different length scales or in different perspectives. Some of those architectures might be of completely irregular physical and geometric characteristics, such as that in particulate reinforced composites, in which reinforcing particulates are usually dispersed at random, or the random fibre distribution in the transverse cross-section of UD composites. Micromechanical modelling of such problems requires the use of representative volume elements (RVEs). The definition of an RVE and its associated considerations have been presented in Chapter 4, followed by a critical review of commonly observed misperceptions and misuse of RVEs in the literature in Chapter 5. Having reached the level of understanding thus far, the need is apparent for RVE formulation and analysis techniques that are logically derived based on the correct concepts and deliver results that are mathematically accurate. The task of this chapter is to offer such an approach to address the problem and elaborate on the steps taken in order to deliver the subject without ambiguity. Exact solution would be obtained if the problem could be analysed over an infinite domain, or the correct boundary conditions had been prescribed if a finite domain is employed as an RVE. However, none of these options is practical in general. The proposed approach is based on the principle of Saint-Venant. In the context of the present application, it states that the results obtained over the regions within the RVE sufficiently distant from the boundary, where inaccurate boundary conditions are prescribed, should be close enough to the exact solution, provided that the forces associated with such incorrect boundary conditions are statically equivalent to those appearing in the infinite domain or those in the exact boundary conditions. However, the principle of Saint-Venant itself does not provide any estimate of how large the distance from the boundary should be. In fact, it varies from problem to problem in general and the characteristic measure of this Representative Volume Elements and Unit Cells ISBN: 978-0-08-102638-0 https://doi.org/10.1016/B978-0-08-102638-0.00009-8
© 2020 Elsevier Ltd. All rights reserved.
319
j
320
Representative Volume Elements and Unit Cells
distance, or decay length, as defined in Chapter 4, will have to be determined in an ad hoc manner for a given type of problems.
9.2 Displacement boundary conditions and traction boundary conditions for an RVE The description of uniform displacement and uniform traction boundary conditions as referred to in Chapter 4 is in fact rather simplistic in a sense that it does not provide instructions of how such boundary conditions are to be implemented. The elaborated description of these boundary conditions will therefore be provided here to benefit potential users. Consider a rectangular area, 0 ya and 0 z b, as the volume element in the yz-plane in a 2D space for the ease of presentation, which can be easily extended to a cuboid in a 3D space. Assuming a uniaxial stress s0y is to be prescribed in the y-direction in order to determine the effective Young’s modulus and Poisson’s ratio in this direction, the traction boundary conditions can be given as follows. py y¼0 ¼ s0y and pz jy¼0 ¼ 0; and pz jy¼a ¼ 0; py y¼a ¼ s0y (9.1) and pz jz¼0 ¼ 0; py z¼0 ¼ 0 and pz jz¼b ¼ 0; py z¼b ¼ 0 where py and pz are the traction components. If one is interested in the determining the effective shear modulus, the traction boundary conditions for an RVE under a pure shear stress s0yz in the yz-plane are to be applied as and pz jy¼0 ¼ s0yz py y¼0 ¼ 0 py y¼a ¼ 0 and pz jy¼a ¼ s0yz (9.2) py z¼0 ¼ s0yz and pz jz¼0 ¼ 0 and pz jz¼b ¼ 0: py z¼b ¼ s0yz For FEM applications, these conditions will have to be topped up with appropriate constraints to rule out rigid body motions, in particular, rigid body rotations. The correct method of constraining rigid body rotations when the domain deviates from orthogonal regularity is elaborated in
321
RVE for media with randomly distributed inclusions
Section 5.5 of Chapter 5. Regarding the application of tractions, an important point is to be made about the vanishing traction components as involved in (9.1) and (9.2). Since traction boundary conditions are natural boundary conditions, they are prescribed when loads are applied. Loads of zero magnitude imply the absence of loads and therefore they do not need to be applied explicitly in an FE analysis. In other words, vanishing traction boundary conditions should be left unattended and this is strictly correct as far as FE modelling is concerned, whilst they will be satisfied approximately through energy minimisation in the same way as equilibrium conditions are satisfied, as well as the nontrivial boundary conditions. In (9.1) and (9.2), what needs to be applied is therefore reduced to the following. py y¼0 ¼ s0y and py y¼a ¼ s0y for uniaxial tension; and (9.3) pz jy¼0 ¼ s0yz py z¼0 ¼ s0yz
and and
pz jy¼a ¼ s0yz py z¼b ¼ s0yz
for pure shear.
(9.4)
Note that if one wishes to apply uniaxial stress s0z in the z-direction, boundary conditions for such case will be similar to those given by (9.3) but will have to be applied on surfaces z¼0 and z¼b. The so-called uniform displacement boundary conditions are not usually prescribed purely in terms of displacements. The boundary conditions for the same loading case as (9.1) and (9.2) are actually defined in a mixed form as follows. vjy¼0 ¼ 0
and
pz jy¼0 ¼ 0
and pz jy¼a ¼ 0 vjx¼a ¼ aε0y py z¼0 ¼ 0 and pz jz¼0 ¼ 0 py z¼b ¼ 0 and pz jz¼b ¼ 0 py y¼0 ¼ 0 and wjy¼0 ¼ 0 py y¼a ¼ 0
and
vjz¼0 ¼ 0
and pz jz¼0 ¼ 0
vjz¼b ¼ bg0yz
and
wjy¼a ¼ 0 pz jz¼b ¼ 0
for uniaxial tension; and
for pure shear.
(9.5)
(9.6)
322
Representative Volume Elements and Unit Cells
These boundary conditions should be topped up with the constraints to eliminate rigid body motions remaining in the system, viz. translation in the z-direction for (9.5) and rigid body rotation for (9.6). Ignoring the vanishing natural boundary conditions, one is left with the following, respectively. vjy¼0 ¼ 0 (
and
vjy¼a ¼ aε0y
for uniaxial tension; and
wjy¼0 ¼ wjy¼a ¼ 0 vjz¼0 ¼ 0
and
vjz¼b ¼ bg0yz
for pure shear
(9.7) (9.8)
One can apply uniaxial strain ε0z in a similar manner as ε0y is prescribed in (9.7). Note that the vanishing traction boundary conditions in (9.5) and (9.6) should not be replaced by vanishing displacements. Otherwise, they will correspond to a completely different loading conditions. It should also be pointed out that although (9.1) and (9.2) or (9.3) and (9.4) are unique forms of prescribing the relevant traction boundary conditions, the way of constraining rigid body motions as a necessary part of the imposition of boundary conditions can take different forms, and rigid body motion constraints have not been included in (9.1) and (9.2) or (9.3) and (9.4). On the other hand, (9.5) and (9.6) or (9.7) and (9.8) are not unique forms of prescribing displacement boundary conditions for the relevant effective stress states. This is due to the fact that any combination of rigid body motions can be superimposed on top of them to present a different appearance without changing the stress state, as has been explained in Section 6.2 of Chapter 6. For instance, one can have 1 vjy¼0 ¼ aε0y 2
and
1 vjy¼a ¼ aε0y 2
for uniaxial tension; and (9.9)
wjy¼0 ¼ 0
and
1 wjy¼a ¼ ag0yz 2
for pure shear (9.10) 1 0 vjz¼0 ¼ 0 and vjz¼b ¼ bgyz 2 They will be to exactly the same effects as (9.7) and (9.8). Any lack of uniqueness is a source of potential confusions and it is the responsibility of the users to apply the displacement boundary conditions correctly.
RVE for media with randomly distributed inclusions
323
In much the same way as in UCs, when (9.7) and (9.8) are employed, average strains ε0y and g0yz serve as Kdofs at which loads can be applied. The concept of Kdofs and their applications have been introduced in Section 6.6 of Chapter 6. Whilst displacement boundary conditions as provided above could be perfectly correct for UCs, as obtained from the symmetry conditions due to the regularity of the architecture at the lower length scale, the traction boundary conditions are usually incorrect unless the material is homogeneous even at its lower length scale. In absence of regularity as is usually the case for RVEs, even the displacement boundary conditions cannot be prescribed correctly, in general, let alone the traction boundary conditions. Any application of such boundary conditions will have to be considered as some kind of approximation. The extent of such an approximation determines the validity of the RVE introduced, as will be the subject of this chapter. Hill (1963) suggested that the volume element was representative if one could obtain sufficiently close results under both the uniform traction and the uniform displacement boundary conditions. As commented in Chapter 4, this is unnecessary, as will be elaborated fully in the next a few sections.
9.3 Decay length for boundary effects Without perfect symmetry, neither uniform tractions nor uniform displacements will give exact boundary conditions for the volume element concerned. However, either of them will represent a statically equivalent uniaxial loading condition if homogeneity is present at the upper length scale. According to the principle of Saint-Venant, either of them should result in perfectly correct stress and strain distributions an appropriate distance into the volume element from its boundary. This distance varies from point to point along the boundary. The maximum value of this distance for all relevant loading conditions is the decay length, that has already been introduced in Chapter 4. The inner part of the volume element a decay length away from its boundary in all directions forms a subdomain of the volume element. As long as the original volume element is large enough to allow such a subdomain to be also representative in terms of volume fractions of constituents, the subdomain will be an RVE, a perfect one. The effective properties of the material should not be extracted directly from the original RVE itself, but from the subdomain where the stress and strain distributions are free from the effects of incorrect boundary conditions, in
324
Representative Volume Elements and Unit Cells
terms of either traction or displacement. The question that remains is what should be the decay length. The principle of Saint-Venant does not provide any clear guideline for it, and in fact, the decay length varies from material to material and from one effective property to another, hence it will have to be determined by the user on a case-by-case basis. In order to estimate the decay length for UD composites, one with regularly packed fibres in its transverse cross-section is considered as an example. In a composite so chosen, symmetry is present and hence correct boundary conditions can be obtained if the border of volume element is selected along symmetry lines as a stack of repeating cells. In order to examine the effect of inaccurate boundary conditions, one side of the stack of such cells is deliberately truncated away from the symmetry line as shown in Fig. 9.1, such that no correct boundary conditions can be prescribed on it. Two analyses were conducted with each stack of cells, one involving uniform traction boundary conditions, as defined by (9.3), whilst another was carried out with uniform displacement (9.7) being prescribed. Both types of boundary conditions are statically equivalent to uniaxial tension. Comparing the stress distributions from the two analyses in the truncated cells, one can observe how far the disturbance introduced by the incorrect boundary conditions propagates as a way of evaluating the decay length. The results shown in Figs. 9.2 and Fig. 9.3 are von Mises stress contours plotted over a squarely packed array and a hexagonally packed array, respectively, where Fig. 9.2(A) and Fig. 9.3(A) were obtained with uniform displacement being applied, whilst in results shown in Fig. 9.2(B) and Fig. 9.3(B) were obtained with uniform traction being prescribed. The von Mises stress is not usually relevant for composites in general. However, in this case, both the fibre and matrix are assumed to be isotropic and the von Mises stress should be meaningful in the context of the present discussion. Only complete cells have been shown in Figs. 9.2 and 9.3, since the stress distributions in the truncated ones are deemed to be wrong. In Figs. 9.2 and 9.3, the columns with title ‘Case-0’ show stress contour over the stacks that did not involve the truncated cells hence had correct boundary conditions applied to them. Such stacks were analysed in order to compare the results with those obtained for truncated stacks. The remaining columns were obtained under incorrect boundary conditions, as each of those stacks had on the top an incomplete cell truncated at different locations as shown in Fig. 9.1. The pattern of contour plots in the column of Case-0 in Fig. 9.2(A) is identical for every cell in the stack, because Case-0 represents a correct solution to the problem. Comparing the patterns in Case-1, Case-2 and Case-
RVE for media with randomly distributed inclusions
325
Fig. 9.1 Arrays of periodic cell of (A) square packing (B) hexagonal packing.
3 with that in Case-0, it can be observed that the further away the cells are from the top border where incorrect boundary conditions are prescribed, i.e. moving from cell No. 1 to cell No. 5, the closer the comparison of the stress contours with that in Case-0. In fact, there is no noticeable difference
326
Representative Volume Elements and Unit Cells
(a) L1 1→
(b) Case-0 No.1
Case-1
Casse-2
Case-3
Case-0 C
L1 →
Case-1
Case-2
Case-3
No.1
L2 2→
No.2 2
L2 →
No.2
L3 →
No.3 3
L3 →
No.3
L4 L →
No.4 4
L4 →
No.4
L5 5→
No.5 5
L5 → No.5
Fig. 9.2 Von Mises stress contours over complete periodic cells from the square packed arrays (A) under prescribed uniform displacement and (B) under prescribed uniform traction.
between stress contours from cell No.2 onwards. Comparing respective cases in Fig. 9.2(A) and (b), it becomes obvious that relatively, the boundary effects decay over a longer distance if a uniform traction is prescribed, although disturbances in the latter case settle down eventually by cell No. 4 or so. The same trends are repeated in Fig. 9.3 for hexagonal packing where the plots have been arranged in exactly the same order as those in Fig. 9.2. In other words, under uniform displacement the decay length is approximately equal to the length of one cell, whilst under uniform traction it can be several cells. Apparently, the boundary conditions corresponding to prescribed uniform displacement represent a better approximation to the exact boundary conditions than those of prescribed uniform traction. Therefore, it is advisable that the displacement boundary condition should be employed. In this case, the typical decay length will be a couple of characteristic lengths. More justifications can be found in (Wongsto and Li, 2005). For regular patterns, the characteristic length is the cell size. In the case of
327
RVE for media with randomly distributed inclusions
(a)
(b) Case-0
L1 o
Case-11
Casse-2
Case-3
L1 o
Case-0
Case-1
C ase-2
Case-3
No o.1
No.1
L2 o
L2 o
No o.2
No.2
L3 o
L3 o
No o.3
No.3
L4 o
L4 o
No o.4
No.4
L5 o
L5 o
No o.5
No.5
Fig. 9.3 Von Mises stress contours over complete periodic cells from the hexagonally packed arrays (A) under prescribed uniform displacement and (B) under prescribed uniform traction. Generating the microstructures of randomly distributed physical and geometric features.
irregular patterns, it can be expected to be the average distance between the centres of two adjacent fibres. This will be further verified in Section 9.5.
9.4 Generation of random patterns In order to facilitate theoretical construction of RVEs, a means of generating a randomly distributed features is desirable. A UD composite having random fibre distribution over its transverse cross-section will be taken as a special case in a 2D space. The process of reproducing the geometry of cross-section is equivalent to generating randomly spaced, nonoverlapping circular discs in a plane within a domain. There could be different approaches to achieve this goal. Digitising microscopic photographs such as that shown in Fig. 4.1 is an obvious one, but the procedure could be computationally costly and time-consuming. Employing a suitable numerical algorithm would be more effective in this respect. One possibility is to adopt a ‘coin-dropping’, or random insertion process, as was
328
Representative Volume Elements and Unit Cells
implemented by Bulsara et al. (1999). Whilst the random nature of the obtained distribution is obvious, it may be difficult to achieve a specified disc (constituent) volume fraction, especially at high volume fractions. Constituent volume fraction is often a dominant parameter in determining the performance of the composite and a reasonable control of it is essential. A different process was employed by Gusev et al. (2000) using a Monte Carlo method. An alternative scheme was proposed in Wongsto and Li (2005), which is capable of reproducing any practical volume fraction and yet maintaining the random nature in distribution. The outline of the procedure is as follows. A set of regularly packed discs of radius R in the yz-plane is generated. An example of such set with hexagonally packed discs is shown in Fig. 9.4. The spacing between discs, 2b, can be calculated making use of the constituent volume fraction, Vf, as 2 pR ffiffi . Next, two frames are introduced and fixed in the plane, Frame b2 ¼ 2p 3V f
1 and Frame 2, with the latter being placed inside the former. The distance from the border of Frame 2 to that of Frame 1 should be equal to a decay length, which, according to the conclusion reached in the previous section, is a few characteristic lengths, b. As in Fig. 9.4, it is approximately 4b. These frames do not have to be perfect squares, but they both have to contain
Frame-1
Frame-2
Fig. 9.4 Hexagonally packed discs at Vf ¼ 65% involving 105 discs.
RVE for media with randomly distributed inclusions
329
sufficient number of discs. They have been chosen to be square, i.e. b¼a, as shown in Fig. 9.4 for the sake of argument. As a next step, each disc is moved in a random manner as illustrated in Fig. 9.5. A random angle q between 0 and 360 is generated first to determine the direction for the disc to shift. Along this direction, a maximum distance of shift can be found, which is the smaller of the distance to the border of Frame 1 and that to the point when the disc comes into contact with another, marked as r in Fig. 9.5. Then a second random number, k, is generated between 0 and 1 and the actual distance of shift is given as kr. This process is performed for every disc to complete one iteration. One iteration is effectively one shake of the layout. A sufficiently large number of iterations will have to be made, e.g. 250, before the original regularity disappears. After that, a check on the volume fractions within both frames will be conducted after each iteration. The iterations terminate when the volume fractions in both frames are close enough to the original one, allowing a small discrepancy, e.g. 2% error, as the criterion for terminating iterations. The procedure can be easily programmed and, given parameters of R, Vf and the number of discs to be involved, a set of coordinates can be obtained as the centres of the randomly distributed discs. As examples of the applications of the scheme outlined above, three cases have been generated at a volume fraction of 65% as shown in Fig. 9.6(AeC). They are of the similar volume fraction to that shown in Figure 4.1 and share all the features present in Figure 4.1. A further case was generated at Vf ¼ 60% as shown in Fig. 9.6(D), for which some experimental data are available (Soden et al., 1998). Computationally, this procedure is very efficient, and any of the cases can be generated literally in a matter of seconds on an ordinary PC. Cases as shown in Fig. 9.6 were meshed and analysed with the results being presented in the next section.
U kU
T
z y
Fig. 9.5 Schematic diagram for perturbation process.
330
Representative Volume Elements and Unit Cells
(a) Vf1 = 65.69%, and Vf2 =64.17%
(b) Vf1 = 65.84% and Vf2 = 64.26%
(c) Vf1 = 66.38% and Vf2 = 64.59%
(d) Vf1 = 61.00% and Vf2 = 59.40%
Fig. 9.6 Random fibre distribution obtained for a case with original volume fraction Vf ¼ 65% after (A) 265 iterations; (B) 286 iterations and (C) 382 iterations. Fibre distribution for a case with Vf ¼ 60% after 262 iterations e plot (D). Volume fractions within Frame 1 and Frame 2 are denoted as V1f and V2f , respectively.
9.5 Strain and stress fields in the RVE and the subdomain The cases as generated in the previous section represent transverse cross-sections of UD composites when randomly distributed discs are viewed as fibres and the surrounding area as the matrix. To proceed with the subsequent micromechanical FE analysis, area delimited by Frame 1 will be treated as the RVE. For the purpose of extracting and processing
RVE for media with randomly distributed inclusions
331
the results afterwards, the areas inside and outside Frame 2 are meshed individually, whilst the compatibility of the mesh will be maintained across the border of Frame 2 so that the part of the mesh inside Frame 2 can be extracted as a subdomain. Examples of these meshes are shown in Figs. 9.7 and 9.8, where respective RVEs correspond to random distributions presented in Fig. 9.6(A) and (D). In Fig. 9.7, a further frame was introduced within Frame 2. This does not compromise the validity of the mesh except that it had an additional consideration during meshing. The purpose of it will be clear when the results are discussed. Any gap between two neighbouring fibres will be meshed into at least two elements across the gap for representativeness of the mesh. In the unlikely event of two touching fibres, compromises will have to be made using triangular elements of very sharp corners. The analysis is carried out based on the generalised plane strain idealisation as has been described in Chapter 6 in the context of UC applications. The quadratic triangular and quadrilateral generalised plane strain elements CPEG6 and CPEG8 of Abaqus (2016) were used for the longitudinal/transverse tension and transverse shear deformation whilst heat conduction elements DC2D6 and DC2D8 were adapted for the anticlastic problem of longitudinal shear based on an analogy as described in Section 6.4.1.1 of Chapter 6.
Fig. 9.7 The mesh for the RVE with fibres distributed at random for the case of Vf¼65%.
332
Representative Volume Elements and Unit Cells
Fig. 9.8 The meshes for one of the RVEs with fibres distributed at random at Vf¼60%.
It has been established in Section 9.3 that imposition of a uniform normal displacement is a more advantageous method of prescribing the loading to the RVE. Therefore, to achieve a macroscopically uniaxial stress state in the y-direction, displacements as defined in (9.7) were applied along the relevant sides of the boundary. This can be conveniently achieved by making use of the Kdof ε0y . As has been established in Chapter 7, a uniform displacement boundary condition can be prescribed by applying a concentrated ‘force’ As0y , where A¼a2 is the planar area of the RVE, to the Kdof instead of a prescribed nodal ‘displacement’ without changing the nature of boundary conditions, provided that displacements on the boundary have been associated with the Kdofs. This ensures a uniform displacement imposed on the boundary whilst leaving its magnitude for the subsequent analysis to determine as the ‘nodal displacement’ at the Kdof concerned. To facilitate the analysis, for all the cases with 65% fibre volume fraction, both the fibre and matrix are assumed to be isotropic and homogeneous with elastic properties Em ¼ 1 GPa, nm ¼ 0.3 and Ef ¼ 10 GPa and nf ¼ 0.2 for matrix and fibre, respectively. For the case with 60% fibre volume fraction, material properties of both fibre and matrix are listed in Table 9.1 (Soden et al., 1998).
333
RVE for media with randomly distributed inclusions
Table 9.1 Constituents properties of the Silenka Eglass 1200 tex and epoxy UD composite.
Matrix Young’s modulus Matrix Poisson’s ratio Fibre Young’s modulus Fibre Poisson’s ratio Fibre volume fraction
3.35 (GPa) 0.35 74 (GPa) 0.2 60%
The analyses were carried out on the domain Frame 1. To obtain the pattern of deformation and stress distribution in the RVE, a load as the ‘concentrated force’ equivalent to 1 MPa average uniaxial stress in the y-direction was applied to it. The results extracted from the subdomain Frame 2 were considered to represent the correct responses of the material, free from the effects of incorrectly prescribed boundary conditions along the border of Frame 1. The analysis and data extraction process have been applied to all the four cases shown in Fig. 9.6. The results are presented and discussed below. The von Mises stress contour plot over the subdomain Frame 2 corresponding to the case with Vf¼65% as shown in Fig. 9.6(A) is displayed in Fig. 9.9. The contour plot is shown in a deformed configuration with the magnitude of deformation being amplified by a factor of 300. The areas with the highest values of stress are always found at interfaces, especially
Fig. 9.9 Contours of von Mises stress in Frame 2 as a part of RVE for Case-1 under macroscopically uniaxial transverse tension.
334
Representative Volume Elements and Unit Cells
in areas where fibres are close to each other in the direction of loading (Zou and Li, 2002). The distortion of the shape of fibre cross-section is small relative to overall deformation because of fibre stiffness being higher than that of the matrix. Therefore, the deformation is mostly taken by the matrix although stresses in the fibres tend to show higher levels than those in the matrix. This explains why observed stress concentrations are so located. As a result of deformation, the edges of Frame 2 do not remain straight. This provides an indication how wrong it would be if Frame 2 was analysed alone and straight edges had been assumed when prescribing displacement boundary conditions. To further emphasise this, a case was generated where the domain Frame 2 was analysed directly, with uniform displacements prescribed along the straight edges of this domain. The stress contour for this case is shown in Fig. 9.10 with the same deformation magnification factor of 300. Comparing this case with results obtained from the more accurate analysis conducted on Frame 1 from which Frame 2 was extracted, as shown in Fig. 9.9, significant difference can be found. The maximum von Mises stress in the incorrect analysis were found to be 9.879 MPa in the fibres and 4.501 MPa in the matrix whilst their counterparts obtained in the more accurate analysis were 14.371 MPa and 10.997 MPa, respectively. A potential but significant risk of employing incorrect analysis is that if one uses the von Mises criterion with a given yield stress for the matrix, the load for the onset of plastic deformation in the matrix will be miscalculated.
Fig. 9.10 Von Mises stress contour over the Frame 2 alone analysed with inaccurate boundary conditions.
RVE for media with randomly distributed inclusions
335
In this particular case, for instance, the respective load was underestimated by approximately a factor of 2.5. Given the erroneous analysis as shown in Fig. 9.10, reasonable accuracy can still be achieved in the subdomain within Frame 2 in Fig. 9.7 as mentioned before. To demonstrate this, the von Mises stress contour plots for this same subdomain but from different analyses, one on Frame 1 and the other on Frame 2, are shown in Fig. 9.11. Both plots share many common features, although the distance to the borders of Frame 2 is not quite equal to a decay length and the innermost frame may not be large enough to be representative. Some of the differences in the patterns of the contour
Fig. 9.11 Frame 3 extracted from (a) Fig. 9.10 and (b) Fig. 9.11, both with a deformation magnification factor of 300.
336
Representative Volume Elements and Unit Cells
plots are due to the difference in the scales of the stress levels as shown in the legends. The consistency in their trends reinforces the justification made to the proposed methodology based on the fact that under approximate boundary conditions, correct solutions can be obtained from an inner subdomain at least one decay length away from the boundary where the approximate boundary conditions are prescribed. The underlying mechanical principle is that of Saint-Venant.
9.6 Post-processing for average stresses, strains and effective properties As was mentioned in Chapter 5, a proper description of the procedure for post-processing the results from the analyses of RVEs and UCs has been a topic often avoided the literature. In the case of RVEs, if stresses and strains are to be extracted from a subdomain, one cannot make use of Kdofs when post-processing the results as the Kdofs are for the complete RVE, not the subdomain. The stresses and strains will have to be averaged one way or another, which involves the integration of stresses and strains over the subdomain. Since implementing the integration procedure of the stresses and strains as a computer code is a demanding job for most users, the following simplifications, which keep the mathematical rigour, can be helpful. One will find Green’s formula Z I I ny vfy vfz ds (9.11) þ dA ¼ fz dy þ fy dz ¼ ½ fy fz vy vz nz A
vA
vA
and the Gauss theorem (divergence theorem) 8 9 > ZZZ = < nx > vfx vfy vfz þ þ dU ¼ % ½ fx fy fz ny dS > vx vy vz ; : > vU nz U
(9.12)
useful, where ½ fx fy and ½ fx fy fz are continuously differentiable vector fields in 2D and 3D, A and U are the domains of interest in 2D and 3D spaces, respectively, with vA and vU being their boundaries, ½ ny nz and ½ nx ny nz are the unit outward normals to vA and vU, s is the arc length of vA and S the surface area of vU. In particular, if one introduces ½ fx fy ¼ ½ gx hx gy hy for 2D cases and ½ fx fy fz ¼ ½ gx hx gy hy gz hz for 3D cases, respectively, where g and h are two separate continuously
337
RVE for media with randomly distributed inclusions
differentiable vector fields like f, Green’s formula and the Gauss theorem give the following integration by parts, respectively Z I Z ny vgy vhy vgz vhz hy þ hz dA ¼ ½ fy fz þ gz ds gy dA vy vz vy vz nz vA
A
A
(9.13) ZZZ U
vgy vgx vgz hx þ hy þ hz dU ¼ vx vy vz
¼ % ½ gx hx vU
8 9 nx > > > > > ZZZ = < > vhy vhx vhz gz hz ny dS þ gy þ gz gx dU > > vx vy vz > > > > : ; U nz (9.14)
gy hy
With these mathematical formulae, the average strains and stresses in 2D case can be derived. Specifically, making use of Green’s formula (9.13) and referring to Fig. 9.12 for symbols involved, average strains are determined as Z Z Z I 1 1 vv 1 vv v0 1 0 εy ¼ εy dA ¼ 0dy þ vdz dA ¼ dA ¼ A A vy A vy vz A A
¼ ¼
1 A 1 A
A
I vdz ¼ vA Zz2
1 A
A
ZC3 vdz þ C2
vjy¼y2 dz
Zz2
1 A
z1
1 A
ZC1 vdz ¼ C4
1 A
vA ZC4
ZC3 vdz C2
(9.15)
z1
z z2 z1
C4
C1
y1
C3
C2
y2
vdz C1
vjy¼y1 dz;
b
1 A
y
a
Fig. 9.12 An RVE and its subdomain.
338
Representative Volume Elements and Unit Cells
1 ε0z ¼ A
Zy2 y1
g0yz ¼
1 A
Z gyz dA ¼ A
(9.16)
1B @ A
Z
ZC2
1B @ A 0
1 A
I vw vv 1 vdy þ wdz ¼ þ dA ¼ vy vz A vA
A
ZC3 vdy þ
C1
¼
wjz¼z1 dy; y1
0 ¼
Zy2
1 wjz¼z2 dy A
wdz C2
Zy2
ZC1 vdy þ
C3
Zy2 vjz¼z2 dy
y1
ZC4
C wdzA ¼
C4
Zz2 vjz¼z1 dy þ
y1
1
Zz2 wjy¼y2 dz
z1
1 C wjy¼y1 dzA;
z1
(9.17) where notation A is used both as the subdomain for post-processing and the area of it for averaging without confusion. To apply (9.15)e(9.17), it is essential that the sides of the subdomain are straight and parallel to the coordinate axes. To evaluate the integrations along the sides of the subdomain in the expressions above consistently, appropriate numerical integration rules should be employed for all elements involved along the side. If the elements are linear, i.e. 3- or 4-noded, the trapezium rule should be applied. For quadratic elements, one should use the Simpson’s rule. These are to reflect the fact that the contributions from different nodes are of different weights depending on the location of the node. It is clear that the coordinates of the nodes on all sides must be available and used appropriately as a part of the post-processing. Similar arguments can be made for average strains in 3D cases, for which the Gauss theorem should be used instead, resulting in the following expressions for average strains 8 9 > ZZZ = < nx > 1 vu v0 v0 1 0 εx ¼ þ þ dU ¼ % ½ u 0 0 ny dS > U vx vy vz U vU ; : > nz U 1 0 ZZ ZZ 1B C (9.18) udS udS A; ¼ @ U Sx2
Sx1
339
RVE for media with randomly distributed inclusions
ε0y ¼
1 U
ZZZ U
v0 vv v0 1 þ þ dU ¼ % ½ 0 v vx vy vz U vU
1 0 ZZ ZZ 1B C vdS vdSA; ¼ @ U Sy2
ε0z ¼
1 U
(9.19)
8 9 nx > > > = < > w ny dS > > > ; : > nz
(9.20)
Sy1
ZZZ U
8 9 nx > > > = < > 0 ny dS > > > ; : > nz
v0 v0 vw 1 þ þ dU ¼ % ½ 0 vx vy vz U vU
0
1 0 ZZ ZZ 1B C wdS wdS A; ¼ @ U Sz2
g0yz
1 ¼ U
Sz1
ZZZ U
v0 vw vv 1 þ þ dU ¼ % ½ 0 vx vy vz U vU
w
8 9 nx > > > > > > > < > = v ny dS > > > > > > > : > ; nz 1
0 ZZ ZZ ZZ ZZ 1B C wdS wdS þ vdS vdSA; ¼ @ U Sy2
Sy1
Sz2
Sz1
(9.21)
g0xz
1 ¼ U
ZZZ U
vw v0 vu 1 þ þ dU ¼ % ½ w vx vy vz U vU
0
8 9 nx > > > > > > > < > = u ny dS > > > > > > > : > ; nz 1
0 ZZ ZZ ZZ ZZ 1B C wdS wdS þ udS udS A; ¼ @ U Sx2
Sx1
Sz2
Sz1
(9.22)
340
g0xy
Representative Volume Elements and Unit Cells
1 ¼ U
ZZZ U
vv vu v0 1 þ þ dU ¼ % ½ v vx vy vz U vU
u
8 9 nx > > > > > > > = < > 0 ny dS > > > > > > > ; : > nz 1
0 ZZ ZZ ZZ ZZ 1B C vdS vdS þ udS udSA; ¼ @ U Sx2
Sx1
Sy2
Sy1
(9.23) where U has been used both as the subdomain and its volume and Sx1 , Sx2 , Sy1 , Sy2 , Sz1 and Sz2 are the six faces of the subdomain, which are flat and parallel to the coordinate planes as indicated in the subscripts. The numerical integrations on the faces of the subdomain will have to be evaluated consistently using a 2D numerical integration according to the tessellations on the faces and the order of the elements involved. It should be noted that users tend to base the calculations of the average stresses and strains on their intuition to obtain some formulae seemingly similar but not exactly the same as rigorously derived above. The point to make through the above derivation is to demonstrate that mathematical rigour is often readily available without introducing additional numerical demand. To convert the volume integration to surface integration for average stresses, similar process can be followed, but the outcomes tend to be slightly more complicated than what one would expect intuitively. Using integration by parts again along with necessary mathematical manipulations, one has ZZZ ZZZ 1 1 vx vx vx 0 sx ¼ sx dU ¼ sx þ sxy þ sxz dU U U vx vy vz U 8 U9 > ZZZ = < nx > 1 1 vsx vsxy vsxz ¼ % x½ sx sxy sxz ny dS þ þ xdU > U vU U vx vy vz ; : > nz U 0 0 1 ZZ ZZ ZZ 1B C 1B ¼ @x2 sx dS x1 sx dS A þ @ xsxy dS U U Sx2
Sx1
Sy2
0 1 ZZ ZZ ZZ C 1B C xsxy dS A þ @ xsxz dS xsxz dS A U Sy1
1
Sz2
Sz1
(9.24)
341
RVE for media with randomly distributed inclusions
vs
vsxz xy x where equilibrium equation vs vx þ vy þ vz ¼ 0 has been made use of. Similarly, 1 0 ZZZ ZZ ZZ 1 1B C s0y ¼ sy dU ¼ @ ysxy dS ysxy dS A U U U
Sx2
0 þ
1B @y2 U
ZZ
ZZ sy dS y1
Sy2
Sy1
1
Sx1
1 0 ZZ ZZ C 1B C sy dS A þ @ ysxz dS ysxz dS A U Sz2
Sz1
(9.25) 1 0 ZZ ZZ 1 1B C s0z ¼ sz dU ¼ @ zsxz dS zsxz dS A U U Sx2 Sx1 1 1 0U 0 ZZ ZZ ZZ ZZ 1B C 1B C þ @ zsyz dS zsyz dS A þ @z2 sz dS z1 sz dS A U U ZZZ
Sy2
Sy1
Sz2
Sz1
(9.26) s0yz
1 ¼ U
ZZZ U
1 syz dU ¼ U
1 ¼ % y½ sxz U vU
syz
ZZZ U
vy vy vy sxz þ syz þ sz dU vx vy vz
8 9 > > > nx > > > = > < sz ny dS > > > > > > > : nz ;
1 U
ZZZ U
vsxz vsyz vsz þ þ ydU vx vy vz
1 1 0 0 ZZ ZZ ZZ ZZ 1B C 1B C ¼ @ ysxz dS ysxz dS A þ @y2 syz dS y1 syz dS A U U Sx2
Sx1
Sy2
Sy1
0 1 ZZ ZZ 1B C þ @ ysz dS ysz dS A U Sz2
Sz1
(9.27)
342
Representative Volume Elements and Unit Cells
s0xz ¼
1 U
1 0 ZZ ZZ 1B C sxz dU ¼ @ zsx dS zsx dS A U
ZZZ U
Sx2
Sx1
1 1 0 0 ZZ ZZ ZZ ZZ 1B C 1B C þ @ zsxy dS zsxy dS A þ @z2 sxz dS z1 sxz dS A U U Sy2
Sy1
Sz2
Sz1
(9.28)
s0xy ¼
1 U
0
ZZZ sxy dU ¼ U
1B @x2 U
ZZ
ZZ sxy dS x1
Sx2
1 C sxy dS A
Sx1
1 1 0 0 ZZ ZZ ZZ ZZ 1B C 1B C xsy dS xsy dS A þ @ xsyz dS xsyz dS A þ @ U U Sy2
Sy1
Sz2
Sz1
(9.29) Apparently, the shear components of the average stress tensor can be derived differently resulting in different expressions as follows that should nevertheless produce the same numerical outcomes. 1 1 0 0 ZZ ZZ ZZ ZZ 1B C 1B C s0yz ¼ @ zsxy dS zsxy dS A þ @ zsy dS zsy dS A U U 0 1B þ @z2 U
Sx2
Szx
ZZ
ZZ syz dS z1
Sz2
1
Sy2
Sy1
C syz dS A
Sz1
(9.27a)
343
RVE for media with randomly distributed inclusions
0 s0xz ¼
1B 2 @x U
ZZ
ZZ sxz dS x1
Sx2
Sx1
1
1 0 ZZ ZZ C 1B C sxz dS A þ @ xsyz dS xsyz dS A U Sy2
1 0 ZZ ZZ 1B C þ @ xsz dS xsz dS A U Sz2
Sy1
Sz1
(9.28a) 1 1 0 0 ZZ ZZ ZZ ZZ 1B C 1B C s0xy ¼ @ ysx dS ysx dS A þ @y2 sxy dS y1 sxy dS A U U Sx2
Sx1
Sy2
1 0 ZZ ZZ 1B C þ @ ysxz dS ysxz dS A U Sz2
Sy1
Sz1
(9.29a) One can easily reduce the above from their 3D form to their 2D counterparts. 0 1 ZZ Zy2 Zy2 1 1B C s0x ¼ sx dA ¼ @x2 sx jx¼x2 dy x1 sx jx¼x1 dyA A U y1 y1 A (9.30) 0 x 1 Z2 Zx2 1 þ @ xsxy jy¼y2 dx ysxy jy¼y1 dxA U x1
s0y ¼
1 A
x1
0
ZZ sy dA ¼ A
0 þ
1@ y2 U
1B @ U
Zy2
Zy2 ysxy jx¼x2 dy
y1
Zx2
y1
Zx2 sy jy¼y2 dx y1
x1
x1
1
sy jy¼y1 dxA
1 C ysxy jx¼x1 dyA (9.31)
344
Representative Volume Elements and Unit Cells
s0xy ¼
1 A
0
ZZ sxy dA ¼
þ
1@ y2 A
Zy2 ysx jx¼x2 dy
y1
A
0
1B @ A
Zy2
y1
Zx2
Zx2 sxy jy¼y2 dx y1
x1
1 C ysx jx¼x1 dyA
1
(9.32)
sxy jy¼y1 dxA
x1
Equally, one can also have 1 0 Zy2 Zy2 1B C s0xy ¼ @x2 sxy jx¼x2 dy x1 sxy jx¼x1 dyA A 0 þ
1@ A
y1
y1
Zx2
Zx2 xsy jy¼y2 dx
x1
1
(9.32a)
xsy jy¼y1 dxA
x1
Relatively, the expressions of average stresses as surface integrations derived above are less likely to be obtained intuitively. For instance, intuition easily leads to the following two expressions as the average direct stress on each of the two surfaces perpendicular to the y-axis. ZZ ZZ y2 y1 y2 y1 0 0 sy ¼ sy dS or sy ¼ sy dS (9.33) U U Sy2
Sy1
If one has achieved complete convergence in terms of RVE size and absolute accuracy, there should not be any difference between either of (9.33) and mathematically derived expression (9.25) in this particular case. However, given the approximate nature of FEM, the accuracy of the derived expressions will be completely consistent with the FEM approximation whilst the intuitive ones lack such consistency. A point to be made here is that although the post-processing is not really exciting, it does not mean that it should be tampered casually, especially when consistent results can be obtained which are not too difficult to evaluate. FEM does not aim at exact solution. However, its accuracy erodes if it is not followed in a consistent manner. The consistency is underlined by the rationality, i.e. the strict definition of physical properties, mathematical rules and logic. For the integrity of science and engineering, it is worth one’s while to search for the rational solutions. Rational ones might demand a bit more skills and patience
RVE for media with randomly distributed inclusions
345
to obtain, but they are not necessarily more complicated to implement. The beauty of truth is that rational ones are often simpler. A good example on the same subject is the post-processing for UCs through the Kdofs as established in Chapter 7. Calling intuitive methods devised for obtaining the desired outcomes the ‘engineering approaches’ could sometimes be an insult to the name. Proper engineering approaches are taken by competent engineers to achieve approximate solutions in absence of rigorous theoretical solutions, whilst keeping necessary control of the errors introduced by the extra assumptions to deliver the approximations. Use of casual approaches that add unnecessary errors to the solutions where theoretically accurate solutions can be obtained can never be justified, especially if the casual approaches are not any simpler than their rigorous counterparts.
9.7 Conclusions It is the subject of this chapter how RVEs should be constructed and analysed and how the results should be extracted through appropriate postprocessing. Practices in this subject as found in the literature tend to be loosely presented without much mathematical rigour. However, through the elaborations and derivations, the coherence in the subject has been firmly established. Together with the conclusion from Chapter 4, an RVE can be constructed from the original configuration without falsifying periodicity. Uniform displacement boundary conditions are prescribed to generate a macroscopically uniaxial stress state or pure shear stress state. A subdomain of a number of characteristic lengths into the RVE should be extracted for post-processing in order to obtain average stresses and strains, from which effective material properties can be evaluated. The evaluation of average stresses and strains can be based on mathematically derived formulae in consistence with the FEM. Intuition is neither necessary nor helpful in this particular case.
References Abaqus Analysis User’s Guide, 2016. Abaqus 2016 HTML Documentation. Bulsara, V.N., Talreja, R., Qu, J., 1999. Damage initiation under transverse loading of unidirectional composites with arbitrarily distributed fibers. Composites Science and Technology 59, 673e682. Gusev, A.A., Hine, P.J., Ward, I.M., 2000. Fiber packing and elastic properties of a transversely random unidirectional glass/epoxy composite. Composites Science and Technology 60, 535e541. Hill, R., 1963. Elastic properties of reinforced solids: some theoretical principles. Journal of the Mechanics and Physics of Solids 11, 357e372.
346
Representative Volume Elements and Unit Cells
Soden, P.D., Hinton, M.J., Kaddour, A.S., 1998. Lamina properties, lay-up configurations and loading conditions for a range of fibre-reinforced composite laminates. Composites Science and Technology 58, 1011e1022. Wongsto, A., Li, S., 2005. Micromechanical FE analysis of UD fibre-reinforced composites with fibres distributed at random over the transverse cross-section. Composites Part A: Applied Science and Manufacturing 36, 1246e1266. Zou, Z., Li, S., 2002. Stresses in an infinite medium with two similar circular cylindrical inclusions. Acta Mechanica 156, 93e108.