A multiscale material model for heterogeneous liquid droplets in solid soft composites

Abstract

Liquid droplets in solid soft composites have been attracting increasing attention in biological applications. In contrary with conventional composites, which are made of solid elastic inclusions, available material models for composites including liquid droplets are for highly idealized configurations and do not include all material real parameters. They are also all deterministic and do not address the uncertainties arising from droplet radius, volume fraction, dispersion and agglomeration. This research revisits the available models for liquid droplets in solid soft composites and presents a multiscale computational material model to determine their elastic moduli, considering nearly all relevant uncertainties and heterogeneities at different length scales. The effects of surface tension at droplets interface, their volume fraction, size, size polydispersity and agglomeration on elastic modulus, are considered. Different micromechanical material models are incorporated into the presented computational framework. The results clearly indicate both softening and stiffening effects of liquid droplets and show that the model can precisely predict the effective properties of liquid droplets in solid soft composites.

Graphical abstract

Keywords

liquid in solid / soft composite / computational modeling / multiscale model / heterogeneity

1 Introduction

Composite materials, as combinations of at least two different materials to create new materials with enhanced properties, have been widely investigated in recent decades. High strength and fracture toughness, improved fatigue resistance, light weight, higher erosion resistance, corrosion resistance and other advanced characteristics, make them suitable choices for challenges in aerospace, structural, military, and transportation industries [1–3].

Fibrous and particulate are two basic types of composite materials. The latter consists of tiny particles or inclusions embedded in a matrix, having simple and low-cost fabrication procedures, but comparable mechanical performance. They show growing application as biological materials, such as composites for dental fillings (with high moldability in the uncured condition and promising mechanical and physical characteristics once solidified) and hydrogel-based composites including drug release vehicles.

Natural biological tissues could be considered to be a liquid-in-solid soft particulate composite-like materials. Artificial polymeric types are also well known and acknowledged for their potential in biomechanics research and bio-related applications. For example, bone and soft tissues with regeneration and repair characteristics show promising applications. Soft composites can also provide platforms for stretchable electronic systems that intimately integrate with the human body. Brown et al. [4] presented a liquid droplet-solid epoxy matrix composite as self-healing biological bone structure. Dong et al. [5] presented PDMS-PVDF composites with self-healing and self-stiffening properties when subjected to external loads. Agrawal et al. [6] analyzed dynamic response of a liquid crystals-polymeric composite as a biocompatible-tissue. Thai et al. [7] proposed an isogeometric approach for flexoelectricity in soft dielectric materials, accounting for Maxwell stresses at the interface. A 3D isogeometric flexoelectricity is presented in Ref. [8].

Several theoretical (including analytical and numerical) and experimental works have been done to obtain the elastic modulus of a solid in solid particulate composites [9–14]. Among them, Eshelby’s inclusion theory, which is based on classical elasticity and is for an ellipsoidal inclusion embedded in an infinite elastic matrix, has been widely used [15]. However, direct use of Eshelby’s inclusion theory fails to precisely predict the elastic modulus of liquid-in-solid composites, because it cannot capture inclusion size or boundary effects [16]. In other words, liquid and solid inclusions are distinctly different in terms of shear and Young’s modulus, which are zero for liquid phase, and surface tension effect at the liquid interface. In fact, the surface energy mismatch at the liquid-solid matrix interface may noticeably impact the elastic modulus of the bulk composite at macroscale.

Analytical homogenization techniques [17–19] of complex shape inclusions are difficult to carry out. Therefore, FEM is employed to get over the inclusion’s geometry and distribution restrictions [20–23]. This technique includes applying special load cases and imposing suitable periodic boundary conditions to the unit cell possessing the characteristic of the whole composite structure. An important advantage of FEM is that there is no restriction on size of the inclusions, their spatial distribution, the number of components or their properties.

There are some theoretical works that predict effective mechanical properties of liquid in solid soft composites. Style et al. [24] has successfully extended Eshelby’s theory into the field of liquid-in-solid soft composite. However, although the surface effect is incorporated, the model does not address the polydispersity of the liquid droplets in the solid. Moreover, the effect of chain-like liquid inclusions or the interactions between close inclusions remains unexplored. A comprehensive review of solid-liquid composites in terms of their fabrication techniques, mechanical properties and applications is presented in [25].

So far, the developed material models for liquid in solid soft composites are for highly idealized configurations, in particular monodisperse type composites which consist of identical liquid droplets. It is believed that polydispersity of particle size may affect the bulk mechanical properties [26,27]. Moreover, it is necessary to take microstructural imperfections into account, since they may lead to approximation errors. Imperfections that influence the mechanical properties are for instance: deviations in droplets shape and size, nonuniformly distributed droplets, local changes in volume fraction and droplet-matrix interface interactions. For more accurate prediction of the bulk properties of these composites, there is a need for models to incorporate the effects of the most influential imperfections. Consequently, this requires new methods to quantify these imperfections.

This research work presents a computational framework to obtain the effective elastic modulus of a typical polydisperse liquid-in-solid soft composite, addressing actual conditions and available uncertainties. The effects of surface tension at droplets’ interfaces, their volume fraction, size, size polydispersity and agglomeration on the elastic modulus, are considered. It also covers uncertainties at different length scales (from micro-, meso- to macroscale) via a random approach. It will be the first framework bridging models from micro- to macroscale for liquid in solid soft composites.

2 Multiscale methodology

Motivated by Refs. [28,29], this research presents a multiscale approach to obtain the effective Young’s modulus of the heterogenous liquid in solid soft composite. As shown in Fig. 1, it utilizes a bottom–up approach, in which the material features at each length scale are assembled from the smaller length scales and then proceed to form the macroscopic properties of the composite material, hierarchically (from micro-, meso- to macroscale).

The presented model covers uncertainties at different length scales. At microscale the liquid droplet is substituted with the equivalent elastic solid inclusion. At macroscale the material is tessellated into constitutive building blocks. A single building block forms the Representative Volume Element (RVE) at mesoscale. The specific amount of the inclusions, determined by the total volume fraction of inclusions in the composite, is randomly distributed over these building blocks. Local agglomeration and polydispersity are addressed at mesoscale. Assuming macroscopic homogeneity of the composite, the effective Young’s modulus of the composite is eventually estimated by averaging the individual values corresponding to each block. A certain number of realizations is done to achieve convergence in results. A schematic representation of this computational multiscale material model is presented in Fig. 2.

2.1 Microscale

Inspired by Refs. [30,31], instead of the liquid droplet, the equivalent solid elastic inclusion at microscale is utilized. In fact, the liquid droplet, with internal pressure and surface tension, is substituted with the equivalent elastic solid one, perfectly bounded at the interface with the surrounding matrix. This technique makes it possible to use the micromechanics rules and equations that are applicable for solid particles.

The Young-Laplace equation [24,31] to model the discontinuity of the traction vector t at the interface is:

(1)\sigma \cdot \mathit{n}=-p\mathit{n}+\gamma \kappa \mathit{n},

where \mathit{n} is the normal to the deformed droplet surface, \sigma \cdot \mathit{n} is the normal stress on the solid side, p is the pressure in the droplet, \kappa is the sum of the local principal curvatures and \gamma is the constant, isotropic strain-independent surface tension. Using Eq. (1), the liquid droplets are found to be then equivalent to solid elastic inclusions of Young’s modulus E_{\text{inc}}, obtained by the following formula [32]:

(2)E_{\text{inc}}=\frac{4}{1+{\upsilon }_{\text{polymer}}+\left(5/3\right)\left(R/L\right)}{E}_{\text{polymer}},

where {\upsilon }_{\text{polymer}} and E_{\text{polymer}} are the Poisson’s ratio and the Young’s modulus of the polymeric host, R is the radius of the liquid droplet and L\equiv {\displaystyle \frac{\gamma }{E_{\text{polymer}}}} is the elastocapillary length scale.

Figure 3 represents the methodology which is implemented at microscale. Now it is possible to follow up the existing rules valid for the particulate-filled polymeric composites. Plenty of theoretical studies have been already performed on this type of composites to obtain their effective mechanical properties, in particular their elastic modulus. Analytical methods have been developed aiming to increase the accuracy of the results or decrease the computational cost. Among these methods Voigt [9], Reuss [10], Hashin-Shtrikman Bounds [11], variational bounding techniques [12], Self-consistent [13] and Mori-Tanaka [14] methods can be counted. Exhibiting isotropic behavior due to random dispersion of equivalent elastic solid inclusions, Mancarella [32], Voigt, Rule of Mixtures (ROM) and the semi-empirical Halpin−Tsai (H−T) [33] equations are used. For the case of spherical particles, the H−T model gives a good estimation of the Young’s modulus of the composite as reported in literatures [34]. Detail derivations for each model are presented in Section 2.2.

2.2 Mesoscale

The agglomeration, polydispersity and volume fraction of liquid droplets are characterized at this scale. The embedded equivalent elastic solid inclusions with different diameters are randomly dispersed in the matrix. They can be either concentrated in local spherical aggregates or uniformly dispersed. In fact, at mesoscale, the RVE consists of randomly dispersed equivalent solid elastic inclusions in the matrix, which form a percolated network of the dispersed inclusions. Some of these percolations are not continuous throughout the material, but are highly localized (illustrated by continuous circles in Fig. 4), so it is considered agglomerated. In other words, any inclusion could be either inside or outside of the spherical aggregates. The RVE at mesoscale is shown in Fig. 4.

The Shi et al. [35] agglomerate model for soft composites is adopted here. According to this model, the total volume V_{\text{inc}}^{\text{total}} of inclusions in the RVE is given by:

(3)V_{\text{inc}}^{\text{total}}=V_{\text{inc}}^{\text{agg}}+V_{\text{inc}}^{\mathrm{m}},

where V_{\text{inc}}^{\text{agg}}and V_{\text{inc}}^{\mathrm{m}} denote the volumes of inclusions dispersed in the agglomerated regions and in the percolated matrix, respectively. Moreover, the following \xi and \zeta parameters are defined to describe the agglomeration of inclusions:

(4)\xi =\frac{V_{\text{agg}}}{V}\mathrm{a}\mathrm{n}\mathrm{d}\zeta =\frac{V_{\text{inc}}^{\text{agg}}}{V_{\text{inc}}^{\text{total}}},

where V_{\text{agg}} is the volume of the spherical agglomerations in the RVE, with the total volume of V. When \xi =1, inclusions are uniformly dispersed in the matrix. The parameter \zeta denotes the volume ratio of inclusions that are dispersed in agglomerations and the total volume of the inclusions. When \zeta =1, all the inclusions are located in the sphere areas. \xi =\zeta, describes the uniform dispersion of all inclusions.

Defining the average volume fraction of inclusions, C_{\text{inc}}, in the composite as:

(5)C_{\text{inc}}=\frac{V_{\text{inc}}^{\text{total}}}{V},

one can express the volume fractions of inclusions in the agglomerations and in the matrix, as:

(6)\frac{V_{\text{inc}}^{\text{agg}}}{V_{\text{agg}}}=\frac{C_{\text{inc}}\zeta }{\xi }\mathrm{a}\mathrm{n}\mathrm{d}\frac{V_{\text{inc}}^{\mathrm{m}}}{V-V_{\text{agg}}}=\frac{C_{\text{inc}}\left(1-\zeta \right)}{1-\xi }.

Using the Voigt model, the effective elastic stiffness of the agglomerations E_{\text{in}}, and the matrix E_{\text{out}}, can be estimated:

\begin{array}{rl}{E}_{\text{matrix}}=& E_{\text{out}}=\frac{3}{8}\left\{\frac{C_{\text{inc}}\left(1-\zeta \right)}{1-\xi }{E}_{\text{inc}}+\left[1-\frac{C_{\text{inc}}\left(1-\zeta \right)}{1-\xi }\right]E_{\text{polymer}}\right\}\\ & +\frac{5}{8}\left\{\frac{\left(1-\xi \right)E_{\text{inc}}{E}_{\text{polymer}}}{\left[\left(1-\xi \right)-C_{\text{inc}}\left(1-\zeta \right)\right]E_{\text{inc}}+C_{\text{inc}}\left(1-\zeta \right)E_{\text{polymer}}}\right\}\end{array}

and

(7)\begin{array}{rl}{E}_{\text{agg}}=& E_{\text{in}}=\frac{3}{8\xi }\left[C_{\text{inc}}\zeta E_{\text{inc}}+\left(\xi -C_{\text{inc}}\zeta \right)E_{\text{polymer}}\right]\\ & +\frac{5}{8}\frac{\xi E_{\text{polymer}}{E}_{\text{inc}}}{\left(\xi -C_{\text{inc}}\zeta \right)E_{\text{inc}}+C_{\text{inc}}\zeta E_{\text{polymer}}},\end{array}

where both the polymer and the inclusions are considered to be isotropic, with Young’s modulus E_{\text{polymer}} and E_{\text{inc}}, respectively.

The H−T equations take the following forms when used inside a single spherical aggregate to obtain its elastic modulus E_{\text{agg}}, or in the percolated network of the dispersed inclusions to obtain its elastic modulus E_{\text{matrix}}:

(8)\begin{array}{rl}& E_{\text{matrix}}=\frac{1+\lambda \eta V_{\text{inc}}^{\mathrm{m}}}{1-\eta V_{\text{inc}}^{\mathrm{m}}}{E}_{\text{polymer}},\\ & E_{\text{agg}}=\frac{1+\lambda \eta V_{\text{inc}}^{\text{agg}}}{1-\eta V_{\text{inc}}^{\text{agg}}}{E}_{\text{polymer}},\end{array}

with \eta ={\displaystyle \frac{\left(E_{\text{inc}}/E_{\text{polymer}}\right)-1}{\left(E_{\text{inc}}/E_{\text{polymer}}\right)+\lambda }}, \lambda =2l/d, where l and d are the length and diameter of the particle.

A simpler model based on twice-use of ROM is also presented. First, the ROM is applied inside a single spherical aggregate to obtain its elastic modulus E_{\text{agg}}, using:

(9)E_{\text{agg}}=V_{\text{inc}}^{\text{agg}}{E}_{\text{inc}}+\left(1-V_{\text{inc}}^{\text{agg}}\right)E_{\text{polymer}},

then, again ROM is applied in the percolated network of the dispersed inclusions to obtain its elastic modulus E_{\text{matrix}}, using:

(10)V_{\text{inc}}^{\mathrm{m}}=\frac{E_{\text{matrix}}-E_{\text{polymer}}}{E_{\text{inc}}-E_{\text{polymer}}},

which can be written in the form of the well-known ROM equation:

(11)\begin{array}{rl}{E}_{\text{matrix}}& =C_{\text{inc}}\left(1-\zeta \right)E_{\text{inc}}+\zeta E_{\text{polymer}}\\ & =V_{\text{inc}}^{\mathrm{m}}{E}_{\text{inc}}+\left(1-V_{\text{inc}}^{\mathrm{m}}\right)E_{\text{polymer}}.\end{array}

For all Voigt, H−T and ROM models, the following formula is used to obtain the effective Young’s modulus of the composite in each constitutive block:

(12)E_{\text{composite}}=\frac{1}{{\displaystyle \frac{1-V_{\text{agg}}}{E_{\text{matrix}}}}+{\displaystyle \frac{V_{\text{agg}}}{E_{\text{agg}}}}}.

For Mancarella model all details are presented in Ref. [32]. Schematic representation of the modeling at mesoscale is presented in Fig. 5.

2.3 Macroscale

Generally, the size, shape, inclusion-matrix interfacial adhesion, and the distribution of the inclusions affect the macroscopic behavior of the composites. The polydispersity in the composite is captured by assuming that the material region is composed of several building blocks, where parameters such as droplet radius, volume fraction, dispersion and agglomeration, can vary locally. The overall isotropic and linearly elastic response of the composite is obtained by averaging the corresponding properties of all constitutive blocks, which come from the mesoscale (Fig. 2).

3 Numerical examples

Glycerol droplets and Silicone gels (of two different stiffness E=3 and 100\text{kPa}) are considered as liquid droplets and incompressible ( \upsilon =0.5) matrix materials, respectively. Other input data and settings are: \gamma =0.0036\mathrm{N}/\mathrm{m}, R=1\times 10^{-6}\mathrm{m}, {\delta }_{\mathrm{R}}=5\mathrm{\%}R, and {\delta }_{\mathrm{C}}=50\mathrm{\%}{C}_{\text{inc}}, where {\delta }_{\mathrm{R}} and {\delta }_{\mathrm{C}} represent the standard deviation of the normal distribution of random variables. A 10\text{mm}\times 10\text{mm} square region is assumed as RVE at macroscale which is tessellated by 250\times 250 square blocks. A total number of 500 realization satisfies the convergence criterion (less than 0.5\mathrm{\%}).

Figure 6(a) depicts the stiffening effect of the small liquid droplets, with the average radius of R=1\text{μm}, in the soft polymeric matrix ( E=3\text{kPa}), where the capillarity effect dominates. In fact, this effect typically appears at length scales \lesssim L\equiv \gamma /E. In this case the capillary length is L=1.2\text{μm}. When the size of liquid droplets increases, R=10\text{μm} as shown in Fig. 6(b), the softening behavior happens, because of the bulk elasticity. In Figs. 6(c) and 6(d) with the stiffer matrix ( E=100\text{kPa}), giving L=0.036\text{μm} which is noticeably smaller than both R=1\text{μ}\text{m} and R=10\text{μ}\text{m}, and leads to the softening behavior. In all cases, the stiffening and softening behaviors become stronger with increase in the droplet’s volume fraction. Moreover, for the same data settings and in all cases, the heterogeneous H−T and ROM show the least effect on the Young’s modulus of the composite. The deterministic and heterogeneous Mancarella [32] models give the steepest effect in all the softening and stiffening cases, respectively. The heterogeneous Voigt model is not applicable for the softening regimes and only converges in stiffening cases.

Figure 7 belongs to the soft polymeric matrix with E=3\text{kPa} and depicts the role of droplet size in the stiffening and softening behaviors of the heterogeneous H−T, Mancarella and ROM models. One can see the stiffening effects of the liquid droplets with decreasing their radius (less than the capillary length, L=1.2\text{μ}\text{m}) and increasing their volume fractions. By increasing the size of droplets and their volume fractions, the softening behavior of the composite happens, as expected from the bulk elasticity, because of replacement of some solid matrix materials with the liquid phase with zero stiffness.

The surface tension \gamma can also vary the value of L. As shown in Fig. 8, by increasing \gamma, the capillary length also increases and for \gamma =0.001 and E=3\text{kPa}, we have L=0.33\text{μm} which is smaller than R=1\text{μ}\text{m}; thus, the softening effect happens. The bigger the increase in \gamma, then the bigger the increase in L, and since R=1\text{μ}\text{m} is constant, the stiffening effect becomes intense.

4 Concluding remarks

A multiscale computational framework, bridging models from microscale to macroscale, for liquid in solid soft composites is presented. It addresses the available uncertainties that result from different causes. It utilizes a bottom–up approach, by which the material features at each length scale are assembled from the smaller length scales, and then proceeding to form the macroscopic properties of the composite material, hierarchically (from micro-, meso- to macroscale). At microscale the liquid droplet is substituted with the equivalent elastic solid inclusion. The agglomeration, polydispersity and volume fraction of liquid droplets are characterized at mesoscale. The effective Young’s modulus of the composite at macroscale is estimated by averaging the values corresponding to the mesoscale. The effect of surface tension at droplets’ interfaces, their volume fraction, size, size polydispersity and agglomeration on the elastic modulus, are captured. The results indicate the stiffening behavior when the size of liquid droplets is smaller than the elastocapillary length, given by the ratio of the surface tension to Young’s modulus of the solid matrix. By increasing the size of droplets, the softening behavior of the composite happens as expected from the bulk elasticity. In all cases, the stiffening and softening behaviors become stronger with increase in droplets’ volume fractions. Conducting experimental investigations, extending the two-phase model to a multiphase model and performing the sensitivity analysis to identify the model key input parameters, which have the most influence on the effective Young’s modulus of the liquid in solid soft composites, present needs for further research.

References

[1] Ghasemi H, Park H S, Rabczuk T. A multi-material level-set based topology optimization of flexoelectric composites. Computer Methods in Applied Mechanics and Engineering, 2018, 332 : 47– 62

[2] Ghasemi H, Kerfriden P, Bordas S P A, Muthu J, Zi G, Rabczuk T. Probabilistic multiconstraints optimization of cooling channels in ceramic matrix composites. Composites. Part B, Engineering, 2015, 81 : 107– 119

[3] Ghasemi H, Brighenti R, Zhuang X, Muthu J, Rabczuk T. Optimization of fiber distribution in fiber reinforced composite by using NURBS functions. Computational Materials Science, 2014, 83 : 463– 473

[4] Brown E, Sottos N R, White S R. Fracture testing of a self-healing polymer composite. Experimental Mechanics, 2002, 42( 4): 372– 379

[5] Dong P, Chipara A C, Loya P, Yang Y, Ge L, Lei S, Li B, Brunetto G, Machado L D, Hong L, Wang Q, Yang B, Guo H, Ringe E, Galvao D S, Vajtai R, Chipara M, Tang M, Lou J, Ajayan P M. Solid–liquid self-adaptive polymeric composite. ACS Applied Materials & Interfaces, 2016, 8( 3): 2142– 2147

[6] Agrawal A, Chipara A C, Shamoo Y, Patra P K, Carey B J, Ajayan P M, Chapman W G, Verduzco R. Dynamic self-stiffening in liquid crystal elastomers. Nature Communications, 2013, 4( 1): 1739–

[7] Thai T Q, Rabczuk T, Zhuang X. A large deformation isogeometric approach for flexoelectricity and soft materials. Computer Methods in Applied Mechanics and Engineering, 2018, 341 : 718– 739

[8] Ghasemi H S.. Park H, Zhuang X, Rabczuk T. Three-dimensional isogeometric analysis of flexoelectricity with MATLAB implementation. Computers, Materials & Continua, 2020, 65( 2): 1157– 1179

[9] Voigt W. On the relationship between the two elastic constants of isotropic bodies. Annals of Physics, 1889, 274(12): 573–587 (in Germen)

[10] Reuss A. Calculation of the yield point of mixed crystals based on the plasticity condition for single crystals. Journal of Applied Mathematics and Mechanics, 1929, 9( 1): 49– 58

[11] Hashin Z, Shtrikman S. A variational approach to the theory of elastic behavior of multiphase materials. Journal of the Mechanics and Physics of Solids, 1963, 11( 2): 127– 140

[12] Paul B. Prediction of elastic constants of multiphase materials. Transactions of the Metallurgical Society of AIME, 1960, 36– 41

[13] Hill R. A self-consistent mechanics of composite materials. Journal of the Mechanics and Physics of Solids, 1965, 13( 4): 213– 222

[14] Mori T, Tanaka K. Average stress in matrix and average elastic energy of materials with mis-fitting inclusions. Acta Metallurgica, 1973, 21( 5): 571– 574

[15] Nasser N, Hori S M. Micromechanics: Overall Properties of Heterogeneous Solids. Amsterdam: Elsevier, 1993

[16] Gao X L, Ma H M. Solution of Eshelby’s inclusion problem with a bounded domain and Eshelby’s tensor for a spherical inclusion in a finite spherical matrix based on a simplified strain gradient elasticity theory. Journal of the Mechanics and Physics of Solids, 2010, 58( 5): 779– 797

[17] Engquist B, Runborg O. Wavelet-based numerical homogenization with applications. In: Barth T J, Chan T, Haimes R, eds. Multiscale and Multiresolution Methods. Berlin: Springer Heidelberg, 2002, 97–148

[18] Brewster M, Beylkin G. A multiresolution strategy for numerical homogenization. Applied and Computational Harmonic Analysis, 1995, 2( 4): 327– 349

[19] Yousefi H, Farjoodi J, Mahmoudzadeh Kani I. Adaptive simulation of wave propagation problems including dislocation sources and random media. Frontiers of Structural and Civil Engineering, 2019, 13( 5): 1054– 1081

[20] Salavati M, Ghasemi H, Rabczuk T. Electromechanical properties of Boron Nitride Nanotube: Atomistic bond potential and equivalent mechanical energy approach. Computational Materials Science, 2018, 149 : 460– 465

[21] Abbès A, Abbès B, Benkabou R, Asroun A. A FEM multiscale homogenization procedure using nanoindentation for high performance concrete. Journal of Applied and Computational Mechanics, 2020, 6( 3): 493– 504

[22] Kamiński M. Homogenization-based finite element analysis of unidirectional composites by classical and multiresolutional techniques. Computer Methods in Applied Mechanics and Engineering, 2005, 194( 18−20): 2147– 2173

[23] Yi S, Xu L, Cheng G, Cai Y. FEM formulation of homogenization method for effective properties of periodic heterogeneous beam and size effect of basic cell in thickness direction. Computers & Structures, 2015, 156 : 1– 11

[24] Style R W, Wettlaufer J S, Dufresne E R. Surface tension and the mechanics of liquid inclusions in compliant solids. Soft Matter, 2015, 11( 4): 672– 679

[25] Style R W, Tutika R, Kim J Y, Bartlett M D. Solid–liquid composites for soft multifunctional materials. Advanced Functional Materials, 2021, 31( 1): 2005804–

[26] Toulemonde P A, Diani J, Gilormini P, Desgardin N. A numerical study of the influence of polydispersity on the behaviour until break of a reinforced hyperelastic material with a cohesive interface. Matériaux & Techniques, 2015, 103( 3): 306–

[27] Zhong X A, Knauss W. Effects of particle interaction and size variation on damage evolution in filled elastomers. Mechanics of Advanced Materials and Structures, 2000, 7( 1): 35– 53

[28] Shokrieh M M, Rafiee R. Stochastic multi-scale modeling of CNT/polymer composites. Computational Materials Science, 2010, 50( 2): 437– 446

[29] Ghasemi H, Rafiee R, Zhuang X, Muthu J, Rabczuk T. Uncertainties propagation in metamodel-based probabilistic optimization of CNT/polymer composite structure using stochastic multi-scale modeling. Computational Materials Science, 2014, 85 : 295– 305

[30] Duan H L, Yi X, Huang Z P, Wang J A. Unified scheme for prediction of effective modulus of multiphase composites with interface effects, Part I: Theoretical framework. Mechanics of Materials, 2007, 39( 1): 81– 93

[31] Style R, Boltyanskiy R, Allen B, Jensen K E, Foote H P, Wettlaufer J S, Dufresne E R. Stiffening solids with liquid inclusions. Nature Physics, 2015, 11( 1): 82– 87

[32] Mancarella F, Style R W, Wettlaufer J S. Surface tension and the Mori–Tanaka theory of non-dilute soft composite solids. Mathematical, Physical and Engineering Sciences, 2016, 472( 2189): 20150853–

[33] Halpin J C, Tsai S W. Effects of Environmental Factors on composite materials. Air Force Materials Lab Wright-Patterson AFB OH, 1969, 67 : 423–

[34] Thorvaldsen T, Johnsen B B, Olsen T, Hansen F K. Investigation of theoretical models for the elastic stiffness of nanoparticle-modified polymer composites. Journal of Nanomaterials, 2015, 281308–

[35] Shi D L, Feng X Q, Huang Y Y, Hwang K C, Gao H. The effect of nanotube waviness and agglomeration on the elastic property of carbon nanotube-reinforced composite. Journal of Engineering Materials and Technology, 2004, 126( 3): 250– 257

RIGHTS & PERMISSIONS

Higher Education Press 2021.