Finite strain elasticity based cohesive zone model for mechanoluminescent composite interface: I. Stiffness of the undamaged interface

An EML–PDMS composite coupon, made as mentioned above, is flash frozen by immersing in liquid N₂ and then fractured into two pieces by bending. Scanning electron micrographs (SEM) of the fracture planes are captured (figure S2(available online at stacks.iop.org/SMS/30/015016/mmedia)). ImageJ software is used to calculate the total pixel area occupied by all particles in a region of interest [18]. The ratio of particle area to the total area of the region is computed as the volume fraction $\left( {f = 0.4733} \right)$ of the composite. Particle size distribution is estimated from separate SEM images of GG45 and GG13 particles (figure S3). The major and minor diameters (d₁, d₂) of the particles in view are measured using ImageJ software. The particle diameter is calculated as $\left( {\left( {{d_1} + {d_2}} \right)/2 + \sqrt {{d_1}{d_2}} } \right)/2$. The resulting size distribution for a mixture of GG45 and GG13 particles in equal parts by weight is shown in figure S4.

Uniaxial tensile test is performed on PDMS elastomer and EML–PDMS composite specimens following ASTM-D412-16 standard for mechanical characterization [19]. Specimen dimensions recommended by the standard are scaled down by a factor of three to fit within the extension limits of current experimental setup. Rectangular specimens measuring 40 mm, 2 mm and 1 mm in length, width and thickness are used for tensile testing. Despite scaling down specimen dimensions, tests were performed at strain-rates close to 0.25 s⁻¹ as specified in the ASTM standard.

Tensile testing is performed on a specially designed setup assembled on an optical microscope stage as shown in figure 2. X-directional rack and pinion drive of the microscope stage (Meiji Techno MA917R) is motorized by a servo motor. One end of the coupon is held by a grip connected to the motorized drive while other end is held by stationary grip connected to a MLP10 load cell from Transducer Techniques in series which measures the axial force. The objective, eye-piece and top–down LED housing of the microscope are removed for the tensile tests as shown in figure 2(a). This is done to accommodate a Windows Lifecam camera, which is placed vertically above the specimen with its plane of view parallel to the plane of loading. Applied strain is measured from videos of the specimens during the tensile test captured by the camera.

Zoom In Zoom Out Reset image size

Red, green and blue (RGB) values-based edge detection is employed to compute longitudinal and lateral strains in every captured frame. The EML–PDMS composite specimens appear white in the grayscale footage, with RGB values above 200. Two markings, parallel to y-axis and along the edges of the grips, are made on the specimen with a fine-tip black Sharpie. The black markings have RGB values below 50 providing high contrast and easy detection against the white specimens. Transparent PDMS specimens required modifications for RGB-based edge detection. Two additional markings along the specimen edges were made and a white background is used to provide high contrast of markings for easy detection.

The average number of pixels within the two markings ($p{x_l}$) along the loading direction are measured in every frame for accurate longitudinal strain measurement. Similarly, the background is painted dark for maximum contrast with the specimen aiding accurate edge definition and width pixel ($p{x_w}$) computation for lateral strain measurement. Each pixel measures 0.055 mm. The load cell measurements are time-synced with the video to map nominal stress to applied nominal strain. The frame right before application of longitudinal strain is picked as the reference frame. $px_l^{ref}$ and $px_w^{ref}$ are measured as the length and width from reference frame. Nominal strains at n^th frame are computed as ${\varepsilon _x} = px_l^n/px_l^{ref}$, ${\varepsilon _y} = px_w^n/px_w^{ref}$, and the Poisson function is calculated as $\nu = - {\varepsilon _y}/{\varepsilon _x}$.

All composite specimens are strained till failure and all PDMS specimens are strained to large strains (>0.3). Small-strain micromechanics models require elastic modulus (E) and Poisson's ratio $\left( \nu \right)$ of both composite and PDMS matrix. E and $\nu$ are obtained from linear fits within 5% strain limits to the engineering stress–strain curves and the lateral strain to longitudinal strain curves respectively. Towards accurate modeling in finite element analysis (FEA), the non-linear behavior of PDMS is fit with a hyperelastic model using ABAQUS FEA software [20]. Stress relaxation test is also performed on PDMS specimens to obtain its viscoelastic Prony series, which is utilized in FEA to account for rate-dependent behavior. The second part of this paper will take the nonlinear behavior of the composite into consideration towards estimation of CZM damage characteristics. The composite will be modeled as a finite-strain material by utilizing the engineering stress–strain data obtained here.

A 2D square RVE measuring 150 $\mu {\text{m}}$ in length and width with 39 circular particles of varying diameter representing the size distribution is built in ABAQUS FEA software as shown in figure S5. The matrix is shown in blue has 7695 elements and particles in gray have 6463 elements. The number of particles was selected to achieve the volume fraction of the composite $\left( f \right)$. A MATLAB code picks the centers of the particles within the RVE at random while satisfying constraints on minimum particle–particle and particle–boundary distances. The constraints avoid meshing difficulties. The first constraint also avoided the need to define particle–particle interaction properties. Top and bottom edges of the RVE were defined with periodic boundary conditions to ensure spatial repeatability.

All particles are defined as linear elastic material with Young's modulus and Poisson ratio of cubic-ZnS from published datasheets [21]. The particle–matrix interfaces are defined as cohesive surfaces with estimated stiffness $k_\sigma ^{int}$ and no damage. Uniaxial tensile test data of pure PDMS specimens along with the estimated Poisson's ratio is fed to ABAQUS FEA to estimate the hyperelastic fit. The second order Ogden model was selected as it had minimum fit error and was stable for all strains.

The stress–strain response of PDMS matrix and EML–PMDS composite both show rate-dependency and the elastic modulus depends on the rate of applied strain. Hence, the rate-dependent viscoelastic behavior of the matrix is taken into consideration in the model for accurate representation and comparison of simulation results with experimental data. Towards this, Prony series fit of the PDMS modulus relaxation is utilized to describe matrix's viscoelastic behavior. Further, the RVE is strained at the same rate as the experimental tensile test, by utilizing the time-series data from experiments.

Monotonic strain loading applied in the x-direction as displacement boundary condition on the right edge while the left edge is constrained in the x-direction. Monotonic strain amplitude applied on the RVE is obtained from the uniaxial tensile test experiments. Two variants of the model are implemented, one with plane strain and another with plane stress formulation mesh elements (CPE3/4R and CPS3/4R, respectively), with all other the material properties, constraints and boundary conditions remaining identical. The two formulations are used to compare results from the 2D model with experimental data.

Tensile testing of pure PDMS and EML–PDMS composite specimens were performed to estimate elastic moduli and Poisson's ratios. Elastic modulus $\left( {\bar E} \right)$ is computed as the slope of a linear fit to experimental stress–strain curves up to 5% strain limit. Due to limitation on camera framerate, experimental data is captured roughly every 0.8% strain for a tensile loading performed at ASTM standard specified strain-rate. This results in only 5–7 data points within the 5% strain limit per loading curve. Hence, to improve accuracy of computed modulus and Poisson's ratio values, 23 pure PDMS specimens and 10 EML–PDMS composite specimens were tested. Figures 3(a) and (e) show the stress–strain curves for PDMS and EML–PDMS specimens respectively. Figures 3(b) and (f) show the variation of lateral strain to longitudinal strain for PDMS and EML–PDMS specimens respectively.

Zoom In Zoom Out Reset image size

Slope of the linear fits within 5% strain gives the elastic modulus $\left( {\bar E} \right)$ and Poisson's ratio $\left( {\bar \nu } \right)$ of the composite to be 6.93 MPa and 0.348 respectively. Figures 3(c) and (d) give the elastic modulus $\left( {{E^m}} \right)$ and Poisson's ratio $\left( {{\nu ^m}} \right)$ of the PDMS matrix to be 1.808 MPa and 0.423, respectively. Figure S1 shows the stress relaxation curves obtained for plain PDMS specimens along with the Prony series fit for modeling modulus relaxation.

Mori–Tanaka based model detailed by Tan et al [14] is utilized to estimate the undamaged stiffness of the particle–matrix interface in the normal direction. The model links the undamaged stiffness of the interface to the composite bulk modulus considering the mechanical properties of the matrix and particle as well as the volume fraction and size distribution.

$$\begin{equation*}\bar K = {\left[ {\frac{1}{{{K^m}}} + \frac{{9\left( {1 - {\nu ^m}} \right)}}{{2{E^m}}}\frac{{f - \mathop \sum \nolimits_N {f_N}{\alpha _N}}}{{1 - f + \mathop \sum \nolimits_N {f_N}{\alpha _N}}}} \right]^{ - 1}}\end{equation*}$$

$$\begin{equation}\hspace{-20pt} where,\,{\alpha _N} = \,\frac{{3\left( {1 - {\nu ^m}} \right)}}{{2{E^m}\left( {\frac{1}{{k_\sigma ^n{a_N}}} + \frac{1}{{3{K^p}}} + \frac{1}{{4{\mu ^m}}}} \right)}}\,\,\end{equation} \tag{ 3 }$$

where E is the elastic modulus, K is the bulk modulus, $\mu$ is the shear modulus, $\nu$ is the Poisson's ratio and f is the volume fraction of particles in the composite. ${a_N}$ is the radius and ${f_N}$ is volume fraction of the of N^th size of particles.$\,{\alpha _N}$ is defined as the ratio of average stress in particle of N^th size to the average stress in matrix. $k_\sigma ^n$ is the undamaged stiffness of the particle–matrix interface in the normal direction. Superscripts m and p refer to matrix and particles respectively and the bar accent indicates composite properties. ${E^m},\,{\nu ^m},\,\bar E$ and $\bar \nu$ are obtained from figure 3, while ${E^p},\,{\nu ^p}$ are obtained from datasheets [21]. Bulk modulus and shear modulus are calculated as $K = \,E/\left[ {3\left( {1 - 2\nu } \right)} \right]$ and $\mu = \,E/\left[ {2\left( {1 + \nu } \right)} \right]$ respectively. All the computed parameters are tabulated in table 1. For these parameters, through equation (3), the undamaged interface stiffness in the normal direction $k_\sigma ^n\,$ is computed to be 20.55 MPa µm⁻¹. From equation (3), ${\alpha _{N,avg}}$ and $\mathop \sum \limits_N {f_N}{\alpha _N}$ corresponding to $k_\sigma ^n = 20.55{\text{\MPa}}\;\mu {{\text{m}}^{ - {\text{1}}}}$ are 1.1995 and 0.5677, respectively.

Table 1. Material properties used for micromechanical models.

Parameter	Unit	Composite	Matrix	Particles
Young's modulus E	MPa	6.93	1.808	98 100
Poisson's Ratio ν	—	0.348	0.423	0.4
Bulk modulus K	MPa	7.5987	3.9134	163 500
Shear modulus $\mu$	MPa	2.5705	0.635	35 036
Density ρ	G cc−1	2.0073	0.97	4.0856

Energy equivalent inhomogeneity model proposed by Nazarenko et al [15] is utilized to obtain the shear or tangential stiffness of the particle–matrix interface ($k_\sigma ^t)$. Specifically, the spring layer model for the interface of spherical particles is utilized. The bulk and shear moduli of an equivalent inhomogeneity incorporating particle and interfacial spring stiffness are given by the following equations for a zero-thickness interface.

$$\begin{equation}{K^{eq}} = \frac{{a\,k_\sigma ^n\,{K^p}}}{{a\,k_\sigma ^n + 3{K^p}}}\,\,;\,\,{\mu ^{eq}} = \frac{{a\,{\mu ^p}\,\left[ {2\,k_\sigma ^n + 3\,k_\sigma ^t} \right]}}{{\left[ {2a\,k_\sigma ^n + 3a\,k_\sigma ^t + 10{\mu ^p}} \right]}}\,\,\,\end{equation} \tag{ 4 }$$

Using method of conditional moments, Nazarenko et al [15] arrive at the following scalar expressions for effective bulk and shear moduli of the composite.

$$\begin{equation}\bar K = f{K^{eq}} + \left( {1 - f} \right){K^m} - \frac{{f\left( {1 - f} \right){{\left[ {{K^{eq}} - {K^m}} \right]}^2}}}{{f{K^m} + \left( {1 - f} \right){K^{eq}} + \frac{4}{3}{{\tilde \mu }^c}}}\,\,\end{equation} \tag{ 5 }$$

$$\begin{equation}\bar \mu = f{\mu ^{eq}} + \left( {1 - f} \right){\mu ^m} + \frac{{4f\left( {1 - f} \right)\,\widetilde {b\,}{{\left[ {{\mu ^{eq}} - {\mu ^m}} \right]}^2}}}{{1 - 4\,\widetilde {b\,}\left[ {f{\mu ^m} + \left( {1 - f} \right){\mu ^{eq}} - {{\tilde \mu }^c}} \right]}}\,\,\end{equation} \tag{ 6 }$$

where

$$\begin{equation*}\widetilde {b\,} = \frac{{ - 3\left( {{{\tilde K}^c} + 2{{\tilde \mu }^c}} \right)}}{{10{{\tilde \mu }^c}\left( {3{{\tilde K}^c} + 4{{\tilde \mu }^c}} \right)}}\end{equation*}$$

and ${\tilde K^c},\,{\tilde \mu ^c}$ are as follows,

$$\begin{equation}{\tilde K^c} = \,\left\{ {\begin{array}{*{20}{c}} {f{K^{eq}} + \left( {1 - f} \right){K^m}\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,}&{}&{if\,{K^{eq}} \leqslant {K^m}} \\ {{{\left[ {f{{\left( {{K^{eq}}} \right)}^{ - 1}} + \left( {1 - f} \right){{\left( {{K^m}} \right)}^{ - 1}}} \right]}^{ - 1}}}&{}&{if\,{K^{eq}} \geqslant {K^m}} \end{array}} \right.\,\,\,\,\end{equation} \tag{ 7 }$$

$$\begin{equation}{\tilde \mu ^c} = \,\left\{ {\begin{array}{*{20}{c}} {f{\mu ^{eq}} + \left( {1 - f} \right){\mu ^m}\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,}&{}&{if\,{\mu ^{eq}} \leqslant {\mu ^m}} \\ {{{\left[ {f{{\left( {{\mu ^{eq}}} \right)}^{ - 1}} + \left( {1 - f} \right){{\left( {{\mu ^m}} \right)}^{ - 1}}} \right]}^{ - 1}}}&{}&{if\,{\mu ^{eq}} \geqslant {\mu ^m}} \end{array}} \right.\,\,\,\,\end{equation} \tag{ 8 }$$

Utilizing the parameters in table 1 and interface stiffness in the normal direction $(k_\sigma ^n)$ computed from equation (3), we can compute $k_\sigma ^t$ that gives $\bar \mu$ from equation (6) equal to the experimentally measured shear modulus of the composite. Doing so, we obtain $k_\sigma ^t$ to be 88.18 MPa $\mu$m⁻¹. The corresponding composite bulk modulus predicted by equation (5) is 8.13 MPa which when compared to the experimentally measured composite bulk modulus of 7.5987 MPa corresponds to a 7% variation.

It is important to report that the interface stiffness values predicted by these models vary over large ranges as they are extremely sensitive to the Poisson's ratio measurements of the composite and the matrix. Variations in the order of 10⁻² in Poisson's ratio can result in 1 order of magnitude variation in interface stiffness values. Hence, considerable number of experimental data sets were utilized to estimate Poisson's ratio of the EML composite and PDMS matrix.

However, an unavoidable consequence of experimental variability in Poisson's ratio measurements is the result that $k_\sigma ^t\,$is about four times larger than $k_\sigma ^n$. Though the interface can be anisotropic, it is unlikely for the tangential stiffness to be higher than normal stiffness. Hashin [16] predicted that the normal stiffness will be substantially larger than the tangential stiffness, at least for the case of a thin interface (non-zero thickness). This result is for a single spherical inclusion in an infinite matrix with imperfect interfacial bonding which is a simpler problem than Nazarenko et al [15] as it does not take into consideration a composite with a distribution of particles with a certain volume fraction. But adopting Hashin's model [16] ensures the EML–PDMS interface is defined by parameters consistent with isotropic behavior. The tangential stiffness $k_\sigma ^t\,$ is derived by Hashin [16] to be a function of $k_\sigma ^n$ as follows.

$$\begin{equation}\frac{{k_\sigma ^n}}{{k_\sigma ^t}} = \frac{{{K^{int}}}}{{{\mu ^{int}}}} + \frac{4}{3}\,\end{equation} \tag{ 9 }$$

The bulk and shear modulus of the interface are difficult to measure or predict. However, assuming the interface to be isotropic, we can estimate the ratio of the moduli in terms of Poisson's ratio. For isotropic materials, $K = \,E/\left[ {3\left( {1 - 2\nu } \right)} \right]$ and $\mu = \,E/\left[ {2\left( {1 + \nu } \right)} \right]$. Hence,

$$\begin{equation}\frac{{{K^{int}}}}{{{\mu ^{int}}}} = \frac{{2\left( {1 + {\nu ^{int}}} \right)}}{{3\left( {1 - 2{\nu ^{int}}} \right)}}\,\,\,\end{equation} \tag{ 10 }$$

Estimating or measuring the Poisson's ratio of the interface ${\nu ^{int}}$ is also not feasible, but we can assume that it will be similar to the matrix and the particle Poisson's ratios. We therefore assume ${\nu ^{int}} = \left( {{\nu ^p} + {\nu ^m}} \right)/2$. From table 1, ${\nu ^{int}} = \left( {0.4 + 0.423} \right)/2\, = \,0.4115$. Therefore, from equation (10), ${K^{int}}/{\mu ^{int}} = 5.316$ and $k_\sigma ^n/k_\sigma ^t = 6.65$ and for $k_\sigma ^n = 20.55{\text{ MPa}}\;\mu {{\text{m}}^{ - 1}}$, we get $k_\sigma ^t = 3.09{\text{ MPa}}\;\mu {{\text{m}}^{ - 1}}$. The two stiffness parameters defining the interface CZM ($k_\sigma ^n\,and\,k_\sigma ^t)$ are fed into ABAQUS finite element software to define the cohesive behavior of the particle–matrix interface. Damage initiation and stabilization are not defined for the interface.

Monotonic strain loading is applied on the RVE in x-direction. Engineering stress is computed as the sum of reaction forces on the right edge nodes divided by the length of the right edge. Engineering stress strain curves obtained from plane strain and plane stress formulations are plotted along with the experimental tensile test data in figure 4.

Zoom In Zoom Out Reset image size

It can be seen that all the experimental stress strain curves lie within the plane strain and plane stress curves up to 5% engineering strain. Further, it is seen that the plane strain formulation matches the experimental data the best. It is important to note that the interface stiffnesses ($k_\sigma ^{int})$ estimated from small strain assumption based theoretical models are able to characterize tensile behavior of the EML composites in the relatively small strain regime. Figure 5 shows the lateral strain variation with longitudinal strain. There is a reasonable correlation between experimental and numerical results here as well. Despite nonlinear behavior of PDMS matrix material, defining the undamaged interface spring by linear elasticity seems to work reasonably well as long as overall strain remains relatively small.

Zoom In Zoom Out Reset image size