- Email: [email protected]
Analysis of the compressive response of Nano Fibrillar Cellulose foams
Mechanics of Materials 80 (2015) 13–26
Contents lists available at ScienceDirect
Mechanics of Materials journal homepage: www.elsevier.com/locate/mechmat
Analysis of the compressive response of Nano Fibrillar Cellulose foams Prashanth Srinivasa 1, Artem Kulachenko ⇑ KTH Royal Institute of Technology, School of Engineering Sciences, Department of Solid Mechanics, SE-100 44 Stockholm, Sweden
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 8 June 2014 Received in revised form 1 September 2014 Available online 5 October 2014
Nano Fibrillar Cellulose (NFC) is fast emerging as a biomaterial with promising applications, one of which is cellular foam. The inner structure of the foam can take various shapes and hierarchical micro-structures depending on the manufacturing parameters. The compressive response of foams developed from these materials is currently a primary criterion for the material development. In this work, we focus on the connection between the non-linear part of the response and the inner structure of the material. We study the effect of internal contact and its contribution to gradual stiffening in the energy absorbing region and accelerated stiffening in the densification region of the large strain compressive response. We use the finite element method in this study and discuss the applicability and efficiency of different modelling techniques by considering well defined geometries and available experimental data. The relative contribution of internal contact is singled out and mapped onto the overall compressive response of the material. The effect of initial non-straightness of the cell walls is studied through superposing differing percentages of the buckling modes on the initial geometry. The initial non-straightness is seen to have a significant effect for only strains up to 1%. The secant modulus measured at slightly higher strains of 4%, demonstrates lesser effect from the non-straightness of cell walls. The simulations capture the compressive response well into the densification regime and there is an order of magnitude agreement in between simulations and experiments. We observed that internal contact is crucial for capturing the trend of compressive response. Ó 2014 Elsevier Ltd. All rights reserved.
Keywords: NFC foams Internal contact Voronoi structures Foam densification Effect of foam porosity
1. Introduction Low density cellular materials like foams have found applications in a wide range of fields. They are used as packaging, shock-absorbing and insulation material. Polymeric foams are leading the way in such applications. Concurrently, a new nature-based material, cellulose nanofibrils (Nano Fibrillar Cellulose, NFC), is fast emerging as an
⇑ Corresponding author. Tel.: +46 87908944. E-mail addresses: [email protected] (A. Kulachenko). 1 Tel.: +46 707948887.
(P.
Srinivasa),
http://dx.doi.org/10.1016/j.mechmat.2014.09.006 0167-6636/Ó 2014 Elsevier Ltd. All rights reserved.
[email protected]
attractive alternative to synthetic polymers with potential applications in biomedical, structural, chemical and numerous other areas (Klemm et al., 2011). The main advantage with using cellulose fibrils is their high stiffness-to-density ratio, wide availability and biocompatibility. Nanofibrils have already demonstrated their potential in reinforcing the cell walls of polymer foams (Svagan et al., 2008), in the manufacturing of high porosity foams from nanofibre suspensions (Sehaqui et al., 2010) and in aerogels which exhibit foam like structure (Aulin et al., 2010). Depending on the inner hierarchical structure and porosity of these materials, the mechanical properties
14
P. Srinivasa, A. Kulachenko / Mechanics of Materials 80 (2015) 13–26
can vary in a wide range. Naturally, providing insights on how the microstructure of foams affects the mechanical response is essential for material engineers as it provides opportunity to establish the mechanical property space of these newly emerging materials. The task of connecting the micro-structure of the material with mechanical response has been a topic of intense research. The mechanical properties of polymeric foam structures and other cellular solids were initially studied through the analysis of a single cell, e.g. hexagonal honeycombs (Gibson and Ashby, 1988). Such analyses proved to be very successful in predicting the elastic properties, but due to a complex nature of the compressive response associated with local cell wall buckling and inner contact, the non-linear response has been difficult to address. Furthermore, there is noticeable difference in response of the nanocellulose (NFC) foam structure under compressive loading as compared to that of the polymeric foams with more ordered cellular structure. The stress plateau or the energy absorbing zone of the compressive response in NFC foams sees gradual stiffening of the response (Fig. 1) as opposed to open cell polymeric foams which exhibit a clear plateau (Gong et al., 2005). This stiffening has been attributed to the better load bearing capability of the NFC cell walls (Sehaqui, 2011). Thus it is imperative that there are microstructural factors at play which should account for the distinct response of the NFC foams under compression. In this work, we will provide insights to the micromechanics of non-linear compressive behaviour of NFC foams. We will begin by modelling a two dimensional foam structure with random Voronoi tessellations. The aim is to establish a numerical model that will be able to capture the stiffening response in the energy absorbing zone as well as the densification at large strains. The suitability of beam element as compared to the plane strain solid element will be analysed for simple hexagonal honeycomb geometry. Finally, the effect of internal contact on the stress response is studied and its contribution to
the gradual stiffening in the energy absorbing and the stiffening region will be analysed. 2. Methods NFC foams are manufactured through the process of freeze drying. When the NFC suspension is freeze dried, nucleation process is initiated for the formation of ice crystals (Svagan et al., 2010). The geometry of the foam like structures is well captured by the Voronoi tessellations. A number of researchers (Roberts and Garboczi, 2002; Silva et al., 1995; Vanderburg et al., 1997; Zhu and Windle, 2002) have studied the Voronoi structure in both two- and three-dimensions. Zhu and Windle (2002) have also studied the effects of cell irregularity on the high strain compression. The stress–strain responses of the random foam structures in the studies above show a clear transition between the linear elastic regime and stress plateau. Furthermore, the stress plateau is truly a plateau with marginal increase in compressive stress levels before densification. However, as we indicated earlier, the experimental results of the compressive stress–strain response of NFC foams suggest that there is a gradual stiffening even in the so called ‘‘stress plateau’’ (Ali and Gibson, 2013; Sehaqui et al., 2010) and the transition between the linear response and stress plateau is not distinct. Studies based on SEM images of NFC foams suggest that the structure is made up of irregular sheets with occasional connectivity to other sheets (Ali and Gibson, 2013). The compressive response is reported to fit well with those predicted by theoretical estimates on open cell foams. The highly irregular nature of the foam structure makes it a good candidate for being modelled through Voronoi tessellations. The connectivity in the out-of-plane direction is non-trivial and might be a limiting factor for a 2D model. However, as the results from simulations will testify, it is an agreeable start. 2.1. Modelling of the Voronoi structure There are two aspects to the modelling of NFC foams as Voronoi tessellations. The first one corresponds to the generation of a periodic structure that lends itself to the application of periodic boundary conditions and secondly, the generation of finite random thickness cell walls for the periodic structure.
8 99.5% 98.1% 95.9% 93.1%
7
Stress [MPa]
6 5 4 3 2 1 0 0
10
20
30
40 50 Strain [%]
60
70
80
Fig. 1. Compressive response of NFC foams of different porosities (Plotted from experimental data provided by Houssine Sehaqui (Sehaqui, 2011)).
2.1.1. Periodic microstructure The generation of periodic structure rests on the property of the tessellations. We consider a set of random points inside a unit square in two-dimensional Euclidean space. Additional points are created around this unit square which are translations of the random points inside the unit square (Okabe et al., 1995). The translations of these points are governed by the transformations in Eqs. (1) and (2).
xo ¼ xi pL
ð1Þ
yo ¼ yi qL
ð2Þ
P. Srinivasa, A. Kulachenko / Mechanics of Materials 80 (2015) 13–26
15
Here, xo , yo are the nucleating points outside of the unit cell of dimensions L L; xi and yi are points within the unit cell and p; q are integers. This is demonstrated in Fig. 2 where the points inside the central square are translated all around the square unit cell. The Voronoi tessellation, of the resulting space which now has an area 9 times that of the unit cell, is shown in Fig. 3. It is seen that the central square is now periodic and can tile the whole two dimensional (2D) space. The tessellation consists of polygons, all of which are convex and have random number of edges. 2.2. Finite thickness cell walls The aspect of having finite random thickness cells is accomplished by dilating these polygons towards their centroidal points by using Eqs. (3) and (4):
xnew ¼ A x þ ð1 AÞ xc
ð3Þ
ynew ¼ A y þ ð1 AÞ yc
ð4Þ Fig. 3. Periodic Voronoi tessellations.
Here xnew and ynew corresponds to the vertices of the deflated polygon, x and y corresponds to the vertices of the initial polygon, xc and yc corresponds to the centroid of the initial polygon, and A corresponds to the factor of dilation. As can be seen, this result in variable thickness for each of the edges of the polygon as the deflation is not uniform. This is because each of the vertices which are at different distance from the centroid is scaled by one single factor. A binary subtraction of the dilated polygons from their respective parent polygons leaves us behind with the necessary finite cell wall thickness, random porous structure as shown in Fig. 4 The variation in porosity can be represented either by varying the thickness of cell walls or by having lesser number of cells in a given space. In this study, we will limit ourselves by varying porosity through altering the cell wall
Fig. 4. Periodic unit cell with finite thickness cell walls.
thickness. The thickness of the cell walls is varied through the deflation factor in order to obtain structures with different porosities. 2.3. Appropriate choice of finite element
Fig. 2. Translation of the tessellation points.
The model of random structures with high porosities consists of slender cell walls which undergo flexural deformation under compressive loads. It is essential to make a reasoned choice on the finite element used to discretise the structure. We intend to compare the effectiveness of 2D plane solid elements against beam elements. Plane elements are occasionally used for similar analyses as they can offer a better representation of the cell wall structure, particularly at the vertexes. We do the comparison through simulating the compression of pre-stressed aluminium foam with hexagonal honeycomb structure, for which we found well defined material data and comprehensive
16
P. Srinivasa, A. Kulachenko / Mechanics of Materials 80 (2015) 13–26
experimental results. The thickness of the aluminium honeycomb cell walls are comparable to that of the random structure that we wish to analyse, making it a good choice for the purpose of choosing a suitable finite element. Both experimental and numerical results have been previously reported in two successive papers (Papka and Kyriakides, 1994, 1998). 2.4. Hexagonal honeycomb structure The geometry consists of aluminium sheets which are bonded along its length at regular intervals as shown in Fig. 5(a). The sheets are then stretched by applying pressure resulting in the formation of hexagonal honeycomb Fig. 5(c). This process of stretching is meant to simulate the actual manufacturing process and contributes to the pre-stressed state in the cell walls. The applied pressure is then released to relax the structure to the new equilibrium state (Fig. 5(d)). Following the relaxation step, large strain compression is carried out through applied displacement on upper and bottom sides of the structure (Fig. 6). The numerical simulations are carried out using a commercial, implicit time integration finite element code (ANSYS). Two different elements are used for the purpose of comparison. First a three node finite strain beam element (BEAM189 in ANSYS) with quadratic displacement behaviour is considered. Each node has three degrees of freedom in translation and three degrees of freedom in rotation. The second element is a quadrilateral 2D plane strain solid element (PLANE183 in ANSYS) with 8 nodes including mid nodes. Each node has two degrees of freedom in translation and incorporates quadratic displacement behaviour. A bilinear isotropic hardening material model is used and the material data is presented in Table 1. The discretisation along the length of cell walls is gradually increased and it is found that convergence is achieved for 20 elements per length of the cell wall. The obtained response is stiffer and the point of first yield was higher than the experimental results. This is attributed to the shear locking phenomenon. Bending energy is proportional to the cube of thickness whereas shear energy is proportional to the first power of thickness. As the thickness
becomes increasingly small, the contribution from bending energy reduces more quickly than from shear energy. This of course can be overcome with mesh refinement and changing the aspect ratio of elements. The cell wall is now discretised along its thickness and it is found that a minimum of 2 elements per thickness is required to achieve mesh convergence. The results are shown in Fig. 7. Finally, a comparison is carried out between beam and solid elements with identical discretisation of twenty elements per length of the cell wall. The simulations were carried out on an Intel (R) Xeon (R) E7-8837 processor (Frequency: 2.66 GHz) and 8 cores were used. The data for computational effort required without and with internal contact is presented in Tables 2 and 3 respectively. The results from both the solid and beam elements along with the experimental results are presented in Fig. 8. The initial linear elastic regime and the final stiffening due to densification are represented equally well by both the beam and solid elements. In the stress plateau region, the beam elements capture localised buckling of the cell walls while the solid element predicts a flatter response. The overall agreement with the experimental results is exceptionally good. The computational effort in terms of the wall time is two times lower for the beam element as compared to the solid element when internal contact is not incorporated. This is surprisingly small difference taking into account the fact that the beam model has almost ten times smaller number of active degrees of freedom. The reason for that is we used the 3D beam elements with constrained out-of-plane degrees of freedom. (The reason 2D beam elements were not used was unavailability of the non-linear stabilization tools for these elements in ANSYS, which were critical for avoiding convergence problems passing over the buckling bifurcation points.) Furthermore, the total number of integration points was actually greater in the model with beam element since the selected beam element uses 5 integration points through thickness as opposed to 4 in the case of two plane elements through the thickness. When the internal contact is considered, the computation time for beam elements is comparable to that of the
Fig. 5. Creation of hexagonal honeycomb structure: (a), (b), (c) – expansion (d) – relaxation.
17
P. Srinivasa, A. Kulachenko / Mechanics of Materials 80 (2015) 13–26
Fig. 6. Different stages of compression of the honeycomb structure.
Table 1 Material parameters for hexagonal honeycomb structure. Young’s modulus [GPa]
Yield stress [MPa]
Tangent modulus [GPa]
Poisson’s ratio [-]
Density [kg/m3]
69.97
292
0.69
0.33
2770
Fig. 7. Mesh convergence.
Table 2 Computational effort comparison for honeycomb model without internal contact. Type of element
Elements per wall length [–]
Elements per wall thickness [–]
Degrees of freedom [–]
Wall time [s]
Cumulative iterations [–]
Beam Solid
20 20
– 2
6671 68303
10,649 17,518
13,702 10,661
18
P. Srinivasa, A. Kulachenko / Mechanics of Materials 80 (2015) 13–26
Table 3 Computational effort comparison for honeycomb model with internal contact. Type of element
Elements per wall length [–]
Elements per wall thickness [–]
Degrees of freedom [–]
Wall time [s]
Cumulative iterations [–]
Beam Solid
20 20
– 2
26,741 68,303
61,817 62,337
13,111 11,358
300
Table 4 Material parameters for the random structure.
Stress [kPa]
250
Young’s modulus [GPa]
Shear modulus [MPa]
Density [kg/m3]
20
769
1600
200 150 100 Beam element Solid element Experimental
50 0
0
10
20
30
40 50 Strain [%]
60
70
80
90
Fig. 8. Elements comparison: beam and solid (plane stress).
solid element. This is because ANSYS incorporates a threedimensional contact search for the beam element, which appeared to be significantly more time-consuming than 2D node-to-surface contact used for beams. We overcome these deficiencies with the beam model by using a custom code (Kulachenko and Uesaka, 2012) for the random structure model of NFC foams, where contact search is limited to two-dimensions and the beams are also two-dimensional. Using the custom code enabled us a possibility to incorporate flexible numerical stabilization strategies and avoid unnecessary overheads. It is also to be noted that the porosity of the honeycomb structure considered here is 96.52% (calculated based on the assumption that the hexagons are regular). The random structure that we will analyse will consider porosities as high as 99.5%. The discretisation required for slender cell walls along its thickness only increases with reduction in thickness. In such a case, plane solid elements will clearly not be viable owing to an increase in computational cost associated with discretising slender cell walls along its thickness to get an acceptable aspect ratio for solid elements.
3. FE model of the random structure The geometry of the random structure is generated in MATLAB as defined in Section 2. The geometry is then transferred to a custom code, which produces the output readable with commercial software. It is a 2 node element with three degrees of freedom (2 translational and 1 rotational) at each node. The beam element is based on engineering theory of bending but also accounts for shear
deformation (Lindström et al., 2013; Przemieniecki, 1968). The longitudinal displacement is assumed to be linearly varying while the transverse displacement is accounted for with cubic Hermite polynomials. The calculations are performed in the local element coordinate system after the finite rotation transformation. A projector-consistent stress-stiffness matrix is used to achieve quadratic rate of convergence (Nour-Omid and Rankin, 1988, 1991). We assume a linear elastic material as the immediate target is to capture geometric non-linearity. The material data is based on those of cellulose and is provided in Table 4. The elastic modulus for that wall material of 20 GPa was chosen based on the considerations that NFC forms sheets that composed the cell wall. The reported elastic modulus for NFC sheets ranges between 10 GPa (Lindström et al., 2013) and 30 GPa (Kulachenko et al., 2012). 3.1. Boundary conditions Compressive loading is one of the most significant loading conditions for foams. The choice of boundary conditions when modelling a representative structure is governed by its ability to yield agreeable average values for macroscopic properties. There have been many studies on the choice of appropriate boundary conditions and it has been shown that representative structures that are periodic are best represented by periodic boundary conditions (Chen et al., 1999; Mesarovic and Padbidri, 2005; Zhu et al., 2000, 2001a,b). We have chosen to model the periodic random representative structure through periodic boundary conditions. In the case of random structure meshed with beam element, the top–bottom and left–right boundaries are meshed equivalently so as to have identical node pairs. For the case of uni-axial compression, the node pairs from the top and bottom are coupled to have identical displacements in the x-direction. The node pairs on the left and right are coupled to have identical displacements in both x and y directions. Similarly, coupling is enforced for rotation with respect to the out of plane axis. Furthermore, one of the nodes is fixed with respect to displacement in x-direction in order to avoid rigid body rotation and translation. Compression is brought about by an applied displacement on the top and bottom boundaries.
P. Srinivasa, A. Kulachenko / Mechanics of Materials 80 (2015) 13–26
3.2. Internal contact modelling The interior of each cell wall is lined with contact and target elements. The contact element used is a 2D node to surface contact element. This element is defined by a single node and is laid on the surface of the beam elements. The target elements are also laid on the interior of the cell walls and all have their normals pointing inwards, that is, towards the centroid. The contact elements do not lock upon contact but is allowed to slide. We have not accounted for friction in this case. The contact algorithm is chosen to be based on penalty formulation. Fig. 9 depicts the contact and target elements (with normals).
19
computed modes are superposed on the initial geometry with a small coefficient. This serves two important purposes. Firstly, it helps overcome the instabilities and convergence issues associated with local buckling of the slender cell walls, especially at high porosities. Secondly, the superposition of these buckling modes on the geometry allows us to introduce some initial curvature to the cell walls. The buckling mode extraction is based on Block Lanczos method. The number of modes chosen to be superposed on the initial geometry is equal to the number of cells in the structure as the computed modes correspond to buckling of individual cell walls. 4. Results
3.3. Size effect 4.1. RVE size To study the size effect, three different sizes, namely, 100, 400 and 900 cells have been chosen. An example of these different sizes can be seen in Fig. 10. In a method proposed by Zhu et al. (2001a,b), and subsequently used by Alsayednoor et al. (2013) in a similar analysis, an irregularity parameter is used to control the initial distribution of seed points for the Voronoi tessellations. The focus of this work is not to study the effects of irregularity and thus, the distribution of seed points is completely random and is Poisson distributed. The different size RVEs are generated by simultaneously increasing the number of seed points and the size of the unit space (refer Section 2.1.1). Moreover, to study the size effect, and arrive at a reasonable representative volume element (RVE) size for studying the effect of porosities, we kept the density unchanged. This implies that, with increase in the number of cells, there is no corresponding decrease in the slenderness of the cell walls. 3.4. Buckling mode superposition Before carrying out the large strain compression, a linear buckling analysis is carried out on the structure. The
Fig. 9. Contact and target elements. (Contact elements are laid over the beam elements on the cell walls. Target elements in the figure are represented by normals pointing towards the centroid of each cell.)
The complete stress–strain curves for three different RVEs with and without contact are presented in Fig. 11. The response of each RVE is computed as a mean of 10 random realisations, and the error bars are symmetric corresponding to two standard deviations in length. The simulations are carried out with the implicit integration method. In the work exploring the RVE size for large strain compressive response of random structures (Alsayednoor et al., 2013) explore the usage of deviation of ‘‘withcontact’’ response from ‘‘without-contact’’, called the contact strain, as a measure of onset of densification. They compare this measure with onset strain of densification as defined by Li et al. (2006) in terms of the maximum of the efficiency function. This is defined as the local maximum of gðÞ, which is given by Eq. (5):
gðÞ ¼
Z 1 rðÞd rðÞ y
ð5Þ
where r is the stress, y is the strain at yield and is the strain. For all the three cases of 100, 400 and 900 cells we observe that the onset strain of densification is much higher than the contact strain. The contact strain decreases quite rapidly with the increase in RVE size. It was shown by Alsayednoor et al. (2013) that the two measures were equivalent and that they converge to the same value as the RVE size is increased. However, this was demonstrated for a given degree of irregularity with an imposed minimum distance between the seed points for the tessellation. In the current case, where there is no restriction imposed on the minimum distance between the seed points in the tessellations and a high degree of irregularity, these two measures do not show any convergence towards each other. However, as can be seen from the plots of with and without contact response, the with-contact model proves to be very important in order to capture the trend of densification at large strains. In Fig. 11, all the responses are superposed on top of each other. We observe that there is a gradual increase in the gradient of plateau stress with increasing RVE sizes. Alsayednoor et al. (2013) have found that increase in irregularity leads to an increase in the positive gradient of the plateau stress. In addition we observe that at very high
20
P. Srinivasa, A. Kulachenko / Mechanics of Materials 80 (2015) 13–26
Fig. 10. An example of the 3 different structures chosen to study the size effects. The linear dimensions specified are in micro metre. (a) 100 cells – 278.01 293.70 (b) 400 cells – 532.80 535.35 (c) 900 cells – 802.86 793.68.
Fig. 11. With and without contact stress–strain curves for RVE size of 100, 400 and 900 cells.
irregularities, an increase in RVE size also results in an increased positive gradient of the plateau stress region. A look at the plot of percentage of area in contact as compared to the total available area may explain the reason for this. In Fig. 12, we plot the percentage of contact area for the three considered RVE sizes. We see that for a given strain, the percentage of area in contact as a function of
total available area increases increased RVE size. Thus the initiation of contact, shown by a black line, and, consequently, its effect on the compressive response begins to appear much earlier as the RVE size is increased. The variability in the stress strain response increases with increase in compressive strain. However, for the 900 cell model, the deviation from mean of the compressive
P. Srinivasa, A. Kulachenko / Mechanics of Materials 80 (2015) 13–26
21
Fig. 12. Percentage of contact area as a function of total available contact area for different RVEs.
stress response is almost negligible up to 70% strains. This suggests that the 900 cell RVE is suitable for the parameter study on effect of porosity. 4.2. Effect of cell wall thickness There is reliable data regarding the thickness of the individual cell walls in the available experimental results. Therefore, before comparing to the experimental data, we need to assess the effects brought by the average wall thickness alone. In the preceding section, as the linear size and the number of cells were increased simultaneously while keeping the relative density constant, the average thickness of the cell walls remained constant across all the three sizes. We now consider the largest RVE in the previous case (900 cells) and allow for variation of thickness while maintaining the relative density constant. Three different structures are considered, each having 900 cells but with average linear size corresponding to 1071 1065 [lm], 802 793 [lm] and 532 533 [lm]. We shall denote them as LS1, LS2 and LS3 respectively. The average thickness of the cell walls now decrease as the linear size is decreased. The results are tabulated below
in Fig. 13. As before, the error bars are symmetric and correspond to one standard deviation on either side of the mean value. The form of stress–strain curves is clearly influenced by thickness. For a given relative density, as the average thickness of the cell walls increase, the slope of the plateau stress region also increases. For the largest linear size with the highest average thickness (LS1), the response after the initial elastic region, is one of monotonic increase with no distinct plateau stress or densification regions. This trend is preserved in both, with contact and without contact responses. This results in an interesting addition to those that were observed by Alsayednoor et al. (2013). They observed that in a model characterised by an increasing degree of irregularity and subsequent decrease of thickness while keeping the relative density constant, led to increase in the size of the positive slope in the plateau region. The results above correspond to a constant and high degree of irregularity. If the degree of irregularity is maintained, and then the thickness is decreased by keeping the number of cells constant but decreasing the linear size, the general trend in behaviour of the plateau region is very clear: an increase
22
P. Srinivasa, A. Kulachenko / Mechanics of Materials 80 (2015) 13–26
Fig. 13. With and without contact stress–strain curves for 900 cell RVE of varying linear sizes.
in thickness increases the slope of the plateau region and decrease the effect brought by inner contact. This is consistent with the experimental data in Sehaqui et al. (2010). A decrease in porosity results not only in stiffer response but also in an increased slope of the stress-plateau region. While there could be more microstructural factors at play, thickness certainly seems to contribute to this trend. 4.3. Influence of the non-straightness of the cell walls on small strain response The non-straightness of the cell wall that we captured by buckling mode superposition, that is carried out before the compression of the random structure, has a marked influence on the small strain response of the structure. We study its influence by controlling the percentage of each buckling mode that is superposed on the original configuration. This is accomplished by first carrying out a linear buckling analysis. The number of modes extracted is equal to the number of cells in the structure. After this, a percentage of the normalised buckling mode shapes are superposed on the original geometry. The results are
plotted in Fig. 14. The number of modes extracted and superposed for each structure is equal to the number of cells. As before, the response of each RVE is computed as a mean of 10 random realisations, and the error bars are symmetric corresponding to two standard deviations in length. For all the 3 RVEs, we observe that there is an initial stiff response corresponding to about 0.01% strain. This initial stiff response is seen to reduce with an increase in the percentage of buckling mode superposition. For a given percentage of buckling mode superposition, stress levels gradually increase with the RVE size. At slightly higher strain of 4%, we observe that the stress levels attained reduces with an increase in degree of non-straightness. In order to quantify this reduction of stiffness with an increase in percentage of buckling mode superposition and RVE size, we measure the secant modulus at 1% and 4% strains. The results are shown in Fig. 15. At 1% strain, the secant modulus decreases drastically with an increase buckling mode superposition. For all the three RVEs that are considered, the secant modulus measured at 1% strain decrease by nearly 40% at 100% buckling mode superposition
P. Srinivasa, A. Kulachenko / Mechanics of Materials 80 (2015) 13–26
23
Fig. 14. Small strain response of RVEs with differing percentage of buckling mode superposition.
Fig. 15. Secant modulus of RVEs measured at 1% and 4% strains.
as compared to 10% buckling mode superposition. The secant modulus computed at 4% strain shows much less variation with an increase in buckling mode superposition.
As compared to 10% buckling mode superposition, the secant modulus at 100% buckling mode superposition reduces by only about 9%.
24
P. Srinivasa, A. Kulachenko / Mechanics of Materials 80 (2015) 13–26
Fig. 16. Comparison with experiments.
The percentage of superposition of the buckling modes appears to have a strong influence at very small strains up to 1%. However at higher strains, the influence of percentage of buckling mode reduces, with the stress levels reaching to 10% of each other irrespective of the percentage of buckling modes that are superposed. This is to say that the form of the stress-strain curve up to 1% strain is strongly influenced by the initial straightness of the cell walls. However, these effects largely even out at strains beyond 4%. The buckling mode superposition proves valuable in achieving convergence with the implicit method in finite element simulation and we demonstrate that it has limited size effects on the nature of stress–strain curve at strains beyond 4%. 4.4. Comparison with experiments for different porosities We compared the results obtained from the 2D model against those from experiments reported in Sehaqui et al. (2010) for 3 different porosities viz. 98.10%, 95.90% and 93.10% in Fig. 16. We assumed the elastic modulus of the cell walls to be equal to 20 GPa, which is probably the upper
estimate based on the tests made on the dense NFC membranes (Kulachenko et al., 2012). The model captures qualitatively the observed trend in experiments but there is only an order of magnitude agreement in the stress values. The gradual stiffening after the initial elastic response observed in the experiments is captured by the twodimensional model with internal contact. As no material non-linearity was introduced in the model of the cell walls, the amount of stiffening captured in the model after the initial elastic response corresponds only to geometrical non-linearity. In a relative scale it captures the trend observed in the experiments rather well. There could be several reasons behind the large difference in the actual stress values between the model and simulation. First is the error in estimating the wall thickness or the average cell size of the cells at a given porosity. We demonstrated, however, that although these parameters change the stiffening part of the response, they cannot explain the order of observed differences. The second possibility is underestimating the elastic properties of the cellulose wall material. We used however the upper estimate for the isotropic nanocellulose network and the elastic
P. Srinivasa, A. Kulachenko / Mechanics of Materials 80 (2015) 13–26
properties of such networks are not size dependent for the sizes larger than a micron (Kulachenko et al., 2012). The introduction of the non-linear constitutive behaviour and damage can only decrease the calculated stress value and move it even further away from the experimental values. The only two remaining possibilities are (1) inadequate description of the geometry and (2) the inability of the considered 2D foam representation to capture the response of the 3D structure. The geometry of the NFC foam inner structure can be affected by the dehydration process or by gravity which can both contribute to the creation of the preferred orientation of the struts. Since the considered foams were traditionally tested in the direction normal to the free surface (owing to the difficulties in sample preparation), it is impossible to rule out such a possibility for the available experimental data. The 2D foam representation can potentially introduce an error in capturing the response however it is unlikely that it will make it stiffer as required to match the experimental results since open cell foams gives the upper estimate of the stiffness and crushing strength at a given density (Gibson and Ashby, 1982). This means that anisotropy in inner structure is most plausible explanation of the observed difference. It can be tested by ether reconstructing the inner structure with, for example, micro-CT equipment or by simply performing the tests in various direction of the sample. 5. Summary and conclusions We studied the compressive response of irregular foam structures at high strains using 2D beam representation with internal contact. We have proved the effectiveness of beam element as opposed to solid element in modelling high porosity structures by reproducing the experimental and numerical results of a previously published work, using both beam and solid elements and concluding that beam elements allow for capturing the local bucking of the cell walls while also being computationally more efficient. The size effect for a given relative density was explored by increasing the linear size and the number of cells simultaneously. After finding the size that provides statistically acceptable results up to nearly 70% compressive strains, we explored the effect of thickness at a given porosity and the effect of porosity itself and related it to the experimental data. For every numerical experiment, we single out the effect of inner contact by performing the simulations with and without contact. Internal contact is shown to be of significance in simulating the densification regime of random structures and is the only way to capture the stiffening response at the compressed state. The cell walls thickness has a significant influence on the general trend in the stress-plateau region. An increase in thickness is seen to increase the slope of the stress-plateau and diminishes the difference brought by the inner contact as the total free surface reduces for the given porosity. The non-straightness of the cell walls significantly affects the initial elastic response, smoothening the
25
transition between the initial elastic region and the plateau. The non-straightness has however a negligible influence on the stress levels reached in the plateau region and subsequent stiffening. The random periodic two dimensional models also capture the trend of the compressive stress–strain response seen in experimental results from NFC foams. While there is an order magnitude agreement in the simulated response of this random model in comparison with the experiments of NFC foams, it is seen that neither the effect of stiffness nor the size dependency can explain the difference between the experimental data and simulation results. Although the considered NFC foam predominantly consists of open cell structures, it has non-trivial out of plane connectivity which requires 3D representation. Acknowledgements Financial support from The Swedish Research Council is gratefully acknowledged by the authors. P. Srinivasa wishes to acknowledge Dr. Houssine Sehaqui for providing the data from experiments on compression of foams. References Ali, Z.M., Gibson, L.J., 2013. The structure and mechanics of nanofibrillar cellulose foams. Soft Matter 9, 1580. http://dx.doi.org/10.1039/ c2sm27197d. Alsayednoor, J., Harrison, P., Guo, Z., 2013. Large strain compressive response of 2-D periodic representative volume element for random foam microstructures. Mech. Mater. 66, 7–20. http://dx.doi.org/ 10.1016/j.mechmat.2013.06.006. Aulin, C., Netrval, J., Wågberg, L., Lindström, T., 2010. Aerogels from nanofibrillated cellulose with tunable oleophobicity. Soft Matter 6, 3298. http://dx.doi.org/10.1039/c001939a. Chen, C., Lu, T.J., Fleck, N.A., 1999. Effect of imperfections on the yielding of two-dimensional foams. J. Mater. Chem. 47, 2235–2272. Gibson, L.J., Ashby, M.F., 1982. The mechanics of three-dimensional cellular materials. Proc. R. Soc. A Math. Phys. Eng. Sci. 382, 43–59. http://dx.doi.org/10.1098/rspa.1982.0088. Gibson, L.J., Ashby, M.F., 1988. Cellular Solids. Pergamon press, New York. Gong, L., Kyriakides, S., Jang, W.-Y., 2005. Compressive response of opencell foams. Part I: Morphology and elastic properties. Int. J. Solids Struct. 42, 1355–1379. http://dx.doi.org/10.1016/ j.ijsolstr.2004.07.023. Klemm, D., Kramer, F., Moritz, S., Lindström, T., Ankerfors, M., Gray, D., Dorris, A., 2011. Nanocelluloses: a new family of nature-based materials. Angew. Chem. Int. Ed. Engl. 50, 5438–5466. http:// dx.doi.org/10.1002/anie.201001273. Kulachenko, A., Uesaka, T., 2012. Direct simulations of fiber network deformation and failure. Mech. Mater. 51, 1–14. http://dx.doi.org/ 10.1016/j.mechmat.2012.03.010. Kulachenko, A., Denoyelle, T., Galland, S., Lindström, S.B., 2012. Elastic properties of cellulose nanopaper. Cellulose 19, 793–807. http:// dx.doi.org/10.1007/s10570-012-9685-5. Li, Q.M., Magkiriadis, I., Harrigan, J.J., 2006. Compressive strain at the onset of densification of cellular solids. J. Cell. Plast. 42, 371–392. http://dx.doi.org/10.1177/0021955X06063519. Lindström, S.B., Kulachenko, A., Jawerth, L.M., Vader, D.A., 2013. Finitestrain, finite-size mechanics of rigidly cross-linked biopolymer networks. Soft Matter 9, 7302. http://dx.doi.org/10.1039/ c3sm50451d. Mesarovic, S.D., Padbidri, J., 2005. Minimal kinematic boundary conditions for simulations of disordered microstructures. Philos. Mag. 85, 65–78. http://dx.doi.org/10.1080/14786430412331313321. Nour-Omid, B., Rankin, C.C., 1988. The use of projectors to improve finite element performance. Comput. Struct. 30, 257–267. Nour-Omid, B., Rankin, C.C., 1991. Finite rotation analysis and consistent linearization using projectors. Comput. Methods Appl. Mech. Eng. 93, 353–384. http://dx.doi.org/10.1016/0045-7825(91)90248-5. Okabe, A., Boots, B., Sugihara, K., 1995. Spatial Tessellations: Concepts and Applications of Voronoi Diagrams, first ed. John Wiley & Sons Inc., West Sussex, England.
26
P. Srinivasa, A. Kulachenko / Mechanics of Materials 80 (2015) 13–26
Papka, S.D., Kyriakides, S., 1994. In-plane compressive response and crushing of honeycomb. J. Mech. Phys. Solids 42, 1499–1532. Papka, S.D., Kyriakides, S., 1998. Experiments and full-scale numerical simulations of in-plane crushing of a honeycomb. Acta Mater. 46, 2765–2776. http://dx.doi.org/10.1016/S1359-6454(97)00453-9. Przemieniecki, J.S., 1968. Theory of Matrix Structural Analysis, first ed. McGraw-Hill Book Co., Inc., New York. Roberts, A.P., Garboczi, E.J., 2002. Elastic properties of model random three-dimensional open-cell solids. J. Mech. Phys. Solids 50, 33–55. http://dx.doi.org/10.1016/S0022-5096(01)00056-4. Sehaqui, H., 2011. Nanofiber networks, aerogels and biocomposites based on nanofibrillated cellulose from wood. Doctoral thesis, KTH – Royal Institute of Technology. Sehaqui, H., Salajková, M., Zhou, Q., Berglund, L.A., 2010. Mechanical performance tailoring of tough ultra-high porosity foams prepared from cellulose I nanofiber suspensions. Soft Matter 6, 1824. http:// dx.doi.org/10.1039/b927505c. Silva, M.J., Hayes, W.C., Gibson, L.J., 1995. The effects of non-periodic microstructure on the elastic properties of two-dimensional cellular solids. Int. J. Mech. Sci. 37, 1161–1177. Svagan, A.J., Samir, M.A.S.A., Berglund, L.A., 2008. Biomimetic foams of high mechanical performance based on nanostructured cell walls
reinforced by native cellulose nanofibrils. Adv. Mater. 20, 1263–1269. http://dx.doi.org/10.1002/adma.200701215. Svagan, A.J., Jensen, P., Dvinskikh, S.V., Furó, I., Berglund, L.A., 2010. Towards tailored hierarchical structures in cellulose nanocomposite biofoams prepared by freezing/freeze-drying. J. Mater. Chem. 20, 6646. http://dx.doi.org/10.1039/c0jm00779j. Vanderburg, M.W.D., Shulmeister, V., Van Der Giessen, E., Marrisen, R., 1997. On the linear elastic properties of regular and random open-cell foam models. J. Cell. Plast. 33, 31–54. http://dx.doi.org/10.1177/ 0021955X9703300103. Zhu, H.X., Windle, A.H., 2002. Effects of cell irregularity on the high strain compression of open-cell foams. Acta Mater. 50, 1041–1052. http:// dx.doi.org/10.1016/S1359-6454(01)00402-5. Zhu, H.X., Hobdell, J.R., Windle, A.H., 2000. Effects of cell irregularity on the elastic properties of open-cell foams. Acta Mater. 48, 4893–4900. http://dx.doi.org/10.1016/S1359-6454(00)00282-2. Zhu, H.X., Hobdell, J.R., Windle, A.H., 2001a. Effects of cell irregularity on the elastic properties of 2D Voronoi honeycombs. J. Mech. Phys. Solids 49, 857–870. Zhu, H.X., Thorpe, S.M., Windle, A.H., 2001b. The geometrical properties of irregular two-dimensional Voronoi tessellations. Philos. Mag. 81, 2765–2783.
Contents lists available at ScienceDirect
Mechanics of Materials journal homepage: www.elsevier.com/locate/mechmat
Analysis of the compressive response of Nano Fibrillar Cellulose foams Prashanth Srinivasa 1, Artem Kulachenko ⇑ KTH Royal Institute of Technology, School of Engineering Sciences, Department of Solid Mechanics, SE-100 44 Stockholm, Sweden
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 8 June 2014 Received in revised form 1 September 2014 Available online 5 October 2014
Nano Fibrillar Cellulose (NFC) is fast emerging as a biomaterial with promising applications, one of which is cellular foam. The inner structure of the foam can take various shapes and hierarchical micro-structures depending on the manufacturing parameters. The compressive response of foams developed from these materials is currently a primary criterion for the material development. In this work, we focus on the connection between the non-linear part of the response and the inner structure of the material. We study the effect of internal contact and its contribution to gradual stiffening in the energy absorbing region and accelerated stiffening in the densification region of the large strain compressive response. We use the finite element method in this study and discuss the applicability and efficiency of different modelling techniques by considering well defined geometries and available experimental data. The relative contribution of internal contact is singled out and mapped onto the overall compressive response of the material. The effect of initial non-straightness of the cell walls is studied through superposing differing percentages of the buckling modes on the initial geometry. The initial non-straightness is seen to have a significant effect for only strains up to 1%. The secant modulus measured at slightly higher strains of 4%, demonstrates lesser effect from the non-straightness of cell walls. The simulations capture the compressive response well into the densification regime and there is an order of magnitude agreement in between simulations and experiments. We observed that internal contact is crucial for capturing the trend of compressive response. Ó 2014 Elsevier Ltd. All rights reserved.
Keywords: NFC foams Internal contact Voronoi structures Foam densification Effect of foam porosity
1. Introduction Low density cellular materials like foams have found applications in a wide range of fields. They are used as packaging, shock-absorbing and insulation material. Polymeric foams are leading the way in such applications. Concurrently, a new nature-based material, cellulose nanofibrils (Nano Fibrillar Cellulose, NFC), is fast emerging as an
⇑ Corresponding author. Tel.: +46 87908944. E-mail addresses: [email protected] (A. Kulachenko). 1 Tel.: +46 707948887.
(P.
Srinivasa),
http://dx.doi.org/10.1016/j.mechmat.2014.09.006 0167-6636/Ó 2014 Elsevier Ltd. All rights reserved.
[email protected]
attractive alternative to synthetic polymers with potential applications in biomedical, structural, chemical and numerous other areas (Klemm et al., 2011). The main advantage with using cellulose fibrils is their high stiffness-to-density ratio, wide availability and biocompatibility. Nanofibrils have already demonstrated their potential in reinforcing the cell walls of polymer foams (Svagan et al., 2008), in the manufacturing of high porosity foams from nanofibre suspensions (Sehaqui et al., 2010) and in aerogels which exhibit foam like structure (Aulin et al., 2010). Depending on the inner hierarchical structure and porosity of these materials, the mechanical properties
14
P. Srinivasa, A. Kulachenko / Mechanics of Materials 80 (2015) 13–26
can vary in a wide range. Naturally, providing insights on how the microstructure of foams affects the mechanical response is essential for material engineers as it provides opportunity to establish the mechanical property space of these newly emerging materials. The task of connecting the micro-structure of the material with mechanical response has been a topic of intense research. The mechanical properties of polymeric foam structures and other cellular solids were initially studied through the analysis of a single cell, e.g. hexagonal honeycombs (Gibson and Ashby, 1988). Such analyses proved to be very successful in predicting the elastic properties, but due to a complex nature of the compressive response associated with local cell wall buckling and inner contact, the non-linear response has been difficult to address. Furthermore, there is noticeable difference in response of the nanocellulose (NFC) foam structure under compressive loading as compared to that of the polymeric foams with more ordered cellular structure. The stress plateau or the energy absorbing zone of the compressive response in NFC foams sees gradual stiffening of the response (Fig. 1) as opposed to open cell polymeric foams which exhibit a clear plateau (Gong et al., 2005). This stiffening has been attributed to the better load bearing capability of the NFC cell walls (Sehaqui, 2011). Thus it is imperative that there are microstructural factors at play which should account for the distinct response of the NFC foams under compression. In this work, we will provide insights to the micromechanics of non-linear compressive behaviour of NFC foams. We will begin by modelling a two dimensional foam structure with random Voronoi tessellations. The aim is to establish a numerical model that will be able to capture the stiffening response in the energy absorbing zone as well as the densification at large strains. The suitability of beam element as compared to the plane strain solid element will be analysed for simple hexagonal honeycomb geometry. Finally, the effect of internal contact on the stress response is studied and its contribution to
the gradual stiffening in the energy absorbing and the stiffening region will be analysed. 2. Methods NFC foams are manufactured through the process of freeze drying. When the NFC suspension is freeze dried, nucleation process is initiated for the formation of ice crystals (Svagan et al., 2010). The geometry of the foam like structures is well captured by the Voronoi tessellations. A number of researchers (Roberts and Garboczi, 2002; Silva et al., 1995; Vanderburg et al., 1997; Zhu and Windle, 2002) have studied the Voronoi structure in both two- and three-dimensions. Zhu and Windle (2002) have also studied the effects of cell irregularity on the high strain compression. The stress–strain responses of the random foam structures in the studies above show a clear transition between the linear elastic regime and stress plateau. Furthermore, the stress plateau is truly a plateau with marginal increase in compressive stress levels before densification. However, as we indicated earlier, the experimental results of the compressive stress–strain response of NFC foams suggest that there is a gradual stiffening even in the so called ‘‘stress plateau’’ (Ali and Gibson, 2013; Sehaqui et al., 2010) and the transition between the linear response and stress plateau is not distinct. Studies based on SEM images of NFC foams suggest that the structure is made up of irregular sheets with occasional connectivity to other sheets (Ali and Gibson, 2013). The compressive response is reported to fit well with those predicted by theoretical estimates on open cell foams. The highly irregular nature of the foam structure makes it a good candidate for being modelled through Voronoi tessellations. The connectivity in the out-of-plane direction is non-trivial and might be a limiting factor for a 2D model. However, as the results from simulations will testify, it is an agreeable start. 2.1. Modelling of the Voronoi structure There are two aspects to the modelling of NFC foams as Voronoi tessellations. The first one corresponds to the generation of a periodic structure that lends itself to the application of periodic boundary conditions and secondly, the generation of finite random thickness cell walls for the periodic structure.
8 99.5% 98.1% 95.9% 93.1%
7
Stress [MPa]
6 5 4 3 2 1 0 0
10
20
30
40 50 Strain [%]
60
70
80
Fig. 1. Compressive response of NFC foams of different porosities (Plotted from experimental data provided by Houssine Sehaqui (Sehaqui, 2011)).
2.1.1. Periodic microstructure The generation of periodic structure rests on the property of the tessellations. We consider a set of random points inside a unit square in two-dimensional Euclidean space. Additional points are created around this unit square which are translations of the random points inside the unit square (Okabe et al., 1995). The translations of these points are governed by the transformations in Eqs. (1) and (2).
xo ¼ xi pL
ð1Þ
yo ¼ yi qL
ð2Þ
P. Srinivasa, A. Kulachenko / Mechanics of Materials 80 (2015) 13–26
15
Here, xo , yo are the nucleating points outside of the unit cell of dimensions L L; xi and yi are points within the unit cell and p; q are integers. This is demonstrated in Fig. 2 where the points inside the central square are translated all around the square unit cell. The Voronoi tessellation, of the resulting space which now has an area 9 times that of the unit cell, is shown in Fig. 3. It is seen that the central square is now periodic and can tile the whole two dimensional (2D) space. The tessellation consists of polygons, all of which are convex and have random number of edges. 2.2. Finite thickness cell walls The aspect of having finite random thickness cells is accomplished by dilating these polygons towards their centroidal points by using Eqs. (3) and (4):
xnew ¼ A x þ ð1 AÞ xc
ð3Þ
ynew ¼ A y þ ð1 AÞ yc
ð4Þ Fig. 3. Periodic Voronoi tessellations.
Here xnew and ynew corresponds to the vertices of the deflated polygon, x and y corresponds to the vertices of the initial polygon, xc and yc corresponds to the centroid of the initial polygon, and A corresponds to the factor of dilation. As can be seen, this result in variable thickness for each of the edges of the polygon as the deflation is not uniform. This is because each of the vertices which are at different distance from the centroid is scaled by one single factor. A binary subtraction of the dilated polygons from their respective parent polygons leaves us behind with the necessary finite cell wall thickness, random porous structure as shown in Fig. 4 The variation in porosity can be represented either by varying the thickness of cell walls or by having lesser number of cells in a given space. In this study, we will limit ourselves by varying porosity through altering the cell wall
Fig. 4. Periodic unit cell with finite thickness cell walls.
thickness. The thickness of the cell walls is varied through the deflation factor in order to obtain structures with different porosities. 2.3. Appropriate choice of finite element
Fig. 2. Translation of the tessellation points.
The model of random structures with high porosities consists of slender cell walls which undergo flexural deformation under compressive loads. It is essential to make a reasoned choice on the finite element used to discretise the structure. We intend to compare the effectiveness of 2D plane solid elements against beam elements. Plane elements are occasionally used for similar analyses as they can offer a better representation of the cell wall structure, particularly at the vertexes. We do the comparison through simulating the compression of pre-stressed aluminium foam with hexagonal honeycomb structure, for which we found well defined material data and comprehensive
16
P. Srinivasa, A. Kulachenko / Mechanics of Materials 80 (2015) 13–26
experimental results. The thickness of the aluminium honeycomb cell walls are comparable to that of the random structure that we wish to analyse, making it a good choice for the purpose of choosing a suitable finite element. Both experimental and numerical results have been previously reported in two successive papers (Papka and Kyriakides, 1994, 1998). 2.4. Hexagonal honeycomb structure The geometry consists of aluminium sheets which are bonded along its length at regular intervals as shown in Fig. 5(a). The sheets are then stretched by applying pressure resulting in the formation of hexagonal honeycomb Fig. 5(c). This process of stretching is meant to simulate the actual manufacturing process and contributes to the pre-stressed state in the cell walls. The applied pressure is then released to relax the structure to the new equilibrium state (Fig. 5(d)). Following the relaxation step, large strain compression is carried out through applied displacement on upper and bottom sides of the structure (Fig. 6). The numerical simulations are carried out using a commercial, implicit time integration finite element code (ANSYS). Two different elements are used for the purpose of comparison. First a three node finite strain beam element (BEAM189 in ANSYS) with quadratic displacement behaviour is considered. Each node has three degrees of freedom in translation and three degrees of freedom in rotation. The second element is a quadrilateral 2D plane strain solid element (PLANE183 in ANSYS) with 8 nodes including mid nodes. Each node has two degrees of freedom in translation and incorporates quadratic displacement behaviour. A bilinear isotropic hardening material model is used and the material data is presented in Table 1. The discretisation along the length of cell walls is gradually increased and it is found that convergence is achieved for 20 elements per length of the cell wall. The obtained response is stiffer and the point of first yield was higher than the experimental results. This is attributed to the shear locking phenomenon. Bending energy is proportional to the cube of thickness whereas shear energy is proportional to the first power of thickness. As the thickness
becomes increasingly small, the contribution from bending energy reduces more quickly than from shear energy. This of course can be overcome with mesh refinement and changing the aspect ratio of elements. The cell wall is now discretised along its thickness and it is found that a minimum of 2 elements per thickness is required to achieve mesh convergence. The results are shown in Fig. 7. Finally, a comparison is carried out between beam and solid elements with identical discretisation of twenty elements per length of the cell wall. The simulations were carried out on an Intel (R) Xeon (R) E7-8837 processor (Frequency: 2.66 GHz) and 8 cores were used. The data for computational effort required without and with internal contact is presented in Tables 2 and 3 respectively. The results from both the solid and beam elements along with the experimental results are presented in Fig. 8. The initial linear elastic regime and the final stiffening due to densification are represented equally well by both the beam and solid elements. In the stress plateau region, the beam elements capture localised buckling of the cell walls while the solid element predicts a flatter response. The overall agreement with the experimental results is exceptionally good. The computational effort in terms of the wall time is two times lower for the beam element as compared to the solid element when internal contact is not incorporated. This is surprisingly small difference taking into account the fact that the beam model has almost ten times smaller number of active degrees of freedom. The reason for that is we used the 3D beam elements with constrained out-of-plane degrees of freedom. (The reason 2D beam elements were not used was unavailability of the non-linear stabilization tools for these elements in ANSYS, which were critical for avoiding convergence problems passing over the buckling bifurcation points.) Furthermore, the total number of integration points was actually greater in the model with beam element since the selected beam element uses 5 integration points through thickness as opposed to 4 in the case of two plane elements through the thickness. When the internal contact is considered, the computation time for beam elements is comparable to that of the
Fig. 5. Creation of hexagonal honeycomb structure: (a), (b), (c) – expansion (d) – relaxation.
17
P. Srinivasa, A. Kulachenko / Mechanics of Materials 80 (2015) 13–26
Fig. 6. Different stages of compression of the honeycomb structure.
Table 1 Material parameters for hexagonal honeycomb structure. Young’s modulus [GPa]
Yield stress [MPa]
Tangent modulus [GPa]
Poisson’s ratio [-]
Density [kg/m3]
69.97
292
0.69
0.33
2770
Fig. 7. Mesh convergence.
Table 2 Computational effort comparison for honeycomb model without internal contact. Type of element
Elements per wall length [–]
Elements per wall thickness [–]
Degrees of freedom [–]
Wall time [s]
Cumulative iterations [–]
Beam Solid
20 20
– 2
6671 68303
10,649 17,518
13,702 10,661
18
P. Srinivasa, A. Kulachenko / Mechanics of Materials 80 (2015) 13–26
Table 3 Computational effort comparison for honeycomb model with internal contact. Type of element
Elements per wall length [–]
Elements per wall thickness [–]
Degrees of freedom [–]
Wall time [s]
Cumulative iterations [–]
Beam Solid
20 20
– 2
26,741 68,303
61,817 62,337
13,111 11,358
300
Table 4 Material parameters for the random structure.
Stress [kPa]
250
Young’s modulus [GPa]
Shear modulus [MPa]
Density [kg/m3]
20
769
1600
200 150 100 Beam element Solid element Experimental
50 0
0
10
20
30
40 50 Strain [%]
60
70
80
90
Fig. 8. Elements comparison: beam and solid (plane stress).
solid element. This is because ANSYS incorporates a threedimensional contact search for the beam element, which appeared to be significantly more time-consuming than 2D node-to-surface contact used for beams. We overcome these deficiencies with the beam model by using a custom code (Kulachenko and Uesaka, 2012) for the random structure model of NFC foams, where contact search is limited to two-dimensions and the beams are also two-dimensional. Using the custom code enabled us a possibility to incorporate flexible numerical stabilization strategies and avoid unnecessary overheads. It is also to be noted that the porosity of the honeycomb structure considered here is 96.52% (calculated based on the assumption that the hexagons are regular). The random structure that we will analyse will consider porosities as high as 99.5%. The discretisation required for slender cell walls along its thickness only increases with reduction in thickness. In such a case, plane solid elements will clearly not be viable owing to an increase in computational cost associated with discretising slender cell walls along its thickness to get an acceptable aspect ratio for solid elements.
3. FE model of the random structure The geometry of the random structure is generated in MATLAB as defined in Section 2. The geometry is then transferred to a custom code, which produces the output readable with commercial software. It is a 2 node element with three degrees of freedom (2 translational and 1 rotational) at each node. The beam element is based on engineering theory of bending but also accounts for shear
deformation (Lindström et al., 2013; Przemieniecki, 1968). The longitudinal displacement is assumed to be linearly varying while the transverse displacement is accounted for with cubic Hermite polynomials. The calculations are performed in the local element coordinate system after the finite rotation transformation. A projector-consistent stress-stiffness matrix is used to achieve quadratic rate of convergence (Nour-Omid and Rankin, 1988, 1991). We assume a linear elastic material as the immediate target is to capture geometric non-linearity. The material data is based on those of cellulose and is provided in Table 4. The elastic modulus for that wall material of 20 GPa was chosen based on the considerations that NFC forms sheets that composed the cell wall. The reported elastic modulus for NFC sheets ranges between 10 GPa (Lindström et al., 2013) and 30 GPa (Kulachenko et al., 2012). 3.1. Boundary conditions Compressive loading is one of the most significant loading conditions for foams. The choice of boundary conditions when modelling a representative structure is governed by its ability to yield agreeable average values for macroscopic properties. There have been many studies on the choice of appropriate boundary conditions and it has been shown that representative structures that are periodic are best represented by periodic boundary conditions (Chen et al., 1999; Mesarovic and Padbidri, 2005; Zhu et al., 2000, 2001a,b). We have chosen to model the periodic random representative structure through periodic boundary conditions. In the case of random structure meshed with beam element, the top–bottom and left–right boundaries are meshed equivalently so as to have identical node pairs. For the case of uni-axial compression, the node pairs from the top and bottom are coupled to have identical displacements in the x-direction. The node pairs on the left and right are coupled to have identical displacements in both x and y directions. Similarly, coupling is enforced for rotation with respect to the out of plane axis. Furthermore, one of the nodes is fixed with respect to displacement in x-direction in order to avoid rigid body rotation and translation. Compression is brought about by an applied displacement on the top and bottom boundaries.
P. Srinivasa, A. Kulachenko / Mechanics of Materials 80 (2015) 13–26
3.2. Internal contact modelling The interior of each cell wall is lined with contact and target elements. The contact element used is a 2D node to surface contact element. This element is defined by a single node and is laid on the surface of the beam elements. The target elements are also laid on the interior of the cell walls and all have their normals pointing inwards, that is, towards the centroid. The contact elements do not lock upon contact but is allowed to slide. We have not accounted for friction in this case. The contact algorithm is chosen to be based on penalty formulation. Fig. 9 depicts the contact and target elements (with normals).
19
computed modes are superposed on the initial geometry with a small coefficient. This serves two important purposes. Firstly, it helps overcome the instabilities and convergence issues associated with local buckling of the slender cell walls, especially at high porosities. Secondly, the superposition of these buckling modes on the geometry allows us to introduce some initial curvature to the cell walls. The buckling mode extraction is based on Block Lanczos method. The number of modes chosen to be superposed on the initial geometry is equal to the number of cells in the structure as the computed modes correspond to buckling of individual cell walls. 4. Results
3.3. Size effect 4.1. RVE size To study the size effect, three different sizes, namely, 100, 400 and 900 cells have been chosen. An example of these different sizes can be seen in Fig. 10. In a method proposed by Zhu et al. (2001a,b), and subsequently used by Alsayednoor et al. (2013) in a similar analysis, an irregularity parameter is used to control the initial distribution of seed points for the Voronoi tessellations. The focus of this work is not to study the effects of irregularity and thus, the distribution of seed points is completely random and is Poisson distributed. The different size RVEs are generated by simultaneously increasing the number of seed points and the size of the unit space (refer Section 2.1.1). Moreover, to study the size effect, and arrive at a reasonable representative volume element (RVE) size for studying the effect of porosities, we kept the density unchanged. This implies that, with increase in the number of cells, there is no corresponding decrease in the slenderness of the cell walls. 3.4. Buckling mode superposition Before carrying out the large strain compression, a linear buckling analysis is carried out on the structure. The
Fig. 9. Contact and target elements. (Contact elements are laid over the beam elements on the cell walls. Target elements in the figure are represented by normals pointing towards the centroid of each cell.)
The complete stress–strain curves for three different RVEs with and without contact are presented in Fig. 11. The response of each RVE is computed as a mean of 10 random realisations, and the error bars are symmetric corresponding to two standard deviations in length. The simulations are carried out with the implicit integration method. In the work exploring the RVE size for large strain compressive response of random structures (Alsayednoor et al., 2013) explore the usage of deviation of ‘‘withcontact’’ response from ‘‘without-contact’’, called the contact strain, as a measure of onset of densification. They compare this measure with onset strain of densification as defined by Li et al. (2006) in terms of the maximum of the efficiency function. This is defined as the local maximum of gðÞ, which is given by Eq. (5):
gðÞ ¼
Z 1 rðÞd rðÞ y
ð5Þ
where r is the stress, y is the strain at yield and is the strain. For all the three cases of 100, 400 and 900 cells we observe that the onset strain of densification is much higher than the contact strain. The contact strain decreases quite rapidly with the increase in RVE size. It was shown by Alsayednoor et al. (2013) that the two measures were equivalent and that they converge to the same value as the RVE size is increased. However, this was demonstrated for a given degree of irregularity with an imposed minimum distance between the seed points for the tessellation. In the current case, where there is no restriction imposed on the minimum distance between the seed points in the tessellations and a high degree of irregularity, these two measures do not show any convergence towards each other. However, as can be seen from the plots of with and without contact response, the with-contact model proves to be very important in order to capture the trend of densification at large strains. In Fig. 11, all the responses are superposed on top of each other. We observe that there is a gradual increase in the gradient of plateau stress with increasing RVE sizes. Alsayednoor et al. (2013) have found that increase in irregularity leads to an increase in the positive gradient of the plateau stress. In addition we observe that at very high
20
P. Srinivasa, A. Kulachenko / Mechanics of Materials 80 (2015) 13–26
Fig. 10. An example of the 3 different structures chosen to study the size effects. The linear dimensions specified are in micro metre. (a) 100 cells – 278.01 293.70 (b) 400 cells – 532.80 535.35 (c) 900 cells – 802.86 793.68.
Fig. 11. With and without contact stress–strain curves for RVE size of 100, 400 and 900 cells.
irregularities, an increase in RVE size also results in an increased positive gradient of the plateau stress region. A look at the plot of percentage of area in contact as compared to the total available area may explain the reason for this. In Fig. 12, we plot the percentage of contact area for the three considered RVE sizes. We see that for a given strain, the percentage of area in contact as a function of
total available area increases increased RVE size. Thus the initiation of contact, shown by a black line, and, consequently, its effect on the compressive response begins to appear much earlier as the RVE size is increased. The variability in the stress strain response increases with increase in compressive strain. However, for the 900 cell model, the deviation from mean of the compressive
P. Srinivasa, A. Kulachenko / Mechanics of Materials 80 (2015) 13–26
21
Fig. 12. Percentage of contact area as a function of total available contact area for different RVEs.
stress response is almost negligible up to 70% strains. This suggests that the 900 cell RVE is suitable for the parameter study on effect of porosity. 4.2. Effect of cell wall thickness There is reliable data regarding the thickness of the individual cell walls in the available experimental results. Therefore, before comparing to the experimental data, we need to assess the effects brought by the average wall thickness alone. In the preceding section, as the linear size and the number of cells were increased simultaneously while keeping the relative density constant, the average thickness of the cell walls remained constant across all the three sizes. We now consider the largest RVE in the previous case (900 cells) and allow for variation of thickness while maintaining the relative density constant. Three different structures are considered, each having 900 cells but with average linear size corresponding to 1071 1065 [lm], 802 793 [lm] and 532 533 [lm]. We shall denote them as LS1, LS2 and LS3 respectively. The average thickness of the cell walls now decrease as the linear size is decreased. The results are tabulated below
in Fig. 13. As before, the error bars are symmetric and correspond to one standard deviation on either side of the mean value. The form of stress–strain curves is clearly influenced by thickness. For a given relative density, as the average thickness of the cell walls increase, the slope of the plateau stress region also increases. For the largest linear size with the highest average thickness (LS1), the response after the initial elastic region, is one of monotonic increase with no distinct plateau stress or densification regions. This trend is preserved in both, with contact and without contact responses. This results in an interesting addition to those that were observed by Alsayednoor et al. (2013). They observed that in a model characterised by an increasing degree of irregularity and subsequent decrease of thickness while keeping the relative density constant, led to increase in the size of the positive slope in the plateau region. The results above correspond to a constant and high degree of irregularity. If the degree of irregularity is maintained, and then the thickness is decreased by keeping the number of cells constant but decreasing the linear size, the general trend in behaviour of the plateau region is very clear: an increase
22
P. Srinivasa, A. Kulachenko / Mechanics of Materials 80 (2015) 13–26
Fig. 13. With and without contact stress–strain curves for 900 cell RVE of varying linear sizes.
in thickness increases the slope of the plateau region and decrease the effect brought by inner contact. This is consistent with the experimental data in Sehaqui et al. (2010). A decrease in porosity results not only in stiffer response but also in an increased slope of the stress-plateau region. While there could be more microstructural factors at play, thickness certainly seems to contribute to this trend. 4.3. Influence of the non-straightness of the cell walls on small strain response The non-straightness of the cell wall that we captured by buckling mode superposition, that is carried out before the compression of the random structure, has a marked influence on the small strain response of the structure. We study its influence by controlling the percentage of each buckling mode that is superposed on the original configuration. This is accomplished by first carrying out a linear buckling analysis. The number of modes extracted is equal to the number of cells in the structure. After this, a percentage of the normalised buckling mode shapes are superposed on the original geometry. The results are
plotted in Fig. 14. The number of modes extracted and superposed for each structure is equal to the number of cells. As before, the response of each RVE is computed as a mean of 10 random realisations, and the error bars are symmetric corresponding to two standard deviations in length. For all the 3 RVEs, we observe that there is an initial stiff response corresponding to about 0.01% strain. This initial stiff response is seen to reduce with an increase in the percentage of buckling mode superposition. For a given percentage of buckling mode superposition, stress levels gradually increase with the RVE size. At slightly higher strain of 4%, we observe that the stress levels attained reduces with an increase in degree of non-straightness. In order to quantify this reduction of stiffness with an increase in percentage of buckling mode superposition and RVE size, we measure the secant modulus at 1% and 4% strains. The results are shown in Fig. 15. At 1% strain, the secant modulus decreases drastically with an increase buckling mode superposition. For all the three RVEs that are considered, the secant modulus measured at 1% strain decrease by nearly 40% at 100% buckling mode superposition
P. Srinivasa, A. Kulachenko / Mechanics of Materials 80 (2015) 13–26
23
Fig. 14. Small strain response of RVEs with differing percentage of buckling mode superposition.
Fig. 15. Secant modulus of RVEs measured at 1% and 4% strains.
as compared to 10% buckling mode superposition. The secant modulus computed at 4% strain shows much less variation with an increase in buckling mode superposition.
As compared to 10% buckling mode superposition, the secant modulus at 100% buckling mode superposition reduces by only about 9%.
24
P. Srinivasa, A. Kulachenko / Mechanics of Materials 80 (2015) 13–26
Fig. 16. Comparison with experiments.
The percentage of superposition of the buckling modes appears to have a strong influence at very small strains up to 1%. However at higher strains, the influence of percentage of buckling mode reduces, with the stress levels reaching to 10% of each other irrespective of the percentage of buckling modes that are superposed. This is to say that the form of the stress-strain curve up to 1% strain is strongly influenced by the initial straightness of the cell walls. However, these effects largely even out at strains beyond 4%. The buckling mode superposition proves valuable in achieving convergence with the implicit method in finite element simulation and we demonstrate that it has limited size effects on the nature of stress–strain curve at strains beyond 4%. 4.4. Comparison with experiments for different porosities We compared the results obtained from the 2D model against those from experiments reported in Sehaqui et al. (2010) for 3 different porosities viz. 98.10%, 95.90% and 93.10% in Fig. 16. We assumed the elastic modulus of the cell walls to be equal to 20 GPa, which is probably the upper
estimate based on the tests made on the dense NFC membranes (Kulachenko et al., 2012). The model captures qualitatively the observed trend in experiments but there is only an order of magnitude agreement in the stress values. The gradual stiffening after the initial elastic response observed in the experiments is captured by the twodimensional model with internal contact. As no material non-linearity was introduced in the model of the cell walls, the amount of stiffening captured in the model after the initial elastic response corresponds only to geometrical non-linearity. In a relative scale it captures the trend observed in the experiments rather well. There could be several reasons behind the large difference in the actual stress values between the model and simulation. First is the error in estimating the wall thickness or the average cell size of the cells at a given porosity. We demonstrated, however, that although these parameters change the stiffening part of the response, they cannot explain the order of observed differences. The second possibility is underestimating the elastic properties of the cellulose wall material. We used however the upper estimate for the isotropic nanocellulose network and the elastic
P. Srinivasa, A. Kulachenko / Mechanics of Materials 80 (2015) 13–26
properties of such networks are not size dependent for the sizes larger than a micron (Kulachenko et al., 2012). The introduction of the non-linear constitutive behaviour and damage can only decrease the calculated stress value and move it even further away from the experimental values. The only two remaining possibilities are (1) inadequate description of the geometry and (2) the inability of the considered 2D foam representation to capture the response of the 3D structure. The geometry of the NFC foam inner structure can be affected by the dehydration process or by gravity which can both contribute to the creation of the preferred orientation of the struts. Since the considered foams were traditionally tested in the direction normal to the free surface (owing to the difficulties in sample preparation), it is impossible to rule out such a possibility for the available experimental data. The 2D foam representation can potentially introduce an error in capturing the response however it is unlikely that it will make it stiffer as required to match the experimental results since open cell foams gives the upper estimate of the stiffness and crushing strength at a given density (Gibson and Ashby, 1982). This means that anisotropy in inner structure is most plausible explanation of the observed difference. It can be tested by ether reconstructing the inner structure with, for example, micro-CT equipment or by simply performing the tests in various direction of the sample. 5. Summary and conclusions We studied the compressive response of irregular foam structures at high strains using 2D beam representation with internal contact. We have proved the effectiveness of beam element as opposed to solid element in modelling high porosity structures by reproducing the experimental and numerical results of a previously published work, using both beam and solid elements and concluding that beam elements allow for capturing the local bucking of the cell walls while also being computationally more efficient. The size effect for a given relative density was explored by increasing the linear size and the number of cells simultaneously. After finding the size that provides statistically acceptable results up to nearly 70% compressive strains, we explored the effect of thickness at a given porosity and the effect of porosity itself and related it to the experimental data. For every numerical experiment, we single out the effect of inner contact by performing the simulations with and without contact. Internal contact is shown to be of significance in simulating the densification regime of random structures and is the only way to capture the stiffening response at the compressed state. The cell walls thickness has a significant influence on the general trend in the stress-plateau region. An increase in thickness is seen to increase the slope of the stress-plateau and diminishes the difference brought by the inner contact as the total free surface reduces for the given porosity. The non-straightness of the cell walls significantly affects the initial elastic response, smoothening the
25
transition between the initial elastic region and the plateau. The non-straightness has however a negligible influence on the stress levels reached in the plateau region and subsequent stiffening. The random periodic two dimensional models also capture the trend of the compressive stress–strain response seen in experimental results from NFC foams. While there is an order magnitude agreement in the simulated response of this random model in comparison with the experiments of NFC foams, it is seen that neither the effect of stiffness nor the size dependency can explain the difference between the experimental data and simulation results. Although the considered NFC foam predominantly consists of open cell structures, it has non-trivial out of plane connectivity which requires 3D representation. Acknowledgements Financial support from The Swedish Research Council is gratefully acknowledged by the authors. P. Srinivasa wishes to acknowledge Dr. Houssine Sehaqui for providing the data from experiments on compression of foams. References Ali, Z.M., Gibson, L.J., 2013. The structure and mechanics of nanofibrillar cellulose foams. Soft Matter 9, 1580. http://dx.doi.org/10.1039/ c2sm27197d. Alsayednoor, J., Harrison, P., Guo, Z., 2013. Large strain compressive response of 2-D periodic representative volume element for random foam microstructures. Mech. Mater. 66, 7–20. http://dx.doi.org/ 10.1016/j.mechmat.2013.06.006. Aulin, C., Netrval, J., Wågberg, L., Lindström, T., 2010. Aerogels from nanofibrillated cellulose with tunable oleophobicity. Soft Matter 6, 3298. http://dx.doi.org/10.1039/c001939a. Chen, C., Lu, T.J., Fleck, N.A., 1999. Effect of imperfections on the yielding of two-dimensional foams. J. Mater. Chem. 47, 2235–2272. Gibson, L.J., Ashby, M.F., 1982. The mechanics of three-dimensional cellular materials. Proc. R. Soc. A Math. Phys. Eng. Sci. 382, 43–59. http://dx.doi.org/10.1098/rspa.1982.0088. Gibson, L.J., Ashby, M.F., 1988. Cellular Solids. Pergamon press, New York. Gong, L., Kyriakides, S., Jang, W.-Y., 2005. Compressive response of opencell foams. Part I: Morphology and elastic properties. Int. J. Solids Struct. 42, 1355–1379. http://dx.doi.org/10.1016/ j.ijsolstr.2004.07.023. Klemm, D., Kramer, F., Moritz, S., Lindström, T., Ankerfors, M., Gray, D., Dorris, A., 2011. Nanocelluloses: a new family of nature-based materials. Angew. Chem. Int. Ed. Engl. 50, 5438–5466. http:// dx.doi.org/10.1002/anie.201001273. Kulachenko, A., Uesaka, T., 2012. Direct simulations of fiber network deformation and failure. Mech. Mater. 51, 1–14. http://dx.doi.org/ 10.1016/j.mechmat.2012.03.010. Kulachenko, A., Denoyelle, T., Galland, S., Lindström, S.B., 2012. Elastic properties of cellulose nanopaper. Cellulose 19, 793–807. http:// dx.doi.org/10.1007/s10570-012-9685-5. Li, Q.M., Magkiriadis, I., Harrigan, J.J., 2006. Compressive strain at the onset of densification of cellular solids. J. Cell. Plast. 42, 371–392. http://dx.doi.org/10.1177/0021955X06063519. Lindström, S.B., Kulachenko, A., Jawerth, L.M., Vader, D.A., 2013. Finitestrain, finite-size mechanics of rigidly cross-linked biopolymer networks. Soft Matter 9, 7302. http://dx.doi.org/10.1039/ c3sm50451d. Mesarovic, S.D., Padbidri, J., 2005. Minimal kinematic boundary conditions for simulations of disordered microstructures. Philos. Mag. 85, 65–78. http://dx.doi.org/10.1080/14786430412331313321. Nour-Omid, B., Rankin, C.C., 1988. The use of projectors to improve finite element performance. Comput. Struct. 30, 257–267. Nour-Omid, B., Rankin, C.C., 1991. Finite rotation analysis and consistent linearization using projectors. Comput. Methods Appl. Mech. Eng. 93, 353–384. http://dx.doi.org/10.1016/0045-7825(91)90248-5. Okabe, A., Boots, B., Sugihara, K., 1995. Spatial Tessellations: Concepts and Applications of Voronoi Diagrams, first ed. John Wiley & Sons Inc., West Sussex, England.
26
P. Srinivasa, A. Kulachenko / Mechanics of Materials 80 (2015) 13–26
Papka, S.D., Kyriakides, S., 1994. In-plane compressive response and crushing of honeycomb. J. Mech. Phys. Solids 42, 1499–1532. Papka, S.D., Kyriakides, S., 1998. Experiments and full-scale numerical simulations of in-plane crushing of a honeycomb. Acta Mater. 46, 2765–2776. http://dx.doi.org/10.1016/S1359-6454(97)00453-9. Przemieniecki, J.S., 1968. Theory of Matrix Structural Analysis, first ed. McGraw-Hill Book Co., Inc., New York. Roberts, A.P., Garboczi, E.J., 2002. Elastic properties of model random three-dimensional open-cell solids. J. Mech. Phys. Solids 50, 33–55. http://dx.doi.org/10.1016/S0022-5096(01)00056-4. Sehaqui, H., 2011. Nanofiber networks, aerogels and biocomposites based on nanofibrillated cellulose from wood. Doctoral thesis, KTH – Royal Institute of Technology. Sehaqui, H., Salajková, M., Zhou, Q., Berglund, L.A., 2010. Mechanical performance tailoring of tough ultra-high porosity foams prepared from cellulose I nanofiber suspensions. Soft Matter 6, 1824. http:// dx.doi.org/10.1039/b927505c. Silva, M.J., Hayes, W.C., Gibson, L.J., 1995. The effects of non-periodic microstructure on the elastic properties of two-dimensional cellular solids. Int. J. Mech. Sci. 37, 1161–1177. Svagan, A.J., Samir, M.A.S.A., Berglund, L.A., 2008. Biomimetic foams of high mechanical performance based on nanostructured cell walls
reinforced by native cellulose nanofibrils. Adv. Mater. 20, 1263–1269. http://dx.doi.org/10.1002/adma.200701215. Svagan, A.J., Jensen, P., Dvinskikh, S.V., Furó, I., Berglund, L.A., 2010. Towards tailored hierarchical structures in cellulose nanocomposite biofoams prepared by freezing/freeze-drying. J. Mater. Chem. 20, 6646. http://dx.doi.org/10.1039/c0jm00779j. Vanderburg, M.W.D., Shulmeister, V., Van Der Giessen, E., Marrisen, R., 1997. On the linear elastic properties of regular and random open-cell foam models. J. Cell. Plast. 33, 31–54. http://dx.doi.org/10.1177/ 0021955X9703300103. Zhu, H.X., Windle, A.H., 2002. Effects of cell irregularity on the high strain compression of open-cell foams. Acta Mater. 50, 1041–1052. http:// dx.doi.org/10.1016/S1359-6454(01)00402-5. Zhu, H.X., Hobdell, J.R., Windle, A.H., 2000. Effects of cell irregularity on the elastic properties of open-cell foams. Acta Mater. 48, 4893–4900. http://dx.doi.org/10.1016/S1359-6454(00)00282-2. Zhu, H.X., Hobdell, J.R., Windle, A.H., 2001a. Effects of cell irregularity on the elastic properties of 2D Voronoi honeycombs. J. Mech. Phys. Solids 49, 857–870. Zhu, H.X., Thorpe, S.M., Windle, A.H., 2001b. The geometrical properties of irregular two-dimensional Voronoi tessellations. Philos. Mag. 81, 2765–2783.