- Email: [email protected]
Apparent elastic properties of random fiber networks
Computational Materials Science 79 (2013) 715–723
Contents lists available at ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
Apparent elastic properties of random fiber networks Yikhan Lee, Iwona Jasiuk ⇑ Department of Mechanical Science and Engineering, University of Illinois at Urbana-Champaign, 1206 W. Green Street, Urbana, IL 61801-2906, USA
a r t i c l e
i n f o
Article history: Received 24 June 2013 Accepted 26 July 2013 Available online 26 August 2013 Keywords: Fiber networks Apparent elastic moduli Computational modeling Scale effect Fiber length Fiber orientation Fiber volume fraction
a b s t r a c t Many natural and synthetic materials have fibrous microstructures. In this paper we study computationally the elastic response of three-dimensional random fiber networks with isotropic or preferential fiber orientations. The fibers, represented as linear elastic and isotropic Timoshenko beams, are generated by a Monte Carlo method and joined through rigid fiber–fiber connections. The fiber networks are subjected to displacement boundary conditions, involving uniform in plane normal or shear strains, and their apparent elastic response is evaluated using a finite element method. We investigate the effects of scale (varying fiber length while keeping the specimen size and fiber diameter fixed), fiber orientation, and void volume fraction on the elastic constitutive response of such fiber networks. We find that the apparent stiffness tensor components C1111 and C1212 decrease and have less scatter with increasing scale (and decreasing fiber aspect ratio), higher orientation angle, and higher void volume fraction. Also, the effect of dangling fibers on the void volume fraction is investigated and the resulting elastic response and three ways of estimating the volume fraction of fibers are discussed. 2013 Elsevier B.V. All rights reserved.
1. Introduction Many biological and synthetic materials exhibit fibrous microstructures. Examples include cytoskeleton of a cell [2,12,13,23], blood clots [4], variety of biological tissues including bone, cartilage, muscle, tendon, and ligament [14,15,35,44], and manmade materials such as paper [3,33], filters, fiberboard, textile felts [29], battery technologies [6,41], and insulation materials. Mechanical properties of such fibrous materials depend on many factors including the shape, volume fraction, orientation, arrangement and mechanical properties of fibers. Spatial randomness and heterogeneity of such fiber networks lead to a scatter in their constitutive coefficients. Another important parameter is the window size (region of observation) used for calculations and/or testing, with considerations of that type going back to Delaunay networks [34]. The window size is typically a statistical volume element (SVE) and its scaling trend towards the representative volume element (RVE) has to be established [32]. Studies of the mechanical behavior of fiber networks may be classified into two main categories: (a) phenomenological models and (b) micromechanical models. The first category of models is based on fitting a mathematical equation to an experimentally obtained stress–strain curve. In this approach there is oftentimes ambiguity on how the model parameters are related to the fiber and network properties. This shortcoming is overcome by the
⇑ Corresponding author. Tel.: +1 217 333 9259; fax: +1 217 244 6534. E-mail address: [email protected] (I. Jasiuk). 0927-0256/$ - see front matter 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.commatsci.2013.07.037
second type of models, which aim to homogenize fiber networks analytically or computationally using micromechanics methods. Several micromechanical models involving matrix-inclusion representations, traditionally used in mechanics of composite materials, have been employed to analyze biological materials with fibrous structures. For example, elastic moduli of a single lamella in bone, which consists of preferentially oriented mineralized collagen fibrils, were predicted using a multi-scale matrix-inclusion model accounting for water, collagen, and mineral interactions [16,21]. Yoon and Cowin [47] modeled a single lamella in bone assuming unidirectional cylindrically shaped collagen fibrils embedded in a hydroxyapatite–water composite matrix. Schwartz et al. [37] studied articular cartilage by modeling it as a composite material containing bilinear fibers embedded in a matrix. Another micromechanics-based approach involves fibrous network models. Pioneering analytical work on modeling of fiber networks is due to Cox [7] who predicted elastic properties and strength of isotropic fiber networks in both two-dimensional (2D) and three-dimensional (3D) settings. Cox’s approach was generalized by Narter et al. [29] to anisotropic fiberwebs. The 2D mechanical models of fiberwebs have been reviewed by Shang [38]. More recent 2D analytical studies include [10,46]. A review of analytical results on spring networks is given by Ostoja-Starzewski [31,32] while the computational framework for modeling of biological tissues is summarized by Zohdi [48]. Computational mechanics models of random fiber networks, accounting for fibers discretely, have been employed by a number of researchers [3,5,6,9,17,18,19,23,25,30,33,36,39,40,42,43,46]. They include a finite element model of elastic response of fiber
716
Y. Lee, I. Jasiuk / Computational Materials Science 79 (2013) 715–723
networks representing a single lamella in bone [25] and studies of deformation and failure of stochastic fibrous networks [6,36,42,43]. Other studies treated network-like biological tissues as semi-flexible polymers [17,18,19,20,30]. It was found that fibers with lower bending stiffness are more non-affine, i.e. under a uniform applied strain the local deformation is non-uniform, and the degree of non-affinity decreases with increasing fiber density and scale of observation; see also [10]. Delince and Delannay [9] compare results from modeling to those obtained by testing networks of metallic fibers. The objective of this paper is to investigate the effects of scale (varying fiber length while keeping the specimen size and fiber diameter fixed), fiber orientation, and void volume fraction on the apparent elastic moduli of quasi 3D fiber networks. The fibers are arranged randomly in plane with either preferential or isotropic orientations. Results are analyzed statistically. The effect of dangling fibers on the apparent elastic moduli is also discussed. This paper complements and extends earlier studies in this area.
2. Methodology 2.1. Fiber network generation There are three basic models for representing fiber networks [31]: 1. Two-force members (truss members) carrying axial forces only, connected by frictionless pins. 2. Beam rods carrying axial and shear forces and bending and torsional moments. The rods can be of Bernoulli–Euler or Timoshenko type and they can be connected either by rigid or flexible connections. 3. Members treated as 3D bodies. In this paper we follow the beam network model 2 outlined above (see also [25,33]) and use the following assumptions: 1. Fibers are modeled as Timoshenko beams and linked by rigid connections. 2. The 3D test box (domain used in calculations) has dimensions 10 10 1 units. 3. Fibers are straight, of circular cross-sections with diameters of 1/3 units, and lengths of 2, 3.33, 5, and 10 units. 4. Fibers are arranged randomly using a Poisson point process by generating positions of fibers’ centers with the in-plane scatter in the fiber orientations. The fibers are either preferentially oriented, having a uniform distribution between 10 and +10, or isotropically arranged with a uniform distribution between 90 and +90. 5. Fibers are generated to fill the space to the desired fiber volume fractions up to 80% (corresponding to 20% void volume fraction). 6. Fibers are assumed to be linear elastic and isotropic having Young’s modulus of 1 GPa and Poisson’s ratio of m = 0.3 and they are capable of carrying both tension and compression.
(measured from the center axis) is less than their diameter, otherwise they are not connected. More details on fiber generation are given in Appendix A. In this paper the effects of scale and fiber aspect ratio (varying fiber length while keeping the test box and fiber diameter fixed), fiber orientation and void volume fraction on the constitutive response of fiber networks are studied. The parameters are as follows: lx, ly, and lz are the lengths of the test box along the x, y, and z axes, respectively (Fig. 1), l is the fiber length, r is the radius of the fiber, hmax is the maximum orientation angle between the xaxis and fiber’s axis where the fibers are oriented between the positive hmax and negative hmax, and vfvoid is the void volume fraction. Given the fiber length l, a dimensionless window size (or size of SVE) can be defined as the ratio of the box length over the fiber length, and denoted as d
d¼
length of window lx ¼ length of fiber l
The parameter d, takes on values 1, 2, 3 and 5, for the studied fiber lengths (since the length of window is 10 units and fiber lengths are 2, 3.33, 5, and 10 units). The parameter used for the fiber orientation study is the orientation angle hmax, and the parameter used for the study of the effect of the void volume fraction is the void volume fraction vfvoid. Following [22], the apparent elastic properties are dependent on the size of a window (or specimen size) and boundary conditions, and this gives rise to a hierarchy of bounds
C R ðS R Þ
1
hS t1 i1 6 hS td0 i1 6 hS td i1 6 C eff 6 hC dd i 6 hC dd0 i
6 hC d1 i C V
ð2Þ
where "d0 < d and the inequality between any two tensors is understood as
C P D () ðC DÞ : a : a P 0;
for any tensor aij – 0
Here C is the fourth order stiffness tensor Cijkl and S is the compliance tensor Sijkl where S1 = C. The superscripts R and V denote Voigt and Reuss bounds, respectively, while the superscripts t and d imply traction and displacement boundary conditions, respectively. Following Eq. (2), the effective properties are bounded from above and below by the apparent elastic moduli obtained using displacement and traction boundary conditions, respectively. The larger is the window size d, the closer are the bounds. When d reaches the size of the RVE, the bounds merge and the effective properties are obtained. Eq. (2) directly applies to fiber networks in which the fiber length and diameter are fixed while the window size changes. In this paper, we fixed the window size and the fiber diameter while varying the fiber length. Those two problems are not the same for a given d since the fiber aspect ratio (fiber length:diameter ratio) is different, leading to a different microstructure. Thus, this paper
y-axis 5 2
Random fiber networks are created in a 3D test box by a Monte Carlo method. For each realization of the network, a different seed number is chosen and the locations of the centers of the fibers and their individual in-plane orientations are generated from a Poisson point field [26]; twenty realizations are generated for each data point. Thus, each generated fiber lies in plane (i.e., is parallel to the two largest surfaces of the test box). Fibers are linked by rigid connections when the nearest distance between two fibers
ð1Þ
7 4 6
3 Test box
z-axis 1 x-axis
Fig. 1. Test box and regions of trimming when the other end of the fiber is in other boxes.
717
Y. Lee, I. Jasiuk / Computational Materials Science 79 (2013) 715–723
E 2ð1 þ mÞ
investigates a combined effect of scale and microstructure, namely as d increases the fiber aspect ratio decreases.
G¼
2.2. Preprocessing in MATLAB
where m is the Poisson ratio, E is the elastic modulus, G is the shear modulus, and r is the radius of the fiber cross-section.
2.2.1. Input parameters The coordinates of the nodes and the elements which constitute the network are imported into a MATLAB program [28]. The information on the fiber network generation and the input parameters is given in Section 2.1. The input parameters are also summarized in Table 1. 2.2.2. Trimming of fibers In the fiber generation process, some fibers might protrude or fall outside of the test box. Therefore it is necessary to trim the fibers so that all of them are contained in the box. Fig. 1 shows the regions where the fibers might fall or protrude into. The procedure that ensures that the fibers do not protrude into negative regions of the box is summarized in the Appendix A. 2.2.3. Removing elements and clusters not contributing to the constitutive response When a loading is applied to the network at external boundaries, some of the randomly generated fibers in the test box may not experience the load, resulting in an overestimated volume fraction of the fibers. This will lead to underestimation of the constitutive response of the fiber network. Thus, for our calculations we remove the elements which do not contribute to the overall response of the network, except when we investigate this effect. The excluded elements include an element or a cluster, which does not experience any load (or is known as dangling) and those that experience loading only at one point. This concept has been verified through Abaqus; [1] whereby, for each condition stated above, the strain energy of the removed elements is negligible compared with the overall strain energy of the fiber network system.
2.3.1. Choosing a beam element There are two types of beams that can be chosen in Abaqus, namely the Euler–Bernoulli beam and the Timoshenko beam. The difference between these two beam elements is that there is a second-order spatial derivative present in the Timoshenko beam, resulting in the lower stiffness of the beam and larger deflection under loading. It is also well known that the Timoshenko beam is a well posed model and thus able to handle stubby beams such as short fiber segments in the network. Thus, we select it over the Euler–Bernoulli beam for modeling fibers. The shear coefficient for a circular cross-section is defined as [24]:
6ð1 þ mÞ2 7 þ 12m þ 4m2
SR2 ¼ kGpr 2 ;
2.3.3. Specifying boundary conditions We apply the in-plane uniaxial and shear displacements (with uniform strains) on the six faces of the test box in order to obtain the apparent stiffness coefficients C1111 and C1212, respectively, of the studied fiber networks. We denote ui as a displacement, xi as a global coordinate, eij as a uniform strain, with subscripts i and j taking values 1, 2, 3, which correspond to x, y, and z directions, respectively. To obtain C1111, the uniaxial displacement (with a uniform normal strain e11) is applied as follows on all six boundaries:
u1 ¼ e11 x1 u2 ¼ 0
For determining C1212, the applied in-plane shear displacement (with a uniform shear strain e12) is taken as:
u1 ¼ e12 x2 u2 ¼ e12 x1
SR3 ¼ 0:25
r11 1 0 C1111 Br C BC B 22 C B 2211 C B B B r33 C B C3311 C B B Br C ¼ BC B 23 C B 2311 C B B @ r31 A @ C3111 r12 C1211 0
C1122
C1133
C1123
C1131
C2222
C2233
C2223
C2231
C3322
C3333
C3323
C3331
C2322 C3122
C2333 C3133
C2323 C3123
C2331 C3131
C1222
C1233
C1223
C1231
e11 1 C B C2212 C CB e22 C CB C C B C3312 e C CB 33 C C B C2312 C CB 2e23 C CB C C3112 A@ 2e31 A C1212 2e12 C1112
10
ð8Þ Since only either normal or shear strains (e11or e12) are applied (in the x–y plane), the constitutive Eq. (8) reduce to the following equations:
r11 ¼ C1111 e11
ð9Þ
r12 ¼ 2C1212 e12
ð10Þ
respectively, where r11 and r12 are normal and shear stresses, respectively. The elastic strain energy W of the system is given by:
ð4Þ
V rij eij 2
ð11Þ
By incorporating Eq. (9) for the uniaxial case, and assuming that the system experiences a global uniform strain, C1111 can be found using Eqs. (12)–(14).
Table 1 Parameters used for generating fiber networks.
Values
ð7Þ
u3 ¼ 0
W¼
Parameters
ð6Þ
u3 ¼ 0
ð3Þ
The shear rigidity (SR) is then imported into Abaqus using the following definitions:
SR1 ¼ kGpr 2 ;
2.3.2. Specifying beam orientations Three dimensional (3D) beam elements are used in the simulations. The orientation of a beam is defined in terms of a local, right handed system. An additional node of the beam axis is defined and included in the element’s connectivity list. Abaqus is then able to compute three vectors through these nodes. The dot and cross product of these vectors give an orientation of the beam.
2.3.4. Calculating C1111 and C1212 stiffness constants The constitutive equations relate stresses to strains as follows:
2.3. Importing into Abaqus
k¼
ð5Þ
Box dimension (lx ly lz)
Radius of fiber r
vfvoid (%)
10 10 1
1/6
10, 20, 40, 60, 80
hmax ()
d
10, 90
1, 2, 3, 5
W¼
V C 1111 e211 2
V ¼ lx ly lz
ð12Þ ð13Þ
718
Y. Lee, I. Jasiuk / Computational Materials Science 79 (2013) 715–723
W¼
N X ei
ð14Þ
i¼1
where V is the volume of the test box, ei is the strain energy of a particular element, and N is the total number of elements in the system. For the shear loading case, C1212 can be found from Eqs. (13)–(15).
W ¼ 2VC 1212 e212
ð15Þ
All the information necessary to calculate C1111 and C1212 is obtained from the output file of Abaqus postprocessing module. For simplicity, the value C1111 is normalized by Cfiber 1111 and C1212 is norfiber fiber malized by Cfiber 1212 , where C1111 and C1212 are the linear elastic and isotropic stiffness components of the fiber. They can be calculated through the stress–strain relationship for an isotropic material. Cfiber 1111 is found from e.g. [8]:
r11
E ¼ ½ð1 mÞe11 þ mðe22 þ e33 Þ ð1 þ mÞð1 2mÞ
ð16Þ
where E, m are the Young’s modulus and Poisson’s ratio. After simplification, using Eq. (16), Cfiber 1111 is given as follows:
Cfiber 1111 ¼
Eð1 mÞ ð1 þ mÞð1 2mÞ
ð17Þ
Cfiber 1212 is the same as the shear modulus of an isotropic fiber, given by Eq. (5). 3. Results and discussion Fiber networks with the preferential and isotropic orientations (one realization for each case) under a uniaxial displacement are shown in Fig. 2(a) and (b), respectively. Note that in our model the fibers can resist compression. Displacements of the fiber endpoints are shown by arrows (in orange in color). Next, the results on the effects of scale and fiber aspect ratio (plotted in terms of window size d), the fiber orientation, and the void (or fiber) volume fraction on the stiffness components C1111 and C1212 are presented. Twenty realizations were considered for each data point. Fig. 3 shows the plot of the normalized apparent stiffness coefficient C1111 =C fiber 1111 versus window size d for fiber networks with preferentially oriented fibers (10 and +10) with the void volume fraction of 20% (i.e. 80% fiber volume fraction), obtained by applying a uniaxial displacement with a uniform strain e11. The normalized apparent stiffness coefficient C1111 =Cfiber 1111
Fig. 3. Non-dimensionalized stiffness component C1111 =Cfiber 1111 versus window size d for the fiber network with a preferential fiber orientation 10 and +10 for 20% void volume fraction (or 80% fiber volume fraction) with standard deviation bar. Twenty realizations were calculated for each window size. In calculations test box size and fiber diameter was fixed while fiber length was varied. Thus, as the window size increases the fiber aspect ratio decreases.
decreases with the increasing d (and decreasing fiber aspect ratio). The fact that the stiffness decreases with increasing window size is consistent with [22]. Following [22], when d becomes large, leading to the RVE size, the normalized apparent stiffness coefficient C1111 =Cfiber 1111 converges from above to an effective stiffness. Recall that the studied problem involves a change in the fiber aspect ratio with the change in d. Thus, our simulations include the effect of scale d combined with changes in the fiber aspect ratio. Also, the scatter in the normalized C1111 =Cfiber 1111 decreases as d increases as expected [32]. The application of traction boundary conditions with a uniform stress would give a lower bound for the apparent elastic stiffness. However, several issues arise when applying the traction boundary conditions to fiber networks. Since fibers are randomly oriented the cross-sectional areas of beam elements are different for each fiber at the surface (boundary) of the box. Also, it is impossible to unambiguously apply a constant stress on the boundary since fibers cross the box boundary at different angles and locations [32,34]. Therefore, we do not apply this second type of boundary conditions. Also, for simplicity of notation we do not include superscripts d. Figs. 4 and 5 show the probability plots to fit the data for C1111 of fiber networks with preferentially oriented fibers (10 and
Fig. 2. Beam network loaded under a uniaxial displacement boundary condition (a) with the preferential fiber orientation (10 and +10), (b) with the isotropic fiber orientation (90 and +90). Displacements of endpoints are shown with arrows (in color in orange). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Y. Lee, I. Jasiuk / Computational Materials Science 79 (2013) 715–723
+10) with 20% void volume fraction (80% fiber volume fraction) for window sizes 1 and 5. Similar curves were obtained for the window sizes 2 and 3. Ten distribution functions: Beta, Gamma, Logistic, Gumbel Min, Rayleigh, Weibull, v and other (see Appendix B) were considered and the Kolmogorov–Smirnov test was used to determine the best distribution function fits for each of the four window sizes. We find that the probability distribution that best fits the data is the Logistic, followed by Gamma and Beta. Note that the commonly used Weibull distribution does not give a good fit to the data. Next, we consider fiber networks with isotropically arranged fibers (90 and +90) and investigate how d affects the apparent constitutive response under a uniform uniaxial strain (Fig. 6). The trends are similar to those for networks with preferential fiber orientations (Fig. 3). Again, the normalized apparent stiffness component C1111 =Cfiber 1111 decreases with increasing d but at lower rate than for the preferentially oriented fiber case. The constitutive response for the isotropically distributed fiber network is lower than that of the preferentially oriented fiber network, in particular for small d. Again, this is expected since in such networks the fibers, which are arranged preferentially in the loading direction and have large aspect ratios, better transfer loads than those in the isotropic networks. We also observe that the number of elements removed in order to achieve a given value of the void (or fiber) volume fraction stay about the same for an isotropic network for different window
719
sizes while the number of the elements removed in a preferentially oriented network increases as the window size increases. Next, the effects of the void (or fiber) volume fraction on the stiffness constant C1111 =Cfiber 1111 are investigated, as shown in Figs. 7 and 8 for the networks with isotropic fiber orientation for d = 2 and 3 and in Fig. 9 for the networks with preferential fiber orientation for d = 2. As expected, the stiffness component decreases with increasing void volume fraction. The initial decrease is linear indicating affine deformations. For the isotropic case the decrease is linear up to vfvoid of about 80% when d = 2 and up to about vfvoid of 60% for d = 3. Note C1111 =Cfiber 1111 does not reach zero because our formulation ensures that all the structural elements are connected to clusters which can sustain loads. Moreover, the scatter in constitutive values decreases with increasing void volume fraction. Also, the stiffness component C1111 =Cfiber 1111 of the fiber network with preferentially oriented fibers is higher than that for an isotropic fiber network as expected since the fibers aligned in the loading direction can carry the load more efficiently. The same study as for C1111 was performed to investigate the effect of scale (with changing fiber aspect ratio), void volume fraction and fiber orientation on the normalized apparent shear stiffness component C1212 =Cfiber 1212 as illustrated in Figs. 10–14. We observe that the beam network with isotropic fiber orientation can withstand a shear loading better than the preferentially oriented beam network, as shown by comparing Figs. 12 and 14 for the window
Fig. 4. Fits for C1111 by several distribution functions for the case of d = 1 (for preferential fiber orientation).
Fig. 5. Fits for C1111 by several distribution functions for the case of d = 5 (for preferential fiber orientation).
720
Y. Lee, I. Jasiuk / Computational Materials Science 79 (2013) 715–723
Fig. 6. C1111 =Cfiber 1111 versus window size d for a fiber network with an isotropic fiber orientation 90 and +90 for 20% void volume fraction (with standard deviation bar).
Fig. 7. C1111 =Cfiber 1111 versus void volume fraction of beam network with an isotropic fiber orientation and d = 2 (with standard deviation bar).
Fig. 9. C1111 =Cfiber 1111 versus void volume fraction of beam network with a preferential fiber orientation and d = 2 (with standard deviation bar).
Fig. 10. Non-dimensionalized stiffness component C1212 =Cfiber 1212 versus window size for a fiber network with a preferential fiber orientation for 20% void volume fraction (with standard deviation bar).
Fig. 8. C1111 =Cfiber 1111 versus void volume fraction of beam network with an isotropic fiber orientation and d = 3 (with standard deviation bar). Fig. 11. C1212 =Cfiber 1212 versus window size for a fiber network with an isotropic fiber orientation for 20% void volume fraction (with standard deviation bar).
size 2. In addition, the preferentially oriented fiber network decays faster than the isotropic fiber network as seen by comparing Figs. 10 and 11. Despite the above mentioned differences, the volume fraction studies (Figs. 12 and 13) of C1212 share similarity with C1111. The stiffness C1212 =Cfiber 1212 decreases linearly with increasing void volume fraction indicating affine deformations. For the isotropic case the decrease is linear up to vfvoid of about 80% when d = 2 and up to about vfvoid of 60% for d = 3 as for the C1111 case.
Next, we investigate the effect of dangling fibers, i.e. fibers that do not carry loads in the network. In order to verify the removal of elements which do not sustain loads, the strain energy densities of the network with dangling fibers and the network with removed elements are compared. A check on the fiber network confirmed that there were no fibers that experience zero strain energy. There are various ways of defining the void (or fiber) volume fraction in a network. When a network is generated, whether it
721
Y. Lee, I. Jasiuk / Computational Materials Science 79 (2013) 715–723
Fig. 12. C1212 =Cfiber 1212 versus void volume fraction of beam network with an isotropic fiber orientation and d = 2 (with standard deviation bar).
Fig. 13. C1212 =Cfiber 1212 N versus void volume fraction of beam network with an isotropic fiber orientation and d = 3 (with standard deviation bar).
the preferentially oriented fiber network in Table 4. For the isotropic fiber network, the normalized C1111 value, fiber volume fraction and percentage of elements removed do not vary significantly (Table 3). This is to be expected since the elements in the isotropic fiber network are more easily connected to each other when compared with the elements in a preferentially oriented fiber network. Therefore, there are fewer dangling fibers in the isotropic fiber network leading to relatively constant values of parameters under investigation. For the fiber network with the preferentially oriented fibers, Cases II and III have the similar normalized C1111 value, with discrepancy of about 2%, compared with a large increase for Case I (Table 4). Note that Case II is Case III plus dangling fibers, and since dangling fibers do not carry any load, we should expect both cases to have similar normalized C1111 value. As for Case III, it has a higher normalized C1111 value since all the fibers are interconnected and contribute to the load sharing capacity of the network. Besides, Case III has a lower fiber volume fraction compared with the other two cases since approximately 16% of dangling fibers have been removed. The above analysis gives three scenarios for different definitions of fiber volume fraction of a network. Occasionally, Case II is being used in modeling of fiber networks. However, the volume fraction might be overestimated since the actual volume fraction might be lower due to dangling fibers. In addition, the percentage of dangling fibers in the network is a variable that has to be taken into account if the isotropic network is to be generated. Therefore, Case I gives a more realistic representation of a network. Another issue that needs to be considered is the effect of an intersected volume on the calculated volume fraction of fibers. The fiber volume fraction vf can be estimated as follows:
Table 2 Different definitions of the volume fraction of a network. Note that Case I was used in data reported in Figs. 3–14. Case I
Case II
Case III
1. Generate the fiber network until 80% fiber volume fraction 2. Remove any dangling fibers 3. Add more fibers and repeat step 2 until fiber volume fraction reaches 80%
1. Generate the network until 80% fiber volume fraction 2. Do not remove any dangling fiber
1. Generate the network until 80% fiber volume fraction 2. Remove any dangling fibers 3. Do not add more fibers
Table 3 Effect of different methods of generating fiber network with an isotropic orientation with a 20% void fraction (80% fiber fraction) for d = 2 on C1111 =Cfiber 1111 . Fig. 14. C1212 =Cfiber 1212 versus void volume fraction of beam network with a preferential fiber orientation and d = 2 (with standard deviation bar).
is through a manufacturing process as in an air filter case, or through a computer simulation, one may find dangling fibers not contributing to load sharing but occupying a definite amount of space in the network. This would overestimate the volume fraction in the network, leading to incorrect representation of the load sharing capability of the network. In this paper, we consider three Cases (I, II and III) that can be used to define the fiber volume fraction of a network, summarized in Table 2. In this analysis the 20% void (80% fiber) volume fraction and d = 2 are used. We illustrate the implication of these three different Cases (I–III) on C1111 for the isotropic fiber network in Table 3 and for
Cases
C1111 =Cfiber 1111
Fiber volume fraction (%)
Elements removed (%)
Case I Case II Case III
0.0968 0.0966 0.0965
80 80 80
0.545 0 0.547
Table 4 Effect of different methods of generating fiber network with a preferential fiber orientation with a 20% void fraction (80% fiber fraction) for d = 2 on C1111 =Cfiber 1111 . Cases
C1111 =Cfiber 1111
Fiber volume fraction (%)
Elements removed (%)
Case I Case II Case III
0.125 0.0792 0.0778
80 80 62.6
8.96 0 16.83
722
vf
¼
Y. Lee, I. Jasiuk / Computational Materials Science 79 (2013) 715–723 n X
pr2 li
ð18Þ
i¼1
where r is the radius and l is the length of each fiber, respectively. This calculation does not take into account the intersected volume at each joint. Thus, the total volume fraction of fibers is exaggerated and therefore the constitutive response is underestimated. This effect is more pronounced when the fibers are aligned, almost parallel to each other, as in the preferential fiber orientation case. Furthermore, when the window sizes get smaller, or the fibers get shorter compared with the length of a test box, the effect of the intersected volume is higher as well. Therefore, there is a need to calculate these intersected volumes in order to obtain more accurate results. More details on the determination of the intersected volume are given in [27]. In summary, we studied numerically the effect of window size (and fiber aspect ratio), fiber orientation and void volume fraction on the apparent elastic stiffness components C1111 and C1212 of fiber networks with isotropic and preferential fiber orientations. In our parametric study we considered only several sets of inputs to investigate the trends; more comprehensive parametric studies can be done in the future. Also, we combined the effect of scale and fiber aspect ratio while these effects could also be studied separately. Also, a nonlinear fiber response can be incorporated in the Abaqus in order to capture nonlinear effects in the constitutive response of fiber networks. Strength of the network can also be investigated through element removal and fiber–fiber debonding once fibers reach a critical stress level. Last but not least, the entire stiffness tensor can be found through applications of additional loadings. The current study can be extended to biological materials, such as an articular cartilage which is known to have three distinct morphological zones with different fiber orientations in each zone [45]. The knowledge how fiber properties and structural characteristics of fiber networks affect their constitutive response is important in studies of healthy and diseased biological tissues and variety of synthetic fibrous materials. 4. Summary We modeled random fiber networks as beam networks and predicted two apparent stiffness components C1111 and C1212 as a function of several parameters that affect their means and scatter. This problem was solved computationally using a finite element method. Two fiber orientations were considered, isotropic and preferential, for various void volume fractions and scales of observation (and fiber aspect ratio). The effect of dangling fibers was also investigated and was shown to play a role in calculations, in particular, for the preferentially oriented network. The calculation of fiber volume fraction was also briefly discussed. Acknowledgement We gratefully acknowledge support of the National Science Foundation (CMMI 0927909 – Dr. Ken Chong). Appendix A In the fiber generation process, some fibers might protrude or fall outside of the test box. Therefore it is necessary to trim the fibers so that all of them are contained in the box. Fig. 1 shows the regions where the fibers might fall or protrude into. The following formulations ensure that the fibers do not protrude into negative regions of the box. For the region 1, the following equations are used:
Y 2new ¼
Z 2new ¼
lx X 1old Y 2old Y 1old
þ Y 1old
ðA:1Þ
þ Z 1old
ðA:2Þ
X 2old X 1old lx X 1old Z 2old Z 1old X 2old X 1old
X 2new ¼ lx
ðA:3Þ
where X, Y, Z are the local coordinates of the fiber endpoints (node), the superscripts 1 and 2 denote the node inside the test box and the node outside the test box, respectively, while the subscripts old and new denote the coordinates before and after trimming, respectively. For the region 2, the following equations are used:
X 2new ¼
Z 2new ¼
ly Y 1old X 2old X 1old Y 2old Y 1old ly Y 1old Z 2old Z 1old Y 2old Y 1old
þ X 1old
ðA:4Þ
þ Z 1old
ðA:5Þ
Y 2new ¼ ly
ðA:6Þ
For the region 3, the following equations are used:
X 2new ¼
Y 2new ¼
lz Z 1old X 2old X 1old Z 2old Z 1old lz Z 1old Y 2old Y 1old
Z 2new ¼ lz
Z 2old Z 1old
þ X 1old
ðA:7Þ
þ Y 1old
ðA:8Þ
ðA:9Þ
For region 4, where we do not know a priori whether a fiber will protrude through region 1 or region 2, we determine it in the following way. We first assume that the fiber will protrude into region 4 through region 1. Eq. (A.1) is then used to find the coordinate Y of the new node. The coordinate Y of this node is then compared with ly. If ly is larger than the new coordinate Y, then the fiber is confirmed to protrude from the test box into region 4 through region 1. Therefore, Eqs. (A.1)–(A.3) are used to find the node of intersection on the face of the test box. On the contrary, if ly is smaller than the new coordinate Y, then Eqs. (A.4)–(A.6) are used to find the new coordinates since the fiber will protrude into region 4 through region 2. The same concept can be applied to fibers protruding into region 5. Eqs. (A.4)–(A.6) are used if the fibers will protrude through region 2, and Eqs. (A.7)–(A.9) are used if the fibers will protrude through region 3. Similarly for region 6, Eqs. (A.1)–(A.3) are used for the fibers protruding through region 1, and Eqs. (A.7)–(A.9) are used for the fibers protruding through region 3. The region 7 requires a more careful treatment since a fiber will generally go through two regions before reaching region 7. First, let us assume that the fiber will protrude from the test box into region 7 through region 1, and, therefore, Eqs. (A.2) and (A.3) are used. If both ly and lz are larger than the new coordinate Y and the new coordinate Z, respectively, then Eqs. (A.1)–(A.3) are used since the fiber will go through region 1 and either region 4 or region 6 before reaching region 7. However, if either ly or lz is smaller than the new coordinate Y or new coordinate Z, respectively, then the fiber does not protrude through region 1, but through region 2 or region 3 instead. In order to determine the protruded region, let us assume that the fiber will protrude through region 2 and therefore Eqs. (A.4) and (A.5) are used. If lz is larger than the new
Y. Lee, I. Jasiuk / Computational Materials Science 79 (2013) 715–723
coordinate Z, then the fiber protrudes through region 2 and then through either region 4 or region 5 before reaching region 7. Therefore, Eqs. (A.4)–(A.6) are used. On the contrary, Eqs. (A.7)–(A.9) are used for the fiber protruding through region 3. Note that some elements might have two nodes at the same place if the elements protrude out of the test box through its faces. In order to avoid this situation, the length of each element is checked so that it is greater than zero, meaning that one node does not overlap with another. This trimming action can be used to study the constitutive response of a certain part of the fiber network as well. The test box can be placed in a specific area of the network, while the rest of the area can be trimmed off.
Appendix B For reference, probability distribution functions used in this study are listed. Beta
f ðx; a1 ; a2 ; d1 ; d2 Þ ¼
a1 1 a2 1 1 1 dxd 2 d1
xd1 d2 d1
ðd2 d1 ÞBða1 ; a2 Þ
ðB:1Þ
Chi
f ðx; a; b; eÞ ¼
a1 a 1 xe 2 212 x e a e2ð b Þ b bC 2
ðB:2Þ
k½kðx eÞj1 e½kðxeÞ CðjÞ
ðB:3Þ
Gamma
f ðx; j; k; eÞ ¼
Gumbel Max
f ðx; l; aÞ ¼ aeðaðlxÞe
aðlxÞ Þ
ðB:4Þ
Gumbel Min
f ðx; l; aÞ ¼ aeðaðxlÞe
aðxlÞ Þ
ðB:5Þ
Logistic
f ðx; a1 ; a2 Þ ¼
1 4a2 cosh
2 xa1 2a2
ðB:6Þ
Normal
1 1 xl 2 f ðx; l; rÞ ¼ pffiffiffiffiffiffiffi e2ð r Þ 2pr
ðB:7Þ
Rayleigh
f ðx; a; eÞ ¼
xe
a2
1 xe 2
e2ð a Þ
ðB:8Þ
Rectangular
f ðx; a1 ; eÞ ¼
1 a1
ðB:9Þ
Weibull
f ðx; l; j; eÞ ¼
xe le
j1
le
xe
j
jeðleÞ
ðB:10Þ
723
References [1] Abaqus, v. 6.7, ed., Providence, Rhode Island, Dassault Systemes Simulia Corp. [2] B. Alberts, A. Johnson, J. Lewis, M. Raff, K. Roberts, P. Walter, Molecular Biology of the Cell, second ed., Garland Science, New York, NY, 2002. [3] M. Alkhagen, S. Toll, Journal of Applied Mechanics – Transactions of the ASME 74 (2007) 723–731. [4] R.D. Averett, B. Menn, E.H. Lee, C.C. Helms, T. Barker, M. Guthold, Biophysical Journal 103 (2012) 1537–1544. [5] P.L. Chandran, V.H. Barocas, Journal of Biomechanical Engineering – Transactions of the ASME 128 (2006) 259–270. [6] X. Cheng, C. Wang, A.M. Sastry, S.B. Choi, Journal of Engineering Materials and Technology 121 (1999) 515–523. [7] H.L. Cox, British Journal of Applied Physics 3 (1952) 72–79. [8] R.R.J. Craig, in: Mechanics of Materials, second ed., Wiley, New York, 2000. [9] M. Delince, F. Delannay, Acta Materialia 52 (2004) 1013–1022. [10] B.A. DiDonna, T.C. Lubensky, Physical Review E 72 (2005) 066619. [12] M.L. Gardel, F. Nakamura, J.H. Hartwig, J.C. Crocker, T.P. Stossel, D.A. Weitz, Proceedings of the National Academy of Sciences of the United States of America 103 (2006) 1762–1767. [13] M.L. Gardel, J.H. Shin, F.C. MacKintosh, L. Mahadevan, P. Matsudaira, D.A. Weitz, Science 304 (2004) 1301–1305. [14] M.M. Giraud-Guille, Calcified Tissue International 42 (1988) 167–180. [15] M.M. Giraud-Guille, L. Besseau, R. Martin, Journal of Biomechanics 36 (2003) 1571–1579. [16] E. Hamed, Y. Lee, I. Jasiuk, Acta Mechanica 213 (2010) 131–154. [17] H. Hatami-Marbini, R.C. Picu, Analytical approach to quantifying the nonaffine behavior of fiber networks, in: Materials Research Society Symposium Proceedings, vol. 1060, 2008a, pp. 33–37. [18] H. Hatami-Marbini, R.C. Picu, Physical Review E 77 (2008) 062103. [19] H. Hatami-Marbini, R.C. Picu, Acta Mechanica 205 (2009) 77–84. [20] D.A. Head, A.J. Levine, F.C. MacKintosh, Physical Review E 68 (2003) 061907. [21] C. Hellmich, J. Barthélémy, L. Dormieux, European Journal of Mechanics, A/ Solids 23 (2004) 783–810. [22] C. Huet, Journal of Mechanics and Physics of Solids 38 (1990) 813–841. [23] E.M. Huisman, T. van Dillen, P.R. Onck, E. van der Giessen, Physical Review Letters 99 (2007) 208103. [24] J.R. Hutchinson, Journal of Applied Mechanics – Transactions of the ASME 68 (2001) 87–92. [25] I. Jasiuk, M. Ostoja-Starzewski, Biomech Model Mechanobiology 3 (2004) 67–74. [26] D. Jeulin, Random Structure Models for Homogenization and Fracture Statistics, in: D. Jeulin, M. Ostoja-Starzewski (Eds.), Mechanics of Random and Multiscale Microstructures CISM Courses and Lectures, Springer, Berlin, Heidelberg, New York, 2001, pp. 33–91. [27] Y.H. Lee, Non-Invasive Characterization of Bone Tissue and Modeling of Bone at Micron and Sub-Micron Levels, in: Mechanical Science and Engineering, University of Illinois, Urbana, 2010, p. 175. [28] MATLAB, TheMathWorks, Natick, Massachusetts, 2008. [29] M.A. Narter, S.K. Batra, D.R. Buchanan, Proceedings of the Royal Society London, Series A – Mathematical Physical and Engineering Sciences 455 (1999) 3543–3563. [30] P.R. Onck, T. Koeman, T. van Dillen, E. van der Giessen, Physical Review Letters 95 (2005) 178102. [31] M. Ostoja-Starzewski, Applied Mechanics Reviews 55 (2002) 35–60. [32] M. Ostoja-Starzewski, Microstructural Randomness and Scaling in Mechanics of Materials, Chapman & Hall/CRC/Taylor & Francis, Boca Raton, FL, 2008. [33] M. Ostoja-Starzewski, D.C. Stahl, Journal of Elasticity 60 (2000) 131–149. [34] M. Ostoja-Starzewski, C. Wang, Acta Mechanica 80 (1989) 61–80. [35] M. Petrty´l, J. Hert, P. Fiala, Journal of Biomechanics 29 (1996) 161–169. [36] A.M. Sastry, C.W. Wang, L. Berhan, Key Engineering Materials 200 (2001) 229– 250. [37] M.H. Schwartz, P.H. Leo, J.L. Lewis, Journal of Biomechanics 27 (1994) 865–873. [38] P.P. Shang, Characterization of Anisotropy of Spunbonded Fiberwebs (nowovens), North Carolina State University, 1985. [39] D.C. Stahl, S.M. Cramer, ASME Journal of Engineering Materials and Technology 120 (1998) 126–130. [40] T. Stylianopoulos, V.H. Barocas, Computer Methods in Applied Mechanics and Engineering 196 (2007) 2981–2990. [41] C. Wang, X. Cheng, A.M. Sastry, S.B. Choi, Journal of Engineering Materials and Technology 121 (1999). 530–513. [42] C.W. Wang, L. Berhan, A.M. Sastry, Journal of Engineering Materials and Technology 122 (2000) 450–459. [43] C.W. Wang, A.M. Sastry, Journal of Engineering Materials and Technology 122 (2000) 460–468. [44] J.W. Weisel, Biophysical Chemistry 112 (2004) 267–276. [45] J.Z. Wu, W. Herzog, Journal of Biomechanics 35 (2002) 931–942. [46] X.-F. Wu, Y.A. Dzenis, Journal of Applied Physics 98 (2005) 093501. [47] Y.J. Yoon, S.C. Cowin, Biomechanics and Modeling in Mechanobiology 7 (2008) 1–11. [48] T.I. Zohdi, Computational Methods in Applied Mechanics and Engineering 196 (2007) 2972–2980.
Contents lists available at ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
Apparent elastic properties of random fiber networks Yikhan Lee, Iwona Jasiuk ⇑ Department of Mechanical Science and Engineering, University of Illinois at Urbana-Champaign, 1206 W. Green Street, Urbana, IL 61801-2906, USA
a r t i c l e
i n f o
Article history: Received 24 June 2013 Accepted 26 July 2013 Available online 26 August 2013 Keywords: Fiber networks Apparent elastic moduli Computational modeling Scale effect Fiber length Fiber orientation Fiber volume fraction
a b s t r a c t Many natural and synthetic materials have fibrous microstructures. In this paper we study computationally the elastic response of three-dimensional random fiber networks with isotropic or preferential fiber orientations. The fibers, represented as linear elastic and isotropic Timoshenko beams, are generated by a Monte Carlo method and joined through rigid fiber–fiber connections. The fiber networks are subjected to displacement boundary conditions, involving uniform in plane normal or shear strains, and their apparent elastic response is evaluated using a finite element method. We investigate the effects of scale (varying fiber length while keeping the specimen size and fiber diameter fixed), fiber orientation, and void volume fraction on the elastic constitutive response of such fiber networks. We find that the apparent stiffness tensor components C1111 and C1212 decrease and have less scatter with increasing scale (and decreasing fiber aspect ratio), higher orientation angle, and higher void volume fraction. Also, the effect of dangling fibers on the void volume fraction is investigated and the resulting elastic response and three ways of estimating the volume fraction of fibers are discussed. 2013 Elsevier B.V. All rights reserved.
1. Introduction Many biological and synthetic materials exhibit fibrous microstructures. Examples include cytoskeleton of a cell [2,12,13,23], blood clots [4], variety of biological tissues including bone, cartilage, muscle, tendon, and ligament [14,15,35,44], and manmade materials such as paper [3,33], filters, fiberboard, textile felts [29], battery technologies [6,41], and insulation materials. Mechanical properties of such fibrous materials depend on many factors including the shape, volume fraction, orientation, arrangement and mechanical properties of fibers. Spatial randomness and heterogeneity of such fiber networks lead to a scatter in their constitutive coefficients. Another important parameter is the window size (region of observation) used for calculations and/or testing, with considerations of that type going back to Delaunay networks [34]. The window size is typically a statistical volume element (SVE) and its scaling trend towards the representative volume element (RVE) has to be established [32]. Studies of the mechanical behavior of fiber networks may be classified into two main categories: (a) phenomenological models and (b) micromechanical models. The first category of models is based on fitting a mathematical equation to an experimentally obtained stress–strain curve. In this approach there is oftentimes ambiguity on how the model parameters are related to the fiber and network properties. This shortcoming is overcome by the
⇑ Corresponding author. Tel.: +1 217 333 9259; fax: +1 217 244 6534. E-mail address: [email protected] (I. Jasiuk). 0927-0256/$ - see front matter 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.commatsci.2013.07.037
second type of models, which aim to homogenize fiber networks analytically or computationally using micromechanics methods. Several micromechanical models involving matrix-inclusion representations, traditionally used in mechanics of composite materials, have been employed to analyze biological materials with fibrous structures. For example, elastic moduli of a single lamella in bone, which consists of preferentially oriented mineralized collagen fibrils, were predicted using a multi-scale matrix-inclusion model accounting for water, collagen, and mineral interactions [16,21]. Yoon and Cowin [47] modeled a single lamella in bone assuming unidirectional cylindrically shaped collagen fibrils embedded in a hydroxyapatite–water composite matrix. Schwartz et al. [37] studied articular cartilage by modeling it as a composite material containing bilinear fibers embedded in a matrix. Another micromechanics-based approach involves fibrous network models. Pioneering analytical work on modeling of fiber networks is due to Cox [7] who predicted elastic properties and strength of isotropic fiber networks in both two-dimensional (2D) and three-dimensional (3D) settings. Cox’s approach was generalized by Narter et al. [29] to anisotropic fiberwebs. The 2D mechanical models of fiberwebs have been reviewed by Shang [38]. More recent 2D analytical studies include [10,46]. A review of analytical results on spring networks is given by Ostoja-Starzewski [31,32] while the computational framework for modeling of biological tissues is summarized by Zohdi [48]. Computational mechanics models of random fiber networks, accounting for fibers discretely, have been employed by a number of researchers [3,5,6,9,17,18,19,23,25,30,33,36,39,40,42,43,46]. They include a finite element model of elastic response of fiber
716
Y. Lee, I. Jasiuk / Computational Materials Science 79 (2013) 715–723
networks representing a single lamella in bone [25] and studies of deformation and failure of stochastic fibrous networks [6,36,42,43]. Other studies treated network-like biological tissues as semi-flexible polymers [17,18,19,20,30]. It was found that fibers with lower bending stiffness are more non-affine, i.e. under a uniform applied strain the local deformation is non-uniform, and the degree of non-affinity decreases with increasing fiber density and scale of observation; see also [10]. Delince and Delannay [9] compare results from modeling to those obtained by testing networks of metallic fibers. The objective of this paper is to investigate the effects of scale (varying fiber length while keeping the specimen size and fiber diameter fixed), fiber orientation, and void volume fraction on the apparent elastic moduli of quasi 3D fiber networks. The fibers are arranged randomly in plane with either preferential or isotropic orientations. Results are analyzed statistically. The effect of dangling fibers on the apparent elastic moduli is also discussed. This paper complements and extends earlier studies in this area.
2. Methodology 2.1. Fiber network generation There are three basic models for representing fiber networks [31]: 1. Two-force members (truss members) carrying axial forces only, connected by frictionless pins. 2. Beam rods carrying axial and shear forces and bending and torsional moments. The rods can be of Bernoulli–Euler or Timoshenko type and they can be connected either by rigid or flexible connections. 3. Members treated as 3D bodies. In this paper we follow the beam network model 2 outlined above (see also [25,33]) and use the following assumptions: 1. Fibers are modeled as Timoshenko beams and linked by rigid connections. 2. The 3D test box (domain used in calculations) has dimensions 10 10 1 units. 3. Fibers are straight, of circular cross-sections with diameters of 1/3 units, and lengths of 2, 3.33, 5, and 10 units. 4. Fibers are arranged randomly using a Poisson point process by generating positions of fibers’ centers with the in-plane scatter in the fiber orientations. The fibers are either preferentially oriented, having a uniform distribution between 10 and +10, or isotropically arranged with a uniform distribution between 90 and +90. 5. Fibers are generated to fill the space to the desired fiber volume fractions up to 80% (corresponding to 20% void volume fraction). 6. Fibers are assumed to be linear elastic and isotropic having Young’s modulus of 1 GPa and Poisson’s ratio of m = 0.3 and they are capable of carrying both tension and compression.
(measured from the center axis) is less than their diameter, otherwise they are not connected. More details on fiber generation are given in Appendix A. In this paper the effects of scale and fiber aspect ratio (varying fiber length while keeping the test box and fiber diameter fixed), fiber orientation and void volume fraction on the constitutive response of fiber networks are studied. The parameters are as follows: lx, ly, and lz are the lengths of the test box along the x, y, and z axes, respectively (Fig. 1), l is the fiber length, r is the radius of the fiber, hmax is the maximum orientation angle between the xaxis and fiber’s axis where the fibers are oriented between the positive hmax and negative hmax, and vfvoid is the void volume fraction. Given the fiber length l, a dimensionless window size (or size of SVE) can be defined as the ratio of the box length over the fiber length, and denoted as d
d¼
length of window lx ¼ length of fiber l
The parameter d, takes on values 1, 2, 3 and 5, for the studied fiber lengths (since the length of window is 10 units and fiber lengths are 2, 3.33, 5, and 10 units). The parameter used for the fiber orientation study is the orientation angle hmax, and the parameter used for the study of the effect of the void volume fraction is the void volume fraction vfvoid. Following [22], the apparent elastic properties are dependent on the size of a window (or specimen size) and boundary conditions, and this gives rise to a hierarchy of bounds
C R ðS R Þ
1
hS t1 i1 6 hS td0 i1 6 hS td i1 6 C eff 6 hC dd i 6 hC dd0 i
6 hC d1 i C V
ð2Þ
where "d0 < d and the inequality between any two tensors is understood as
C P D () ðC DÞ : a : a P 0;
for any tensor aij – 0
Here C is the fourth order stiffness tensor Cijkl and S is the compliance tensor Sijkl where S1 = C. The superscripts R and V denote Voigt and Reuss bounds, respectively, while the superscripts t and d imply traction and displacement boundary conditions, respectively. Following Eq. (2), the effective properties are bounded from above and below by the apparent elastic moduli obtained using displacement and traction boundary conditions, respectively. The larger is the window size d, the closer are the bounds. When d reaches the size of the RVE, the bounds merge and the effective properties are obtained. Eq. (2) directly applies to fiber networks in which the fiber length and diameter are fixed while the window size changes. In this paper, we fixed the window size and the fiber diameter while varying the fiber length. Those two problems are not the same for a given d since the fiber aspect ratio (fiber length:diameter ratio) is different, leading to a different microstructure. Thus, this paper
y-axis 5 2
Random fiber networks are created in a 3D test box by a Monte Carlo method. For each realization of the network, a different seed number is chosen and the locations of the centers of the fibers and their individual in-plane orientations are generated from a Poisson point field [26]; twenty realizations are generated for each data point. Thus, each generated fiber lies in plane (i.e., is parallel to the two largest surfaces of the test box). Fibers are linked by rigid connections when the nearest distance between two fibers
ð1Þ
7 4 6
3 Test box
z-axis 1 x-axis
Fig. 1. Test box and regions of trimming when the other end of the fiber is in other boxes.
717
Y. Lee, I. Jasiuk / Computational Materials Science 79 (2013) 715–723
E 2ð1 þ mÞ
investigates a combined effect of scale and microstructure, namely as d increases the fiber aspect ratio decreases.
G¼
2.2. Preprocessing in MATLAB
where m is the Poisson ratio, E is the elastic modulus, G is the shear modulus, and r is the radius of the fiber cross-section.
2.2.1. Input parameters The coordinates of the nodes and the elements which constitute the network are imported into a MATLAB program [28]. The information on the fiber network generation and the input parameters is given in Section 2.1. The input parameters are also summarized in Table 1. 2.2.2. Trimming of fibers In the fiber generation process, some fibers might protrude or fall outside of the test box. Therefore it is necessary to trim the fibers so that all of them are contained in the box. Fig. 1 shows the regions where the fibers might fall or protrude into. The procedure that ensures that the fibers do not protrude into negative regions of the box is summarized in the Appendix A. 2.2.3. Removing elements and clusters not contributing to the constitutive response When a loading is applied to the network at external boundaries, some of the randomly generated fibers in the test box may not experience the load, resulting in an overestimated volume fraction of the fibers. This will lead to underestimation of the constitutive response of the fiber network. Thus, for our calculations we remove the elements which do not contribute to the overall response of the network, except when we investigate this effect. The excluded elements include an element or a cluster, which does not experience any load (or is known as dangling) and those that experience loading only at one point. This concept has been verified through Abaqus; [1] whereby, for each condition stated above, the strain energy of the removed elements is negligible compared with the overall strain energy of the fiber network system.
2.3.1. Choosing a beam element There are two types of beams that can be chosen in Abaqus, namely the Euler–Bernoulli beam and the Timoshenko beam. The difference between these two beam elements is that there is a second-order spatial derivative present in the Timoshenko beam, resulting in the lower stiffness of the beam and larger deflection under loading. It is also well known that the Timoshenko beam is a well posed model and thus able to handle stubby beams such as short fiber segments in the network. Thus, we select it over the Euler–Bernoulli beam for modeling fibers. The shear coefficient for a circular cross-section is defined as [24]:
6ð1 þ mÞ2 7 þ 12m þ 4m2
SR2 ¼ kGpr 2 ;
2.3.3. Specifying boundary conditions We apply the in-plane uniaxial and shear displacements (with uniform strains) on the six faces of the test box in order to obtain the apparent stiffness coefficients C1111 and C1212, respectively, of the studied fiber networks. We denote ui as a displacement, xi as a global coordinate, eij as a uniform strain, with subscripts i and j taking values 1, 2, 3, which correspond to x, y, and z directions, respectively. To obtain C1111, the uniaxial displacement (with a uniform normal strain e11) is applied as follows on all six boundaries:
u1 ¼ e11 x1 u2 ¼ 0
For determining C1212, the applied in-plane shear displacement (with a uniform shear strain e12) is taken as:
u1 ¼ e12 x2 u2 ¼ e12 x1
SR3 ¼ 0:25
r11 1 0 C1111 Br C BC B 22 C B 2211 C B B B r33 C B C3311 C B B Br C ¼ BC B 23 C B 2311 C B B @ r31 A @ C3111 r12 C1211 0
C1122
C1133
C1123
C1131
C2222
C2233
C2223
C2231
C3322
C3333
C3323
C3331
C2322 C3122
C2333 C3133
C2323 C3123
C2331 C3131
C1222
C1233
C1223
C1231
e11 1 C B C2212 C CB e22 C CB C C B C3312 e C CB 33 C C B C2312 C CB 2e23 C CB C C3112 A@ 2e31 A C1212 2e12 C1112
10
ð8Þ Since only either normal or shear strains (e11or e12) are applied (in the x–y plane), the constitutive Eq. (8) reduce to the following equations:
r11 ¼ C1111 e11
ð9Þ
r12 ¼ 2C1212 e12
ð10Þ
respectively, where r11 and r12 are normal and shear stresses, respectively. The elastic strain energy W of the system is given by:
ð4Þ
V rij eij 2
ð11Þ
By incorporating Eq. (9) for the uniaxial case, and assuming that the system experiences a global uniform strain, C1111 can be found using Eqs. (12)–(14).
Table 1 Parameters used for generating fiber networks.
Values
ð7Þ
u3 ¼ 0
W¼
Parameters
ð6Þ
u3 ¼ 0
ð3Þ
The shear rigidity (SR) is then imported into Abaqus using the following definitions:
SR1 ¼ kGpr 2 ;
2.3.2. Specifying beam orientations Three dimensional (3D) beam elements are used in the simulations. The orientation of a beam is defined in terms of a local, right handed system. An additional node of the beam axis is defined and included in the element’s connectivity list. Abaqus is then able to compute three vectors through these nodes. The dot and cross product of these vectors give an orientation of the beam.
2.3.4. Calculating C1111 and C1212 stiffness constants The constitutive equations relate stresses to strains as follows:
2.3. Importing into Abaqus
k¼
ð5Þ
Box dimension (lx ly lz)
Radius of fiber r
vfvoid (%)
10 10 1
1/6
10, 20, 40, 60, 80
hmax ()
d
10, 90
1, 2, 3, 5
W¼
V C 1111 e211 2
V ¼ lx ly lz
ð12Þ ð13Þ
718
Y. Lee, I. Jasiuk / Computational Materials Science 79 (2013) 715–723
W¼
N X ei
ð14Þ
i¼1
where V is the volume of the test box, ei is the strain energy of a particular element, and N is the total number of elements in the system. For the shear loading case, C1212 can be found from Eqs. (13)–(15).
W ¼ 2VC 1212 e212
ð15Þ
All the information necessary to calculate C1111 and C1212 is obtained from the output file of Abaqus postprocessing module. For simplicity, the value C1111 is normalized by Cfiber 1111 and C1212 is norfiber fiber malized by Cfiber 1212 , where C1111 and C1212 are the linear elastic and isotropic stiffness components of the fiber. They can be calculated through the stress–strain relationship for an isotropic material. Cfiber 1111 is found from e.g. [8]:
r11
E ¼ ½ð1 mÞe11 þ mðe22 þ e33 Þ ð1 þ mÞð1 2mÞ
ð16Þ
where E, m are the Young’s modulus and Poisson’s ratio. After simplification, using Eq. (16), Cfiber 1111 is given as follows:
Cfiber 1111 ¼
Eð1 mÞ ð1 þ mÞð1 2mÞ
ð17Þ
Cfiber 1212 is the same as the shear modulus of an isotropic fiber, given by Eq. (5). 3. Results and discussion Fiber networks with the preferential and isotropic orientations (one realization for each case) under a uniaxial displacement are shown in Fig. 2(a) and (b), respectively. Note that in our model the fibers can resist compression. Displacements of the fiber endpoints are shown by arrows (in orange in color). Next, the results on the effects of scale and fiber aspect ratio (plotted in terms of window size d), the fiber orientation, and the void (or fiber) volume fraction on the stiffness components C1111 and C1212 are presented. Twenty realizations were considered for each data point. Fig. 3 shows the plot of the normalized apparent stiffness coefficient C1111 =C fiber 1111 versus window size d for fiber networks with preferentially oriented fibers (10 and +10) with the void volume fraction of 20% (i.e. 80% fiber volume fraction), obtained by applying a uniaxial displacement with a uniform strain e11. The normalized apparent stiffness coefficient C1111 =Cfiber 1111
Fig. 3. Non-dimensionalized stiffness component C1111 =Cfiber 1111 versus window size d for the fiber network with a preferential fiber orientation 10 and +10 for 20% void volume fraction (or 80% fiber volume fraction) with standard deviation bar. Twenty realizations were calculated for each window size. In calculations test box size and fiber diameter was fixed while fiber length was varied. Thus, as the window size increases the fiber aspect ratio decreases.
decreases with the increasing d (and decreasing fiber aspect ratio). The fact that the stiffness decreases with increasing window size is consistent with [22]. Following [22], when d becomes large, leading to the RVE size, the normalized apparent stiffness coefficient C1111 =Cfiber 1111 converges from above to an effective stiffness. Recall that the studied problem involves a change in the fiber aspect ratio with the change in d. Thus, our simulations include the effect of scale d combined with changes in the fiber aspect ratio. Also, the scatter in the normalized C1111 =Cfiber 1111 decreases as d increases as expected [32]. The application of traction boundary conditions with a uniform stress would give a lower bound for the apparent elastic stiffness. However, several issues arise when applying the traction boundary conditions to fiber networks. Since fibers are randomly oriented the cross-sectional areas of beam elements are different for each fiber at the surface (boundary) of the box. Also, it is impossible to unambiguously apply a constant stress on the boundary since fibers cross the box boundary at different angles and locations [32,34]. Therefore, we do not apply this second type of boundary conditions. Also, for simplicity of notation we do not include superscripts d. Figs. 4 and 5 show the probability plots to fit the data for C1111 of fiber networks with preferentially oriented fibers (10 and
Fig. 2. Beam network loaded under a uniaxial displacement boundary condition (a) with the preferential fiber orientation (10 and +10), (b) with the isotropic fiber orientation (90 and +90). Displacements of endpoints are shown with arrows (in color in orange). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Y. Lee, I. Jasiuk / Computational Materials Science 79 (2013) 715–723
+10) with 20% void volume fraction (80% fiber volume fraction) for window sizes 1 and 5. Similar curves were obtained for the window sizes 2 and 3. Ten distribution functions: Beta, Gamma, Logistic, Gumbel Min, Rayleigh, Weibull, v and other (see Appendix B) were considered and the Kolmogorov–Smirnov test was used to determine the best distribution function fits for each of the four window sizes. We find that the probability distribution that best fits the data is the Logistic, followed by Gamma and Beta. Note that the commonly used Weibull distribution does not give a good fit to the data. Next, we consider fiber networks with isotropically arranged fibers (90 and +90) and investigate how d affects the apparent constitutive response under a uniform uniaxial strain (Fig. 6). The trends are similar to those for networks with preferential fiber orientations (Fig. 3). Again, the normalized apparent stiffness component C1111 =Cfiber 1111 decreases with increasing d but at lower rate than for the preferentially oriented fiber case. The constitutive response for the isotropically distributed fiber network is lower than that of the preferentially oriented fiber network, in particular for small d. Again, this is expected since in such networks the fibers, which are arranged preferentially in the loading direction and have large aspect ratios, better transfer loads than those in the isotropic networks. We also observe that the number of elements removed in order to achieve a given value of the void (or fiber) volume fraction stay about the same for an isotropic network for different window
719
sizes while the number of the elements removed in a preferentially oriented network increases as the window size increases. Next, the effects of the void (or fiber) volume fraction on the stiffness constant C1111 =Cfiber 1111 are investigated, as shown in Figs. 7 and 8 for the networks with isotropic fiber orientation for d = 2 and 3 and in Fig. 9 for the networks with preferential fiber orientation for d = 2. As expected, the stiffness component decreases with increasing void volume fraction. The initial decrease is linear indicating affine deformations. For the isotropic case the decrease is linear up to vfvoid of about 80% when d = 2 and up to about vfvoid of 60% for d = 3. Note C1111 =Cfiber 1111 does not reach zero because our formulation ensures that all the structural elements are connected to clusters which can sustain loads. Moreover, the scatter in constitutive values decreases with increasing void volume fraction. Also, the stiffness component C1111 =Cfiber 1111 of the fiber network with preferentially oriented fibers is higher than that for an isotropic fiber network as expected since the fibers aligned in the loading direction can carry the load more efficiently. The same study as for C1111 was performed to investigate the effect of scale (with changing fiber aspect ratio), void volume fraction and fiber orientation on the normalized apparent shear stiffness component C1212 =Cfiber 1212 as illustrated in Figs. 10–14. We observe that the beam network with isotropic fiber orientation can withstand a shear loading better than the preferentially oriented beam network, as shown by comparing Figs. 12 and 14 for the window
Fig. 4. Fits for C1111 by several distribution functions for the case of d = 1 (for preferential fiber orientation).
Fig. 5. Fits for C1111 by several distribution functions for the case of d = 5 (for preferential fiber orientation).
720
Y. Lee, I. Jasiuk / Computational Materials Science 79 (2013) 715–723
Fig. 6. C1111 =Cfiber 1111 versus window size d for a fiber network with an isotropic fiber orientation 90 and +90 for 20% void volume fraction (with standard deviation bar).
Fig. 7. C1111 =Cfiber 1111 versus void volume fraction of beam network with an isotropic fiber orientation and d = 2 (with standard deviation bar).
Fig. 9. C1111 =Cfiber 1111 versus void volume fraction of beam network with a preferential fiber orientation and d = 2 (with standard deviation bar).
Fig. 10. Non-dimensionalized stiffness component C1212 =Cfiber 1212 versus window size for a fiber network with a preferential fiber orientation for 20% void volume fraction (with standard deviation bar).
Fig. 8. C1111 =Cfiber 1111 versus void volume fraction of beam network with an isotropic fiber orientation and d = 3 (with standard deviation bar). Fig. 11. C1212 =Cfiber 1212 versus window size for a fiber network with an isotropic fiber orientation for 20% void volume fraction (with standard deviation bar).
size 2. In addition, the preferentially oriented fiber network decays faster than the isotropic fiber network as seen by comparing Figs. 10 and 11. Despite the above mentioned differences, the volume fraction studies (Figs. 12 and 13) of C1212 share similarity with C1111. The stiffness C1212 =Cfiber 1212 decreases linearly with increasing void volume fraction indicating affine deformations. For the isotropic case the decrease is linear up to vfvoid of about 80% when d = 2 and up to about vfvoid of 60% for d = 3 as for the C1111 case.
Next, we investigate the effect of dangling fibers, i.e. fibers that do not carry loads in the network. In order to verify the removal of elements which do not sustain loads, the strain energy densities of the network with dangling fibers and the network with removed elements are compared. A check on the fiber network confirmed that there were no fibers that experience zero strain energy. There are various ways of defining the void (or fiber) volume fraction in a network. When a network is generated, whether it
721
Y. Lee, I. Jasiuk / Computational Materials Science 79 (2013) 715–723
Fig. 12. C1212 =Cfiber 1212 versus void volume fraction of beam network with an isotropic fiber orientation and d = 2 (with standard deviation bar).
Fig. 13. C1212 =Cfiber 1212 N versus void volume fraction of beam network with an isotropic fiber orientation and d = 3 (with standard deviation bar).
the preferentially oriented fiber network in Table 4. For the isotropic fiber network, the normalized C1111 value, fiber volume fraction and percentage of elements removed do not vary significantly (Table 3). This is to be expected since the elements in the isotropic fiber network are more easily connected to each other when compared with the elements in a preferentially oriented fiber network. Therefore, there are fewer dangling fibers in the isotropic fiber network leading to relatively constant values of parameters under investigation. For the fiber network with the preferentially oriented fibers, Cases II and III have the similar normalized C1111 value, with discrepancy of about 2%, compared with a large increase for Case I (Table 4). Note that Case II is Case III plus dangling fibers, and since dangling fibers do not carry any load, we should expect both cases to have similar normalized C1111 value. As for Case III, it has a higher normalized C1111 value since all the fibers are interconnected and contribute to the load sharing capacity of the network. Besides, Case III has a lower fiber volume fraction compared with the other two cases since approximately 16% of dangling fibers have been removed. The above analysis gives three scenarios for different definitions of fiber volume fraction of a network. Occasionally, Case II is being used in modeling of fiber networks. However, the volume fraction might be overestimated since the actual volume fraction might be lower due to dangling fibers. In addition, the percentage of dangling fibers in the network is a variable that has to be taken into account if the isotropic network is to be generated. Therefore, Case I gives a more realistic representation of a network. Another issue that needs to be considered is the effect of an intersected volume on the calculated volume fraction of fibers. The fiber volume fraction vf can be estimated as follows:
Table 2 Different definitions of the volume fraction of a network. Note that Case I was used in data reported in Figs. 3–14. Case I
Case II
Case III
1. Generate the fiber network until 80% fiber volume fraction 2. Remove any dangling fibers 3. Add more fibers and repeat step 2 until fiber volume fraction reaches 80%
1. Generate the network until 80% fiber volume fraction 2. Do not remove any dangling fiber
1. Generate the network until 80% fiber volume fraction 2. Remove any dangling fibers 3. Do not add more fibers
Table 3 Effect of different methods of generating fiber network with an isotropic orientation with a 20% void fraction (80% fiber fraction) for d = 2 on C1111 =Cfiber 1111 . Fig. 14. C1212 =Cfiber 1212 versus void volume fraction of beam network with a preferential fiber orientation and d = 2 (with standard deviation bar).
is through a manufacturing process as in an air filter case, or through a computer simulation, one may find dangling fibers not contributing to load sharing but occupying a definite amount of space in the network. This would overestimate the volume fraction in the network, leading to incorrect representation of the load sharing capability of the network. In this paper, we consider three Cases (I, II and III) that can be used to define the fiber volume fraction of a network, summarized in Table 2. In this analysis the 20% void (80% fiber) volume fraction and d = 2 are used. We illustrate the implication of these three different Cases (I–III) on C1111 for the isotropic fiber network in Table 3 and for
Cases
C1111 =Cfiber 1111
Fiber volume fraction (%)
Elements removed (%)
Case I Case II Case III
0.0968 0.0966 0.0965
80 80 80
0.545 0 0.547
Table 4 Effect of different methods of generating fiber network with a preferential fiber orientation with a 20% void fraction (80% fiber fraction) for d = 2 on C1111 =Cfiber 1111 . Cases
C1111 =Cfiber 1111
Fiber volume fraction (%)
Elements removed (%)
Case I Case II Case III
0.125 0.0792 0.0778
80 80 62.6
8.96 0 16.83
722
vf
¼
Y. Lee, I. Jasiuk / Computational Materials Science 79 (2013) 715–723 n X
pr2 li
ð18Þ
i¼1
where r is the radius and l is the length of each fiber, respectively. This calculation does not take into account the intersected volume at each joint. Thus, the total volume fraction of fibers is exaggerated and therefore the constitutive response is underestimated. This effect is more pronounced when the fibers are aligned, almost parallel to each other, as in the preferential fiber orientation case. Furthermore, when the window sizes get smaller, or the fibers get shorter compared with the length of a test box, the effect of the intersected volume is higher as well. Therefore, there is a need to calculate these intersected volumes in order to obtain more accurate results. More details on the determination of the intersected volume are given in [27]. In summary, we studied numerically the effect of window size (and fiber aspect ratio), fiber orientation and void volume fraction on the apparent elastic stiffness components C1111 and C1212 of fiber networks with isotropic and preferential fiber orientations. In our parametric study we considered only several sets of inputs to investigate the trends; more comprehensive parametric studies can be done in the future. Also, we combined the effect of scale and fiber aspect ratio while these effects could also be studied separately. Also, a nonlinear fiber response can be incorporated in the Abaqus in order to capture nonlinear effects in the constitutive response of fiber networks. Strength of the network can also be investigated through element removal and fiber–fiber debonding once fibers reach a critical stress level. Last but not least, the entire stiffness tensor can be found through applications of additional loadings. The current study can be extended to biological materials, such as an articular cartilage which is known to have three distinct morphological zones with different fiber orientations in each zone [45]. The knowledge how fiber properties and structural characteristics of fiber networks affect their constitutive response is important in studies of healthy and diseased biological tissues and variety of synthetic fibrous materials. 4. Summary We modeled random fiber networks as beam networks and predicted two apparent stiffness components C1111 and C1212 as a function of several parameters that affect their means and scatter. This problem was solved computationally using a finite element method. Two fiber orientations were considered, isotropic and preferential, for various void volume fractions and scales of observation (and fiber aspect ratio). The effect of dangling fibers was also investigated and was shown to play a role in calculations, in particular, for the preferentially oriented network. The calculation of fiber volume fraction was also briefly discussed. Acknowledgement We gratefully acknowledge support of the National Science Foundation (CMMI 0927909 – Dr. Ken Chong). Appendix A In the fiber generation process, some fibers might protrude or fall outside of the test box. Therefore it is necessary to trim the fibers so that all of them are contained in the box. Fig. 1 shows the regions where the fibers might fall or protrude into. The following formulations ensure that the fibers do not protrude into negative regions of the box. For the region 1, the following equations are used:
Y 2new ¼
Z 2new ¼
lx X 1old Y 2old Y 1old
þ Y 1old
ðA:1Þ
þ Z 1old
ðA:2Þ
X 2old X 1old lx X 1old Z 2old Z 1old X 2old X 1old
X 2new ¼ lx
ðA:3Þ
where X, Y, Z are the local coordinates of the fiber endpoints (node), the superscripts 1 and 2 denote the node inside the test box and the node outside the test box, respectively, while the subscripts old and new denote the coordinates before and after trimming, respectively. For the region 2, the following equations are used:
X 2new ¼
Z 2new ¼
ly Y 1old X 2old X 1old Y 2old Y 1old ly Y 1old Z 2old Z 1old Y 2old Y 1old
þ X 1old
ðA:4Þ
þ Z 1old
ðA:5Þ
Y 2new ¼ ly
ðA:6Þ
For the region 3, the following equations are used:
X 2new ¼
Y 2new ¼
lz Z 1old X 2old X 1old Z 2old Z 1old lz Z 1old Y 2old Y 1old
Z 2new ¼ lz
Z 2old Z 1old
þ X 1old
ðA:7Þ
þ Y 1old
ðA:8Þ
ðA:9Þ
For region 4, where we do not know a priori whether a fiber will protrude through region 1 or region 2, we determine it in the following way. We first assume that the fiber will protrude into region 4 through region 1. Eq. (A.1) is then used to find the coordinate Y of the new node. The coordinate Y of this node is then compared with ly. If ly is larger than the new coordinate Y, then the fiber is confirmed to protrude from the test box into region 4 through region 1. Therefore, Eqs. (A.1)–(A.3) are used to find the node of intersection on the face of the test box. On the contrary, if ly is smaller than the new coordinate Y, then Eqs. (A.4)–(A.6) are used to find the new coordinates since the fiber will protrude into region 4 through region 2. The same concept can be applied to fibers protruding into region 5. Eqs. (A.4)–(A.6) are used if the fibers will protrude through region 2, and Eqs. (A.7)–(A.9) are used if the fibers will protrude through region 3. Similarly for region 6, Eqs. (A.1)–(A.3) are used for the fibers protruding through region 1, and Eqs. (A.7)–(A.9) are used for the fibers protruding through region 3. The region 7 requires a more careful treatment since a fiber will generally go through two regions before reaching region 7. First, let us assume that the fiber will protrude from the test box into region 7 through region 1, and, therefore, Eqs. (A.2) and (A.3) are used. If both ly and lz are larger than the new coordinate Y and the new coordinate Z, respectively, then Eqs. (A.1)–(A.3) are used since the fiber will go through region 1 and either region 4 or region 6 before reaching region 7. However, if either ly or lz is smaller than the new coordinate Y or new coordinate Z, respectively, then the fiber does not protrude through region 1, but through region 2 or region 3 instead. In order to determine the protruded region, let us assume that the fiber will protrude through region 2 and therefore Eqs. (A.4) and (A.5) are used. If lz is larger than the new
Y. Lee, I. Jasiuk / Computational Materials Science 79 (2013) 715–723
coordinate Z, then the fiber protrudes through region 2 and then through either region 4 or region 5 before reaching region 7. Therefore, Eqs. (A.4)–(A.6) are used. On the contrary, Eqs. (A.7)–(A.9) are used for the fiber protruding through region 3. Note that some elements might have two nodes at the same place if the elements protrude out of the test box through its faces. In order to avoid this situation, the length of each element is checked so that it is greater than zero, meaning that one node does not overlap with another. This trimming action can be used to study the constitutive response of a certain part of the fiber network as well. The test box can be placed in a specific area of the network, while the rest of the area can be trimmed off.
Appendix B For reference, probability distribution functions used in this study are listed. Beta
f ðx; a1 ; a2 ; d1 ; d2 Þ ¼
a1 1 a2 1 1 1 dxd 2 d1
xd1 d2 d1
ðd2 d1 ÞBða1 ; a2 Þ
ðB:1Þ
Chi
f ðx; a; b; eÞ ¼
a1 a 1 xe 2 212 x e a e2ð b Þ b bC 2
ðB:2Þ
k½kðx eÞj1 e½kðxeÞ CðjÞ
ðB:3Þ
Gamma
f ðx; j; k; eÞ ¼
Gumbel Max
f ðx; l; aÞ ¼ aeðaðlxÞe
aðlxÞ Þ
ðB:4Þ
Gumbel Min
f ðx; l; aÞ ¼ aeðaðxlÞe
aðxlÞ Þ
ðB:5Þ
Logistic
f ðx; a1 ; a2 Þ ¼
1 4a2 cosh
2 xa1 2a2
ðB:6Þ
Normal
1 1 xl 2 f ðx; l; rÞ ¼ pffiffiffiffiffiffiffi e2ð r Þ 2pr
ðB:7Þ
Rayleigh
f ðx; a; eÞ ¼
xe
a2
1 xe 2
e2ð a Þ
ðB:8Þ
Rectangular
f ðx; a1 ; eÞ ¼
1 a1
ðB:9Þ
Weibull
f ðx; l; j; eÞ ¼
xe le
j1
le
xe
j
jeðleÞ
ðB:10Þ
723
References [1] Abaqus, v. 6.7, ed., Providence, Rhode Island, Dassault Systemes Simulia Corp. [2] B. Alberts, A. Johnson, J. Lewis, M. Raff, K. Roberts, P. Walter, Molecular Biology of the Cell, second ed., Garland Science, New York, NY, 2002. [3] M. Alkhagen, S. Toll, Journal of Applied Mechanics – Transactions of the ASME 74 (2007) 723–731. [4] R.D. Averett, B. Menn, E.H. Lee, C.C. Helms, T. Barker, M. Guthold, Biophysical Journal 103 (2012) 1537–1544. [5] P.L. Chandran, V.H. Barocas, Journal of Biomechanical Engineering – Transactions of the ASME 128 (2006) 259–270. [6] X. Cheng, C. Wang, A.M. Sastry, S.B. Choi, Journal of Engineering Materials and Technology 121 (1999) 515–523. [7] H.L. Cox, British Journal of Applied Physics 3 (1952) 72–79. [8] R.R.J. Craig, in: Mechanics of Materials, second ed., Wiley, New York, 2000. [9] M. Delince, F. Delannay, Acta Materialia 52 (2004) 1013–1022. [10] B.A. DiDonna, T.C. Lubensky, Physical Review E 72 (2005) 066619. [12] M.L. Gardel, F. Nakamura, J.H. Hartwig, J.C. Crocker, T.P. Stossel, D.A. Weitz, Proceedings of the National Academy of Sciences of the United States of America 103 (2006) 1762–1767. [13] M.L. Gardel, J.H. Shin, F.C. MacKintosh, L. Mahadevan, P. Matsudaira, D.A. Weitz, Science 304 (2004) 1301–1305. [14] M.M. Giraud-Guille, Calcified Tissue International 42 (1988) 167–180. [15] M.M. Giraud-Guille, L. Besseau, R. Martin, Journal of Biomechanics 36 (2003) 1571–1579. [16] E. Hamed, Y. Lee, I. Jasiuk, Acta Mechanica 213 (2010) 131–154. [17] H. Hatami-Marbini, R.C. Picu, Analytical approach to quantifying the nonaffine behavior of fiber networks, in: Materials Research Society Symposium Proceedings, vol. 1060, 2008a, pp. 33–37. [18] H. Hatami-Marbini, R.C. Picu, Physical Review E 77 (2008) 062103. [19] H. Hatami-Marbini, R.C. Picu, Acta Mechanica 205 (2009) 77–84. [20] D.A. Head, A.J. Levine, F.C. MacKintosh, Physical Review E 68 (2003) 061907. [21] C. Hellmich, J. Barthélémy, L. Dormieux, European Journal of Mechanics, A/ Solids 23 (2004) 783–810. [22] C. Huet, Journal of Mechanics and Physics of Solids 38 (1990) 813–841. [23] E.M. Huisman, T. van Dillen, P.R. Onck, E. van der Giessen, Physical Review Letters 99 (2007) 208103. [24] J.R. Hutchinson, Journal of Applied Mechanics – Transactions of the ASME 68 (2001) 87–92. [25] I. Jasiuk, M. Ostoja-Starzewski, Biomech Model Mechanobiology 3 (2004) 67–74. [26] D. Jeulin, Random Structure Models for Homogenization and Fracture Statistics, in: D. Jeulin, M. Ostoja-Starzewski (Eds.), Mechanics of Random and Multiscale Microstructures CISM Courses and Lectures, Springer, Berlin, Heidelberg, New York, 2001, pp. 33–91. [27] Y.H. Lee, Non-Invasive Characterization of Bone Tissue and Modeling of Bone at Micron and Sub-Micron Levels, in: Mechanical Science and Engineering, University of Illinois, Urbana, 2010, p. 175. [28] MATLAB, TheMathWorks, Natick, Massachusetts, 2008. [29] M.A. Narter, S.K. Batra, D.R. Buchanan, Proceedings of the Royal Society London, Series A – Mathematical Physical and Engineering Sciences 455 (1999) 3543–3563. [30] P.R. Onck, T. Koeman, T. van Dillen, E. van der Giessen, Physical Review Letters 95 (2005) 178102. [31] M. Ostoja-Starzewski, Applied Mechanics Reviews 55 (2002) 35–60. [32] M. Ostoja-Starzewski, Microstructural Randomness and Scaling in Mechanics of Materials, Chapman & Hall/CRC/Taylor & Francis, Boca Raton, FL, 2008. [33] M. Ostoja-Starzewski, D.C. Stahl, Journal of Elasticity 60 (2000) 131–149. [34] M. Ostoja-Starzewski, C. Wang, Acta Mechanica 80 (1989) 61–80. [35] M. Petrty´l, J. Hert, P. Fiala, Journal of Biomechanics 29 (1996) 161–169. [36] A.M. Sastry, C.W. Wang, L. Berhan, Key Engineering Materials 200 (2001) 229– 250. [37] M.H. Schwartz, P.H. Leo, J.L. Lewis, Journal of Biomechanics 27 (1994) 865–873. [38] P.P. Shang, Characterization of Anisotropy of Spunbonded Fiberwebs (nowovens), North Carolina State University, 1985. [39] D.C. Stahl, S.M. Cramer, ASME Journal of Engineering Materials and Technology 120 (1998) 126–130. [40] T. Stylianopoulos, V.H. Barocas, Computer Methods in Applied Mechanics and Engineering 196 (2007) 2981–2990. [41] C. Wang, X. Cheng, A.M. Sastry, S.B. Choi, Journal of Engineering Materials and Technology 121 (1999). 530–513. [42] C.W. Wang, L. Berhan, A.M. Sastry, Journal of Engineering Materials and Technology 122 (2000) 450–459. [43] C.W. Wang, A.M. Sastry, Journal of Engineering Materials and Technology 122 (2000) 460–468. [44] J.W. Weisel, Biophysical Chemistry 112 (2004) 267–276. [45] J.Z. Wu, W. Herzog, Journal of Biomechanics 35 (2002) 931–942. [46] X.-F. Wu, Y.A. Dzenis, Journal of Applied Physics 98 (2005) 093501. [47] Y.J. Yoon, S.C. Cowin, Biomechanics and Modeling in Mechanobiology 7 (2008) 1–11. [48] T.I. Zohdi, Computational Methods in Applied Mechanics and Engineering 196 (2007) 2972–2980.