Modeling the Nonlinear Deformation of Highly Porous Cellular Plastics Filled with Clay Nanoplatelets

Abstract

Rigid low-density plastic foams subjected to mechanical loads typically exhibit a nonlinear deformation stage preceding failure. At moderate strains, when the geometrical nonlinearity is negligible, such foam response is predominantly caused by the nonlinearity of deformation of their principal structural elements—foam struts. Orientational averaging of stresses in foam struts enables estimation of the stresses taken up by foams at a given applied strain. Based on a structural model of highly porous anisotropic cellular plastics filled with clay nanoplatelets and the orientational averaging, a method for calculating their nonlinear deformation is derived in terms of structural parameters of the porous material, the mechanical properties of the monolithic polymer, and filler particles and their spatial orientation. The method is applied to predicting the tensile stress-strain diagrams of organoclay-filled low-density rigid polyurethane foams, and reasonable agreement with experimental data is demonstrated.

1. Introduction

The aim of this research was to develop a structural model for calculating the nonlinear deformation of highly porous cellular plastics with a porosity Π (0 < Π < 1) exceeding 90%. Such cellular plastics (i.e., plastic foams) are used in various fields of engineering, e.g., as damping and heat- and sound-insulating materials. They are not meant to carry high mechanical loads, but their stiffness and strength have to be sufficiently high for intended applications. Filling the plastic foams with micro- or nanoparticles enables improving their functional properties without interfering adversely with foam morphology.

When anisometric filler particles are used in production of low-density composite foams, they tend to align with the principal flow direction of the liquid chemical formulation during foaming as well as with the stretch directions of cell walls and struts, as demonstrated for polymer foams with rod- [1,2,3,4,5,6] and plate-like [7,8,9] mico- and nanofillers. Interactions between dispersed carbon fibers and growing cells in high-pressure foam injection molding experiments using a polystyrene/carbon fiber/carbon dioxide system was investigated in [3], by an in situ visualization technique. It was found that the fibers in close proximity to the growing cells exhibited both translational and rotational displacements. An analytical model was developed in [4] to describe the instantaneous location and angle of rod-like conductive fillers as affected by cell growth during the foaming of conductive polypropylene composites. The reorientation of nanoclay particles during the foaming process has been considered, and its effect on the stiffness of the cell walls was modeled in [9].

The nonlinear mechanical response of such low-density polymer foams can be described by either applying a constitutive material model needing calibration based on foam tests (see, e.g., in [10]), or by a structural model relating foam properties to their morphology and mechanical properties of the respective monolithic polymer. In low-density closed-cell polyurethane (PUR) foams, most of the polymer (more than 90%, see, e.g., in [11]) is contained within cell struts. Such morphology enables accurate prediction of the elastic properties of anisotropic foams by considering the deformation of a regular elementary cell composed of struts, such as parallelepiped [12,13], Kelvin [14,15], and generalized Kelvin [16,17,18] cells. Alternatively, models based on orientation distribution of foam struts are used [19,20]. Concerning the nonlinear deformation, the behavior of a large network of cells under compression has been simulated numerically by a finite element analysis employing either regular [21] or disordered cell geometry. In the latter case, both random generated cellular microstructure [22] and the actual morphology obtained by 3D scanning of foams [23] have been considered. The application of strut-based models also has been mostly confined to the compressive response, when the nonlinearity is caused by a progressive buckling of the struts [19,24]. In this study, we consider the effect of nonlinear deformation of foam strut material on the foam response.

The aim of the present study is to create a mathematical model for describing the nonlinear deformation of highly porous cellular plastics filled with clay nanoplatelets, taking into account the influence of the spatial alignment of filler particles in foam struts. To this end, we utilize the method developed in [20,25,26] for predicting the elastic properties of foams by means of orientational averaging of the linear elastic characteristics of their structural elements, and extend it to nonlinear deformation of structural elements. The paper is organized as follows. Load-carrying and structural elements of low-density foams, calculation of stiffness tensor of the composite foam strut material, and evaluation of the nonlinear deformation of foams by orientation averaging of the stress–strain response of the structural elements are presented in Section 2. Analytical approximations of the nonlinear axial stress–strain relation of the load-carrying element of foams are described in Section 3, and application of the model to prediction of the mechanical response in uniaxial tension of rigid low-density PUR foams with different loadings of nanoclay filler is demonstrated in Section 4.

2. Structural Model

We will consider a cellular plastic of porosity Π ≥ 90% filled with clay particles in the form of nanoplatelets. Owing to the high porosity, we assume that all mass of the polymer is concentrated in rod-like struts and knots of the foams, as seen in the micrograph Figure 1a, and cell walls are too thin to take up a significant mechanical load. The effect of gas pressure in the closed cells on foam deformation can be neglected at small strains [20]. Then the external loads are taken up only by the struts and knots.

The load-carrying element of such a cellular plastic, as in [20,25], is taken in the form of a straight shaft with thickened end parts, depicted in Figure 1b,c. The cross section of the shaft is an equilateral triangle of side $a_s$. The thickened parts are rectangular parallelepipeds with a side $a_k$ and height $l_k/2$. The load-carrying elements are connected to each other by their ends in knots and form a spatial structure able to resist various mechanical actions. This element is placed in a coaxial rectangular parallelepiped of length $l=l_s+l_k$ and a square cross section with a side $a={\left[\frac{V_l}{\left(1-\Pi \right)l}\right]}^{\frac{1}{2}}$, where $V_l$ is volume of the load-carrying element. In such a way, we arrive at the structural element of the foam, Figure 1b.

The load-carrying element takes up axial forces P, transverse forces Q, and moments $M_A=M_B=Ql/2$ (see Figure 1c), and, as a result, it changes its length and transverse dimensions, undergoes shear deformations, and bends. The transverse deformation of the load-carrying element practically does not affect that of the structural element, but bending loads add to it a supplementary shear deformation. Considering that the second moments of cross-sectional area of the triangular shaft about the $x_2^{\prime }$ and $x_3^{\prime }$ axes coincide, the stiffness tensor ${\stackrel{*}{\mathit{C}}}^{lce}$ of the load-carrying element has only three nonzero components in the $x_i^{\prime }$ axes—${\stackrel{*}{C}}_{1111}^{lce}$ and ${\stackrel{*}{C}}_{1212}^{lce}={\stackrel{*}{C}}_{1313}^{lce}$. Their magnitude depends on the geometrical dimensions of the load-carrying element and the elastic characteristics of its material [20,25,26].

2.1. Stiffness of the Composite Material of the Load-Carrying Element

The polymer material of load-carrying elements is uniformly filled with clay nanoplatelets around the longitudinal axis $x_1^{\prime }$ with a volume fraction $V_f$. As a result, the element is transversely isotropic. The stiffness tensor $\mathit{C}$ of its material is found by spatially averaging the transversely isotropic stiffness tensor $\stackrel{*}{\mathit{C}}$ of the calculation element (see, e.g., in [27])

$$\mathit{C}={{\displaystyle \int }}_{O}f\left(g\right)\stackrel{*}{\mathit{C}}\left(g\right)dg, g\in O,$$

(1)

where O is the group of 3D orientations of the elements. Their orientation distribution function $f\left(g\right)$ is taken in the form:

$$f\left(g\right)=f\left(\mathbf{\xi }\right)=\frac{2n+1}{{\lambda }_1+{2\lambda }_{23}}\left[{\lambda }_1{\xi }_1^{2n}+2\frac{\left(2n-1\right)!!}{\left(2n\right)!!}{\lambda }_{23}{(1-{\xi }_1^2)}^n\right], n\ge 0$$

(2)

where ${\lambda }_1$ and ${\lambda }_{23}$ are nonnegative parameters such that ${\lambda }_1+{\lambda }_{23}=1$. The unit vector $\mathbf{\xi }$ indicates the spatial directions of calculation elements, with ${\xi }_1=\cos \theta$ , where $\theta$ is the angle between it and the axis $x_1^{\prime }$ of the load-carrying element, see Figure 1d.

For example, if ${\lambda }_{23}$ = 0, then ${\lambda }_1=1$ and the distribution function Equation (2) has only one parameter, $n$; performing integration (1), we have

$$\begin{array}{c}\begin{array}{c}{C}_{1111}=k\left[\left(2n+3\right)\left(2n+1\right){\stackrel{*}{C}}_{1111}+3\left({\stackrel{*}{C}}_{2222}+{\stackrel{*}{C}}_{3333}\right)+2{\stackrel{*}{C}}_{2233}+4{\stackrel{*}{C}}_{2323}+\right.\\ \left.2\left(2n+1\right)\left({\stackrel{*}{C}}_{1122}+{\stackrel{*}{C}}_{3311}+2{\stackrel{*}{C}}_{1212}+2{\stackrel{*}{C}}_{3131}\right)\right],\end{array}\\ \begin{array}{c}{C}_{2222}=C_{3333}=k\left[3{\stackrel{*}{C}}_{1111}+\frac{1}{2}\left(n+2\right)\left(n+1\right)(3{\stackrel{*}{C}}_{2222}+3{\stackrel{*}{C}}_{3333}+2{\stackrel{*}{C}}_{2233}+\right.\\ \left.4{\stackrel{*}{C}}_{2323})+2\left(n+1\right)\left({\stackrel{*}{C}}_{1122}+{\stackrel{*}{C}}_{3311}+2{\stackrel{*}{C}}_{1212}+2{\stackrel{*}{C}}_{3131}\right)\right],\end{array}\\ \begin{array}{c}{C}_{1122}=C_{3311}=k\left[\left(2n+1\right)\left({\stackrel{*}{C}}_{1111}-2{\stackrel{*}{C}}_{1212}+{\stackrel{*}{C}}_{3131}\right)+\left(n+1\right)({\stackrel{*}{C}}_{2222}+\right.\\ \left.{\stackrel{*}{C}}_{3333})+2\left(n+2\right){\stackrel{*}{C}}_{2233}-2{\stackrel{*}{C}}_{2323}+\left(2n^2+5n+4\right)\left({\stackrel{*}{C}}_{1122}+{\stackrel{*}{C}}_{3311}\right)\right],\end{array}\\ \begin{array}{c}{C}_{2233}=k\left[{\stackrel{*}{C}}_{1111}-2{\stackrel{*}{C}}_{1212}-2{\stackrel{*}{C}}_{3131}+\frac{1}{2}\left(n+2\right)\left(n+1\right)\left({\stackrel{*}{C}}_{2222}+{\stackrel{*}{C}}_{3333}\right)+\right.\\ \left(3n^2+9n+4\right){\stackrel{*}{C}}_{2233}-2\left(n+3n+1\right){\stackrel{*}{C}}_{2323}+2\left(n+2\right)\left.\left({\stackrel{*}{C}}_{1122}+{\stackrel{*}{C}}_{3311}\right)\right],\end{array}\\ \begin{array}{c}{C}_{1212}=C_{3131}=k\left[\left(2n+1\right)\left({\stackrel{*}{C}}_{1111}-{\stackrel{*}{C}}_{1122}-{\stackrel{*}{C}}_{3311}\right)+\left(n+1\right)({\stackrel{*}{C}}_{2222}+\right.\\ {\stackrel{*}{C}}_{3333})-{\stackrel{*}{C}}_{2233}+\left(2n+3\right){\stackrel{*}{C}}_{2323}+\left(2n^2+3n+3\right)\left.\left({\stackrel{*}{C}}_{1212}+{\stackrel{*}{C}}_{3131}\right)\right]\end{array}\\ \begin{array}{c}{C}_{2323}=k\left[{\stackrel{*}{C}}_{1111}-{\stackrel{*}{C}}_{1122}-{\stackrel{*}{C}}_{3311}+\frac{1}{2}\left(n+2\right)\left(n+1\right)\left({\stackrel{*}{C}}_{2222}+{\stackrel{*}{C}}_{3333}\right)\right.+\\ \left(n^2+3n+1\right){\stackrel{*}{C}}_{2233}+\left(2n+6n+3\right){\stackrel{*}{C}}_{2323}+\left(2n+3\right)\left.\left({\stackrel{*}{C}}_{1212}+{\stackrel{*}{C}}_{3131}\right)\right],\end{array}\\ k=\frac{1}{\left(2n+5\right)\left(2n+3\right)}.\end{array}$$

(3)

Using the Halpin–Pagano–Tsai equations modified for the case of platelike filler particles [28], the stiffness ${\stackrel{*}{C}}_{ijkl}$ is expressed in the form

$$\begin{array}{c}{\stackrel{*}{C}}_{1111}=\frac{1-{\stackrel{*}{\nu }}_{23}^2}{{\stackrel{*}{E}}_2^2{\stackrel{*}{\mathsf{\Delta }}}_1}, {\stackrel{*}{C}}_{2222}={\stackrel{*}{C}}_{3333}=\frac{1-{\stackrel{*}{\nu }}_{12}{\stackrel{*}{\nu }}_{21}}{{\stackrel{*}{E}}_1{\stackrel{*}{E}}_2{\stackrel{*}{\mathsf{\Delta }}}_1},\\ {\stackrel{*}{C}}_{1122}={\stackrel{*}{C}}_{3311}=\frac{{\stackrel{*}{\nu }}_{21}(1+{\stackrel{*}{\nu }}_{23})}{{\stackrel{*}{E}}_1{\stackrel{*}{E}}_2{\stackrel{*}{\mathsf{\Delta }}}_1}=\frac{{\stackrel{*}{\nu }}_{12}(1+{\stackrel{*}{\nu }}_{23})}{{\stackrel{*}{E}}_2^2{\stackrel{*}{\mathsf{\Delta }}}_1}, {\stackrel{*}{C}}_{2233}=\frac{{\stackrel{*}{\nu }}_{23}+{\stackrel{*}{\nu }}_{12}{\stackrel{*}{\nu }}_{21}}{{\stackrel{*}{E}}_1{\stackrel{*}{E}}_2{\stackrel{*}{\mathsf{\Delta }}}_1},\\ {\stackrel{*}{C}}_{1212}={\stackrel{*}{C}}_{3131}={\stackrel{*}{G}}_{12}={\stackrel{*}{G}}_{31}, {\stackrel{*}{C}}_{2323}={\stackrel{*}{G}}_{23}=\frac{{\stackrel{*}{E}}_2}{2(1+{\stackrel{*}{\nu }}_{23})},\\ {\stackrel{*}{\mathsf{\Delta }}}_1=\frac{1}{{\stackrel{*}{E}}_1{\stackrel{*}{E}}_2^2}\left(1+{\stackrel{*}{\nu }}_{23}\right)\left(1-{\stackrel{*}{\nu }}_{23}-2{\stackrel{*}{\nu }}_{12}{\stackrel{*}{\nu }}_{21}\right),\end{array}$$

(4)

with

$$\begin{array}{c}{\stackrel{*}{E}}_1=\frac{1+{\xi }^{\prime }{\eta }^{\prime }{V}_f}{1-{\eta }^{\prime }{V}_f}{E}_m, {\stackrel{*}{E}}_2={\stackrel{*}{E}}_3=\frac{1+{\xi }^{″}{\eta }^{″}{V}_f}{1-{\eta }^{″}{V}_f}{E}_m,\\ {\stackrel{*}{G}}_{12}={\stackrel{*}{G}}_{31}=\frac{1+{\chi }^{\prime }{\lambda }^{\prime }{V}_f}{1-{\lambda }^{\prime }{V}_f}{G}_m, {\stackrel{*}{G}}_{23}=\frac{1+{\chi }^{″}{\lambda }^{″}{V}_f}{1-{\lambda }^{″}{V}_f}{G}_m,\\ {\stackrel{*}{\nu }}_{12}={\stackrel{*}{\nu }}_{13}={\nu }_f{V}_f+{\nu }_m\left(1-V_f\right).\end{array}$$

(5)

Here,

$$\begin{array}{c}{\eta }^{\prime }=\frac{E_f/E_m-1}{E_f/E_m+{\xi }^{\prime }}, {\eta }^{″}=\frac{E_f/E_m-1}{E_f/E_m+{\xi }^{″}},\\ {\lambda }^{″}=\frac{G_f/G_m-1}{G_f/G_m+{\chi }^{″}}, {\lambda }^{″}=\frac{G_f/G_m-1}{G_f/G_m+{\chi }^{″}},\\ {\xi }^{\prime }=2, {\chi }^{\prime }=1, {\xi }^{″}={\chi }^{″}=\frac{2D}{h}\end{array}$$

(6)

In the above relations, $V_f$ is the volume fraction of filler (clay nanoplatelets), ${\nu }_f$ and ${\nu }_m$ are the Poisson ratios, $E_f$ and $E_m$ are the elastic moduli, and $G_f$ and $G_m$ are the shear moduli of the filler and polymer matrix, respectively. D and h are the diameter and thickness, respectively, of clay nanoplatelets.

2.2. Stresses in the Structural Element

The stresses $\stackrel{*}{{\sigma }_{11}}$, $\stackrel{*}{{\sigma }_{12}}$, and $\stackrel{*}{{\sigma }_{13}}$ in the structural element, considering the porosity Π of the foam, are found as follows:

$$\begin{array}{c}\stackrel{*}{{\sigma }_{11}}={(1-\Pi )}^{\alpha }\frac{1-{\stackrel{*}{\nu }}_{12}}{\left(1+{\stackrel{*}{\nu }}_{12}\right)\left(1-2{\stackrel{*}{\nu }}_{13}\right)}{\left(\frac{l_k}{l}\frac{F_l}{F_k}+\frac{l_s}{l}\frac{F_l}{F_s}\right)}^{-1}\stackrel{*}{\sigma },\\ \stackrel{*}{{\sigma }_{12}}=\left(1-\Pi \right){\stackrel{*}{G}}_0\stackrel{*}{{\epsilon }_{12}}, \stackrel{*}{{\sigma }_{13}}=\left(1-\Pi \right){\stackrel{*}{G}}_0\stackrel{*}{{\epsilon }_{13}},\end{array}$$

(7)

where $\stackrel{*}{\sigma }=f\left({\stackrel{*}{\epsilon }}_{11}\right)$ is the stress of the material of load-carrying element (i.e., clay nanoplatelet-filled polymer) as a function of strain in loading along the axis $x_1^{\prime }$, and [20,25]

$${\stackrel{*}{G}}_0=\frac{2}{{\stackrel{*}{E}}_1}·\frac{l_s}{l}{\left(\frac{l_s}{a_s}\right)}^2\frac{F_l}{F_s}+\frac{\left[1+{\left(\frac{l_s}{l}\right)}^3\right]\frac{F_l}{F_s}}{{\left(\frac{a_k}{l_k}·\frac{l_k}{l}\right)}^2{\stackrel{*}{E}}_1}+\frac{\chi }{{\stackrel{*}{E}}_1}+\frac{1}{{\stackrel{*}{G}}_{12}}{\left(\frac{l_k}{l}·\frac{F_l}{F_k}+\frac{l_s}{l_k}·\frac{F_l}{F_s}\right)}^{-1}.$$

(8)

In Equation (8), ${\stackrel{*}{E}}_1$ and ${\stackrel{*}{G}}_{12}$ are the elastic and shear moduli of the platelet-filled material of the load-carrying element (determined from stiffness tensor $\mathit{C}$ of the composite given by Equation (1)) and

$$\begin{array}{c}\frac{F_l}{F_k}=\left[\frac{\sqrt{3}}{4}{\left(\frac{a_s}{l_s}\right)}^2{\left(\frac{l_s}{l}\right)}^3+{\left(\frac{a_k}{l_k}\right)}^2{\left(\frac{l_k}{l}\right)}^3\right]{\left(\frac{a_k}{l_k}·\frac{l_k}{l}\right)}^{-2},\\ \frac{F_l}{F_s}=\left[{\left(\frac{a_s}{l_s}\right)}^2{\left(\frac{l_s}{l}\right)}^3+\frac{4}{\sqrt{3}}{\left(\frac{a_s}{l_k}\right)}^2{\left(\frac{l_k}{l}\right)}^3\right]{\left(\frac{a_s}{l_s}·\frac{l_s}{l}\right)}^{-2},\\ \frac{a_k}{a_s}=\frac{\frac{a_k}{l_s}}{\frac{l_s}{l}·\frac{a_s}{l_s}}, \frac{l_k}{l_s}=\frac{\frac{l_k}{l}}{1-\frac{l_k}{l}}, \frac{l_s}{l}=1-\frac{l_k}{l},\\ \chi =\frac{1}{2}{\chi }_0{\stackrel{*}{E}}_1{V}_l, {\chi }_0\ge 0.\end{array}$$

(9)

It is seen that the stresses $\stackrel{*}{{\sigma }_{11}}$, $\stackrel{*}{{\sigma }_{12}}$, and $\stackrel{*}{{\sigma }_{13}}$ depend on the non-dimensional parameters $\frac{a_s}{l_s}$, $\frac{a_k}{l_k}$, $\frac{l_k}{l}$, $\chi$, and $\Pi =1-\frac{\rho }{{\rho }_m}$, where $\rho$ and ${\rho }_m$ are densities of the porous material and matrix, respectively, the power $\alpha$ is introduced to consider the fact that, under linear deformations, the open-cell foam can change its volume and, thus, the density $\rho$.

Taking into account that approximately $\frac{a_k}{a_s}\approx$ 1.5 and $a_k\approx l_k$, we obtain that

$$\frac{l_k}{l}\approx \frac{\frac{a_s}{l_s}}{\frac{a_s}{l_s}+0.77}, \frac{l_s}{l}\approx \frac{0.77}{\frac{a_s}{l_s}+0.77}.$$

The parameter $\frac{a_s}{l_s}$ is found by numerically solving the nonlinear equation [20,25]

$$1-\Pi =\frac{\rho }{{\rho }_m}=\frac{{\left(\frac{a_s}{l_s}\right)}^2}{{\left(\frac{a_s}{l_s}+0.77\right)}^3}\left(\frac{3\sqrt{3}}{4}{0.77}^3+\frac{a_s}{l_s}\right).$$

2.3. Evaluation of Foam Stresses and Foam Stiffness

Now that the composition of the structural element of foam is known, the stresses ${\sigma }_{ij}$ of the deformed porous material are calculated by orientationally averaging the stresses $\stackrel{*}{{\sigma }_{11}}$, $\stackrel{*}{{\sigma }_{12}}$, and $\stackrel{*}{{\sigma }_{13}}$, Equation (7), of the transversely isotropic structural elements with account of their orientational distribution function $f_p\left(\mathbf{\xi }\right)$.

Assuming that [20]

$$\begin{gathered} f_p\left(\mathbf{\xi }\right)=\frac{1}{k}\left[1+\left(k^3-1\right){\xi }_1^{2n}\right], {\xi }_1=\cos \theta , \\ n=\frac{k\left(k+1\right)}{2}, k>0, \end{gathered}$$

(10)

where $k$ is the extension degree of pores, we write that

$$\begin{array}{ll}{\sigma }_{ij}=& {{\displaystyle \int }}_O{f}_p\left(\mathbf{\xi }\right)\left[l_{1i}{l}_{ij}\stackrel{*}{{\sigma }_{11}}+\left(l_{1i}{l}_{2j}+l_{1j}{l}_{2i}\right)\stackrel{*}{{\sigma }_{12}}\right]dg\\ & \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=\frac{1}{S}{\displaystyle \int }{f}_p\left(\mathbf{\xi }\right)\left\{l_{1i}{l}_{ij}\stackrel{*}{{\sigma }_{11}}\right.\\ & \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+\left.\frac{1}{2}\left[\left(l_{1i}{l}_{2j}+l_{1j}{l}_{2i}\right)\stackrel{*}{{\sigma }_{12}}+\left(l_{1i}{l}_{3j}+l_{1j}{l}_{3i}\right)\stackrel{*}{{\sigma }_{13}}\right]\right\}ds,\end{array}$$

(11)

where $\stackrel{*}{{\sigma }_{11}}$ is a nonlinear function, calculated by Equation (7) using $\stackrel{*}{\sigma }=f\left({\stackrel{*}{\epsilon }}_{11}\right)$ found from experiments,

$$\begin{gathered} \stackrel{*}{{\sigma }_{12}}=2{\stackrel{*}{C}}_{1212}^p{\stackrel{*}{\epsilon }}_{12}, \stackrel{*}{{\sigma }_{13}}=2{\stackrel{*}{C}}_{1313}^p{\stackrel{*}{\epsilon }}_{13}, {\stackrel{*}{C}}_{1212}^p={\stackrel{*}{C}}_{1313}^p, \\ {l}_{ij}=\cos \left(x_i^{\prime },x_j\right), ds=\frac{1}{4\pi }. \end{gathered}$$

The evaluation of integral (11) can be simplified assuming that the Euler angles ${\phi }_2=0$ and ${\phi }_1=\frac{\pi }{2}-\phi$, where $\phi$ is the spherical angle. Then,

$$\begin{array}{c}{l}_{11}=\cos \theta , l_{12}=\cos \phi \sin \theta , l_{13}=\sin \phi \sin \theta , l_{21}=0, l_{22}=\sin \phi ,\\ l_{23}=-\cos \phi , l_{31}=-\sin \theta , l_{32}=\cos \phi \cos \theta , l_{33}=\sin \phi \cos \theta ,\\ 0\le \phi <2\pi , 0\le \theta <\pi ,\end{array}$$

(12)

and in the linear case, we have the following expressions for components of the stiffness tensor ${\mathit{C}}_p$ of the porous material [20]

$$\begin{array}{c}{C}_{1111}^p=\frac{1}{k}\left\{\left[\frac{1}{5}+\frac{\left(2n+1\right)\left(k-1\right)}{2n+5}\right]\right.{\stackrel{*}{C}}_{1111}^p+\left.8\left[\frac{1}{15}+\frac{\left(2n+1\right)\left(k-1\right)}{\left(2n+5\right)\left(2n+3\right)}\right]{\stackrel{*}{C}}_{1212}^p\right\},\\ C_{2222}^p=C_{3333}^p=\frac{1}{k}\left\{\left[\frac{1}{5}+\frac{3\left(k-1\right)}{\left(2n+5\right)\left(2n+3\right)}\right]\right.{\stackrel{*}{C}}_{1111}^p+\left.8\left[\frac{1}{15}+\frac{\left(2n+1\right)\left(k-1\right)}{\left(2n+5\right)\left(2n+3\right)}\right]{\stackrel{*}{C}}_{1212}^p\right\},\\ C_{1122}^p=C_{3311}^p=\frac{1}{k}\left\{\left[\frac{1}{15}+\frac{\left(2n+1\right)\left(k-1\right)}{\left(2n+5\right)\left(2n+3\right)}\right]\left({\stackrel{*}{C}}_{1111}^p-4{\stackrel{*}{C}}_{1212}^p\right)\right\},\\ C_{2233}^p=\frac{1}{k}\left\{\left[\frac{1}{15}+\frac{k-1}{\left(2n+5\right)\left(2n+3\right)}\right]\left({\stackrel{*}{C}}_{1111}^p-4{\stackrel{*}{C}}_{1212}^p\right)\right\},\\ C_{1212}^p=C_{3131}^p=\frac{1}{k}\left\{\left[\frac{1}{15}+\frac{\left(2n+1\right)\left(k-1\right)}{\left(2n+5\right)\left(2n+3\right)}\right]\right.{\stackrel{*}{C}}_{1111}^p+\left.2\left[\frac{1}{5}+\frac{\left(2n^2+3n+3\right)\left(k-1\right)}{\left(2n+5\right)\left(2n+3\right)}\right]{\stackrel{*}{C}}_{1212}^p\right\},\\ C_{2323}^p=\frac{1}{k}\left\{\left[\frac{1}{15}+\frac{k-1}{\left(2n+5\right)\left(2n+3\right)}\right]\right.{\stackrel{*}{C}}_{1111}^p+\left.2\left(\frac{1}{5}+\frac{k-1}{2n+5}\right){\stackrel{*}{C}}_{1212}^p\right\},\end{array}$$

(13)

$C_{2222}^p-C_{2233}^p=2C_{2323}^p$, where ${\stackrel{*}{C}}_{1111}^p$ and ${\stackrel{*}{C}}_{1212}^p$ are the stiffnesses of structural elements,

$$\begin{array}{c}{E}_1^p=\frac{\mathsf{\Delta }}{{\left(C_{2222}^p\right)}^2-{\left(C_{2233}^p\right)}^2}=C_{1111}^p-\frac{2{\left(C_{1122}^p\right)}^2}{C_{2222}^p+C_{2233}^p},\\ E_2^p=E_3^p=\frac{\mathsf{\Delta }}{C_{1111}^p{C}_{2222}^p-{\left(C_{1122}^p\right)}^2},\\ {\nu }_{12}^p={\nu }_{13}^p=\frac{C_{1122}^p\left(C_{2222}^p-C_{2233}^p\right)}{C_{1111}^p{C}_{2233}^p-{\left(C_{1122}^p\right)}^2}, {\nu }_{31}^p={\nu }_{21}^p=\frac{C_{1122}^p}{C_{2222}^p+C_{2233}^p},\\ {\nu }_{23}^p={\nu }_{32}^p=\frac{C_{1111}^p{C}_{2233}^p-{\left(C_{1122}^p\right)}^2}{C_{1111}^p{C}_{2222}^p-{\left(C_{1122}^p\right)}^2},\\ G_{12}^p=G_{31}^p=C_{1212}^p=C_{3131}^p, G_{23}^p=C_{2323}^p=\frac{1}{2}\left(C_{2222}^p-C_{2233}^p\right),\\ \mathsf{\Delta }=\left(C_{2222}^p-C_{2233}^p\right)\left[C_{1111}^p\left(C_{2222}^p+C_{2233}^p\right)-2{\left(C_{1122}^p\right)}^2\right].\end{array}$$

(14)

When the liquid-filled polymer composition has hardened, the porous transversely isotropic material is linearly elastic at small strains; therefore, its load-carrying element is also assumed linearly elastic in this strain range and its stress–strain relation is described by the straight line $\stackrel{*}{\sigma }={\stackrel{*}{E}}_1\stackrel{*}{\epsilon }$, where ${\stackrel{*}{E}}_1$ is the elastic modulus.

3. Nonlinearity

If the load-carrying element is nonlinearly elastic outside the strain interval $[\overline{\epsilon }, \stackrel{+}{\epsilon }]$, its nonlinear deformation is described by an appropriate nonlinear function. For this aim, e.g., parabolic and ellipsoidal functions can be used.

3.1. Parabolic Functions

If $\stackrel{*}{\epsilon }>\stackrel{+}{\epsilon }$, then

$$\stackrel{*}{\sigma }=\stackrel{+}{{\sigma }_0}+\left(\stackrel{+}{\sigma }-\stackrel{+}{{\sigma }_0}\right){\left[{\left(\frac{\stackrel{*}{\epsilon }-\stackrel{+}{{\epsilon }_0}}{\stackrel{+}{\epsilon }-\stackrel{+}{{\epsilon }_0}}\right)}^2\right]}^{\frac{\stackrel{+}{n}}{2}},$$

(15)

if $\stackrel{*}{\epsilon }<\overline{\epsilon }$, then

$$\stackrel{*}{\sigma }=\overline{{\sigma }_0}+\left(\overline{\sigma }-\overline{{\sigma }_0}\right){\left[{\left(\frac{\stackrel{*}{\epsilon +}\overline{{\epsilon }_0}}{\overline{\epsilon }+\overline{{\epsilon }_0}}\right)}^2\right]}^{\frac{\overline{n}}{2}},$$

(16)

$$\left\{\begin{array}{ll}0\le \stackrel{+}{{\epsilon }_0}<\stackrel{+}{\epsilon },& \overline{\epsilon }<\overline{{\epsilon }_0}\le 0,\\ 0\le \stackrel{+}{{\sigma }_0}<\stackrel{+}{\sigma },& \overline{\sigma }<\overline{{\sigma }_0}\le 0.\end{array}\right.$$

The points $\left(\stackrel{\pm }{{\epsilon }_0}, \stackrel{\pm }{{\sigma }_0}\right)$ are vertices of the parabolas. The parabolas join the straight line at the points $\left(\stackrel{\pm }{\epsilon }, \stackrel{\pm }{\sigma }\right)$.

If it is required that the straight line be tangential to the parabolas at the points $\left(\stackrel{\pm }{\sigma }, \stackrel{\pm }{\epsilon }\right)$, then the powers $\stackrel{\pm }{n}$ have to satisfy the relation

$$\stackrel{\pm }{n}=\frac{\stackrel{\pm }{\epsilon }-\stackrel{\pm }{{\epsilon }_0}}{\stackrel{\pm }{\sigma }-\stackrel{\pm }{{\sigma }_0}}{\stackrel{*}{E}}_1.$$

(17)

This case is illustrated in Figure 2.

If $\stackrel{\pm }{{\epsilon }_0}=0$ and ${\stackrel{\pm }{\sigma }}_0=0$, then $\stackrel{\pm }{n }=1$, and the parabolas degenerate into the straight line $\stackrel{*}{\sigma }={\stackrel{*}{E}}_1\stackrel{*}{\epsilon }$. In general, the parameters $\stackrel{\pm }{{\epsilon }_0}$ and $\stackrel{\pm }{{\sigma }_0}$ depend on the volume fraction $V_f$ of clay nanoparticles in the polymer. Then, they can be expressed as

$$\stackrel{\pm }{{\epsilon }_0}={\stackrel{\pm }{\epsilon }}_{00}\exp \left(\stackrel{\pm }{n_{\epsilon }}{V}_f\right), \stackrel{\pm }{{\sigma }_0}={\stackrel{\pm }{\sigma }}_{00}\exp \left(\stackrel{\pm }{n_{\sigma }}{V}_f\right).$$

(18)

3.2. Elliptic Functions

If $|\stackrel{+}{\epsilon }|>|\stackrel{\pm }{\epsilon }|$, then

$$\stackrel{\pm }{\sigma }=\stackrel{\pm }{b_c}\pm \frac{\stackrel{\pm }{b}}{\stackrel{\pm }{a}}\sqrt{{\stackrel{\pm }{a}}^2-{\left(\stackrel{*}{\epsilon }-{\stackrel{\pm }{a}}_c\right)}^2}$$

(19)

are ellipses with semiaxes $\stackrel{\pm }{a}$ and $\stackrel{\pm }{b}$ and centers $({\stackrel{\pm }{a}}_c, {\stackrel{\pm }{b}}_c)$.

They pass through the points $\left(\stackrel{\pm }{\epsilon }, \stackrel{\pm }{\sigma }\right)$ if

$$\stackrel{\pm }{b_c}=\stackrel{\pm }{\sigma }\mp \frac{\stackrel{\pm }{b}}{\stackrel{\pm }{a}}\sqrt{{\stackrel{\pm }{a}}^2-{\left(\stackrel{\pm }{\epsilon }-{\stackrel{\pm }{a}}_c\right)}^2}.$$

(20)

If it is required that the straight line $\stackrel{*}{\sigma }={\stackrel{*}{E}}_1\stackrel{\pm }{\epsilon }$ be the tangent to the ellipses at the points $\left(\stackrel{\pm }{\epsilon }, \stackrel{\pm }{\sigma }\right)$, then the parameters $\stackrel{\pm }{b}$ have to obey the relation

$$\stackrel{\pm }{b}=\stackrel{\pm }{\sigma }\mp \frac{\stackrel{\pm }{a}\sqrt{{\stackrel{\pm }{a}}^2-{\left(\stackrel{\pm }{\epsilon }-{\stackrel{\pm }{a}}_c\right)}^2}}{\stackrel{\pm }{\epsilon }-{\stackrel{\pm }{a}}_c}{\stackrel{*}{E}}_1.$$

(21)

If $V_f>0$, then $\stackrel{\pm }{a}$, ${\stackrel{\pm }{a}}_c$, and ${\stackrel{*}{E}}_1$ are functions of $V_f$, e.g.,

$$\stackrel{\pm }{a}={\stackrel{\pm }{a}}_0\exp \left(\stackrel{\pm }{k_a}{V}_f\right),\hspace{1em}{\stackrel{\pm }{a}}_c={\stackrel{\pm }{a}}_{c\left(0\right)}\exp \left(\stackrel{\pm }{k_c}{V}_f\right).$$

(22)

If the data for nonlinear deformation of a structural element with $V_f>0$ are not available, then the function $\stackrel{*}{\sigma }(\stackrel{*}{\epsilon })$ can be approximately predicted employing the parabolic or elliptic function indicated, assuming that the parameters $\stackrel{\pm }{\epsilon }$ do not depend on $V_f$ and ${\stackrel{\pm }{n}}_{\epsilon }={\stackrel{\pm }{n}}_{\sigma }=0$ in relations (18) or ${\stackrel{\pm }{k}}_a={\stackrel{\pm }{k}}_b=0$ in relations (22) and considering the fact that the value of $V_f$ is typically rather small. Then, $V_f$ will affect mainly the elastic modulus ${\stackrel{*}{E}}_1$.

4. Comparison of Theoretical Calculations with Experimental Data

As an example of nonlinear deformation of a transversely isotropic cellular plastic, we calculated its tension diagram along the $x_2$ axis at ${\epsilon }_{22}>0$, ${\epsilon }_{11}=-{\nu }_{12}{\epsilon }_{22}$, ${\epsilon }_{33}=-{\nu }_{32}{\epsilon }_{22}$, and ${\sigma }_{11}={\sigma }_{33}=0$. Then,

$$\begin{array}{c}\stackrel{*}{{\epsilon }_{11}}=\left(-{\nu }_{12}{l}_{11}^2+l_{12}^2-{\nu }_{32}{l}_{13}^2\right){\epsilon }_{22},\\ \stackrel{*}{{\epsilon }_{13}}=\left(l_{12}{l}_{22}-{\nu }_{32}{l}_{23}\right){\epsilon }_{22},\\ \stackrel{*}{{\epsilon }_{12}}=0.\end{array}$$

(23)

Young’s modulus $E_m$ = 2300 MPa and the nonlinear deformation parameters, Equation (19), of the monolithic polymer were obtained from experimental stress–strain diagrams of the neat PUR polymer in tension [29,30], while the Poisson’s ratio was taken as ${\nu }_m$ = 0.32 [31]. The volume fraction $V_f$ of the filler was evaluated from the weight fraction W_f as $V_f=W_f{\rho }_m/\left[W_f{\rho }_m+\left(1-W_f\right){\rho }_f\right]$, using the densities of monolithic PUR ${\rho }_m$ reported in [29] and montmorillonite clay ${\rho }_f$ [32]. For filled foams, it was assumed that Cloisite^® 30B clay filler had fully exfoliated during the rigid PUR foam production; geometrical dimensions and mechanical characteristics of the clay nanoplatelets used in calculation were taken from [31]. Introduction of clay filler led to a slight reduction of cell size, while the geometrical anisotropy of the foam cells remained the same within the experimental scatter [33] (see also SEM images of neat [30] and filled [33] foams), therefore the average value of cell shape anisotropy ratio (i.e., the degree of cell extension) k = 1.5 [30,33] was used in modeling.

Predicted deformation diagrams are compared with the experimental data, reported in [34], of neat and Cloisite^® 30B clay-filled foams with the same composition of the PUR matrix as the monolithic polymer [30]. The formulation and manufacturing procedure of foams, as well as preparation and testing of foam specimens, are described in [34]. For completeness and ease of reference, we briefly recapitulate the relevant information. Foam specimens of dog-bone shape, with a rectangular test section of 85 mm length, 22 mm width and 20 mm thickness, where cut from slices of the free-rise foam blocks so that the mechanical response in the direction normal to foam rise could be characterized. Metallic plates with hooks were glued to the ends of the specimens, and short chains were attached to the hooks to enable gripping. Such a gripping system [29] provided alignment of the specimen with the line of action of the applied load, eliminating bending and ensuring pure tension during the test. A clip-on extensometer with 50 mm base length was used for strain measurement in the loading direction. Tensile tests were performed at a displacement rate of 8 mm/min.

The test results in tension transverse to the foam rise direction for foams of ~40 kg/m³ density, with 0, 1, 2, and 5 wt.% loading of the clay filler [34], are shown in Figure 3 together with the model prediction obtained using the elliptic nonlinear function (19) with ${\stackrel{\pm }{k}}_a={\stackrel{\pm }{k}}_b=0$. Three foam specimens were tested for each filler loading level; the respective experimental stress–strain diagrams are plotted by dashed lines in Figure 3.

Two cases of the spatial orientation of clay nanoplatelets are considered in modeling—random orientation and perfect alignment of the platelets along strut axis. The alignment occurs during foaming, due to suction of the liquid chemical formulation and filler particles through the Gibbs’ channels between the growing bubbles and stretching of the filled polymeric struts along their longitudinal axis [35,36].

Numerical calculations showed that, as expected, alignment of the anisometric filler in the axial direction of foam struts increased foam stiffness compared with the case of random orientation, see Figure 3b–d. The effect of alignment became more pronounced with increasing clay loading. The same qualitative tends are seen also in Figure 4 presenting comparison of test results of rigid PUR foams with a slightly greater, ~50 kg/m³, density and 0, 1, 2, and 3 wt.% loading of Cloisite^® 30B [34] with theoretical prediction performed as described above.

TEM microscopy has revealed that clay nanoplatelet-filled cellular plastics, e.g., polypropylene foams [6], represent an intermediate case, when a certain amount of the anisometric nanoparticles remains randomly oriented in the central part of load-carrying elements, especially knots. The stress–strain diagrams of such materials are expected to lie between the diagrams of materials comprising irregularly oriented and perfectly aligned nanoplatelets.

It is seen in Figure 3b–d that the calculated ${\sigma }_{22}\left({\epsilon }_{22}\right)$ curves for filled foams tend to pass above the experimental diagrams. Such an overestimation of foam stiffness by the model is likely to be caused by a partial exfoliation of the clay contained in foam struts, as opposed to complete exfoliation assumed in calculations. It has been shown in, e.g., [37] that, for a given loading of clay filler, reduction in the fraction of exfoliated particles leads to a substantial decrease of the stiffness of the filled composite. In addition, the nonlinearity of deformation appears overestimated at the higher loadings of the filler, Figure 3c,d and Figure 4d. This can be explained by the fact that, when describing nonlinear deformation of foams, the effect of $V_f$ was taken into account only approximately due to the simplifying assumption used that ${\stackrel{\pm }{k}}_a={\stackrel{\pm }{k}}_b=0$. To fully reflect in the model the effect of filler on the nonlinearity of foam deformation, it is necessary to determine the coefficients ${\stackrel{\pm }{k}}_a$ and ${\stackrel{\pm }{k}}_b$ in relations (22) for some non-zero value of $V_f$, thus calibrating the model.

5. Conclusions

A structural model for a nonlinearly deforming highly porous cellular plastic filled with clay nanoplatelets is derived taking into account the spatial alignment of the anisometric particles with respect to the strut axis attained during foaming. The model relates nonlinearity in foam deformation to nonlinearly elastic response of foam strut material. Reasonable agreement of model prediction with experimental data is demonstrated for tensile deformation of neat and clay nanoplatelet-filled highly porous rigid polyurethane foams. For neat foams, this suggests that the axial stress-strain response of a foam strut is close to that of a monolithic polymer specimen of a macroscopic size. For composite foams, alignment of the clay platelets with the axial direction of foam struts increased the predicted foam stiffness compared with that for random orientation of platelets. This effect became more pronounced at greater clay loadings. The theoretical deformation diagrams of composite foams where close to or above the experimental ones. As complete exfoliation of the clay particles was assumed in the model, the predicted stiffness of foams being higher than the experimental is likely to stem from only partial exfoliation of clay achieved during manufacture of filled foams.