A Numerical Method Charactering the Electromechanical Properties of Particle Reinforced Composite Based on Statistics.

A novel model for a network of polymer chains is proposed considering the distribution of polymer chains inside the composite in this work. Some factors that influence the distribution of polymer chains are quantitatively investigated, such as external surface geometry, internal filler, and local deformation. Furthermore, the Maxwell stress induced by an electric field is characterized by the statistics of local charge density, as the basic analyzing electromechanical properties of materials. In particular, taking the non-uniform distribution of polymer chains into account, the electromechanical properties of two materials-VHB 4910 and CaCu₃Ti₄O12-polydimethylsiloxane (CCTO-PDMS)-are investigated to validate the applicability of the proposed model. The comparison between simulation results and experimental results from existing literature shows that the model was successfully employed to predict the electromechanical properties of polymer composites.

1. Introduction

Structures consisting of soft materials can undergo large deformation under the effects of external electric fields and force fields, and can be used for different applications such as actuators, artificial muscles, and energy harvesters [1,2]. The experimental investigations of planar actuators have shown that the Maxwell stress induced by the external electric field is generally proportional to the voltage value, and in inverse proportion to the thickness [3]. Methods to reduce the driving voltage (generally kV) have been investigated, including reducing thickness (thin membrane), expanding length (pre-stretch), and increasing dielectric permittivity (adding filler). However, the applications of soft materials are limited due to their failure mechanism both for pure mechanical (loss of tension, rupture by stretch and aging) and electromechanical loading (dielectric breakdown, wrinkling/buckling, and loss and instability) [4]. Therefore, there have been many attempts to develop polymer–filler composites or new structures to overcome the drawbacks mentioned above [5,6].

As for composites with different geometries or conductive fillers, the material properties are quite complex and related to the microstructures [7,8]. The morphology examined by scanning electron microscopy (SEM) indicates that around particles or near surfaces, the polymer chains will re-distribute according to the manufacture process and the interaction between the fillers and the matrix [9,10]. Finite element method at microscale, principle of energy conservation at macroscale, and modification of parameters are generally used to solve this problem in these composites [11,12,13].

In recent years, conductive fillers have been added to polymer-based composites to develop high-performance dielectric polymers (high dielectric permittivity and low loss tangent), such as carbon nanotubes [14,15], graphite [16,17], metal particles [18,19], and conductive oxides [20,21,22]. Generally, higher dielectric permittivities can be achieved in composites by adding these particles, however, unexpected loss tangents are also found, especially under high temperature and frequency. Among those particles, CaCu3Ti4O12 (CCTO) shows a high dielectric permittivity which is almost unchanging over a range of temperature and frequency greater than 104 [23,24,25].

The macro-mechanical properties should originate from the microstructure of polymers, including spatial distribution, interfacial reaction, and chemical modification [26,27]. Based on the analysis of microstructure [28], several multi-scale models have been proposed, such as the three-chain model, eight-chain model, and tetrahedral model [29,30,31]. Beyond that, hyperelastic models aiming to capture the limiting stretch of polymer chains are always used to describe the mechanical behavior of soft polymers [32,33]. A micro-mechanical model for finite rubber elasticity based on a non-affine network and topological constraints was proposed by Miehe and Göktepe [34]. The finite rubber viscoelasticity and anisotropic Mullins-type damage is further discussed in their follow-up articles [35,36]. By studying experimental data and the microstructure of carbon black-filled rubbers, the anisotropic Mullins effect is investigated by Dargazany and Itskov [37]. In this work, a damage polymer–filler network is introduced to qualitatively illustrate chain sliding on or debonding from aggregates. The model is further developed in their work: chain length statistics are incorporated into full network rubber [38]; a closed form of the Rayleigh exact distribution function for non-Gaussian chains is considered [39]; a cross link-particle (CP) network is added to the elastic rubber (cross link- cross link; CC) and polymer-filler (particle-paritcle; PP) network [40]. Moreover, several studies in the literature focusing on the electromechanical properties of polymers revealed the relationship between network structure and electromechanical loading/response [41,42].

The foundation of the models mentioned above is the hypothesis of the random distribution of the polymer chain. The objective of this work is to develop a theoretical model to quantify the structure–function relationship by modifying the distribution of the polymer chains. The distribution of the polymer chain is determined by the geometry characteristic and the interaction between matrix and fillers. This work could shed some light on the development of polymer-based composites.

2. Improved Visco-Hyperelastic Model

2.1. Strain Energy Potential and Viscoelasticity

A system with internal structures , (with normal vector ) and outer space is introduced to illustrate the model. At the reference state t = 0, the undeformed stress-free configuration will be maintained. At time t, force field and electric field are applied on the outer space and , respectively. Meanwhile, the system has changed to current configuration . The deformation is described by function , which maps a reference particle X in to its deformed position in . The associated deformation gradient will be denoted by , and identifies its determinant, as shown in Figure 1.

Figure 1

Illustration of the system in reference state t = 0 and current state t.

In arbitrary subregion , the integral of the flux density over a closed surface equals the charge enclosed based on Gauss’s law: where and are the flux density and charge, respectively. The constitutive equations for current density I, electric field E, and charge Q are expressed as: where is the electric mobility, V is the applied voltage, and are the dielectric permittivity for material and free space, respectively. We should note that Equations (2)–(5) are suitable for materials that have relatively stable performances. However, dielectric permittivities of most of soft materials are reported as a function of the deformation and frequency f, which can be expressed as .

A surface charge with density is formed under the effect of voltage V. Thus, the Maxwell stress is induced by q and the electric field [43]. With force Pf, a new shape is formed ; this is beyond the explanation of pure force working on . Similarly, a new surface charge density is formed:

The local electric breakdown will happen when the electric field of two points is larger than the dielectric breakdown strength [44]:

In order to ensure global electric breakdown, Equation (7) is extended to two points on the outer surface A and B : the general electric breakdown of the system happened if a line with two ends A and B across the body are found.

2.2. Distribution of Polymer Chain

In this section, the effective polymer chain length L = nl (L is the total length of a polymer chain; n is the segment number; l is the length of a single segment) is introduced to investigate the distribution of polymer chains. We start by considering the random distribution, which means the polymer chain is inside the system and far from any boundaries (both inside and outside). In this situation, a spherical coordinate system is used for a better description, and the possible rotation angle of a polymer chain around point X can be expressed as: where and are the space rotation angles. Space rotation angle will be restricted when the point is near boundary. As shown in Figure 2, a plane with X is intended to illustrate the restricted behavior of a surface with different radius of curvature (r, , ). Particularly, the space rotation angle of point X can be obtained: where is the distance from point X to a point A on the boundary.

Figure 2

Effect of surface on the distribution of polymer chains on a plane: (a) ; (b) .

The probability of random walking in one dimension without restriction can be expressed as a Gaussian distribution: where x = ml is the walking distance and m is the testing step measured from the origin ().

The random walking in three dimensional space can be described as: n set of rotational angles and have arbitrary value ( and , i = 1, 2, …, n) as shown in Figure 3. In contrast, the values of and under the restricted condition can be expressed by Equation (10). It is interesting to note that internal structures have different effects on the distribution of the polymer chains according to the interaction of different materials and the manufacturing operation: first, with gas-like material in an internal surface , the polymer chain will automatically avoid the region, resulting in a conditional distribution of polymer chains in this region; second, with particle reinforced in an internal surface, the polymer chains show cross-linking characteristic around this region, meaning that several chains will share the same point. The length of polymer chains will be changed due to this reinforcement. Generally, we assume that chains are walking through this region, and segments are cross-linked around the particle; thus, the modified segment number of the chain can be expressed as: where is the internal segment number of the reinforced chain.

Figure 3

The distribution of polymer chains: (a) random distribution; (b) restricted distribution.

In Equation (12), it is obvious that the segment number is larger than that of random distribution. For the convenience of simulation, it is suggested that a parameter k (k > 1) should be used in this situation. The details of the numerical method are shown in Figure 4. According to Figure 4, n set of and (i = 1, 2, …, n) will first be automatically generated. Combining the coordinate of X, the coordinate of the point will be obtained. Following the same procedure, the new coordinate of point will be obtained if the coordinate () meets the condition of constraint as shown in Figure 4. Otherwise, the values of and should be re-calculated according to Equation (10). Repeating the same step, the coordinates of the n points can be obtained.

Figure 4

The numerical method used to calculate the probability distribution.

Considering the restricted behavior, the probability tested from point X to another point is obtained as:

2.3. The Mechanical Model of Soft Materials

The relationship between stress and strain can be summarized as follows [45,46]: where and are the material parameters. Equation (14) is the general expression for the modified Kelvin–Voigt model. Applying the hydrostatic pressure into Equation (14), we can rewrite the stress–strain relationship as: where is the Kronecker delta ( = 1 if i = j, and = 0 otherwise). In Equation (15), only normal stress is considered (shear stress can be ignored) as the main mode; shear stress is inappropriate for soft materials.

Many studies have attempted to separate the time-dependent behavior from mechanical models by introducing a time function, as can be seen in [46,47]. Similarly, the equivalent modulus (a function of time and strain) can be used for the characterization of a viscoelastic material: where denotes the equivalent elastic modulus and is the unit tensor.

When we consider np models connected in parallel, the relationship between stress and strain of the point X integrating probability can be expressed as: where P is the matrix expression of probability as shown in Equation (13).

2.4. Re-Distribution of Charge on the Surface

In this section, two distribution laws of the charges on the surface of a material are investigated. Firstly, the homogeneous distribution which is usually used in the planar actuator is shown in Figure 5a. Three points O, A, and B with distance are chosen to study the distribution law. Two dielectric permittivities and are used to characterize the dielectric property inside and outside the material, respectively. The charges on three points at homogeneous state are the same:

Figure 5

Distribution of charges on the outer surface of material: (a) homogeneous; (b) inhomogeneous.

Secondly, an assumption is made over a small time: the charge cannot move from a point to another point along a smooth surface, as shown in Figure 5b. Taking point O as an example, the components of force along the tangential direction () and the vertical direction () are as follows: where and are the angles, . The deformation is changing due to force and force field, which will result in charge redistribution. It is obvious that the ratio is a function of (or ) in a small area. When Equation (19) is extended to the general situation of the whole material, the following expression can be obtained:

Along with Equations (1)–(6), the charge redistribution behavior can be obtained at any time.

3. The Effect of External Structures on the Electromechanical Properties

3.1. Determination of the Material Parameters

A sample with dimension (, ) and loading in z-direction is taken to illustrate the influence of l and n on the mechanical properties. As shown in Figure 6, the equivalent elastic modulus is increased with increasing segment number and decreasing with . Considering , a large ratio is found under condition n = 10,000 and = 10 (E is the elastic modulus equaling random distribution), as shown in Figure 6a. However, considering a dimension , the ratio is much smaller than that of —especially = 10,000, as shown in Figure 6b–e.

Figure 6

Effect of geometric effect on the equivalent elastic modulus: (a) = 1; (b) = 2; (c) = 10; (d) = 500; (e) = 10,000.

3.2. Simulation of Rectangle Sample

As an example, tensile tests of VHB 4910 (Acrylic, 3M Company, USA) specimens were conducted. The specimen (25 mm × 1 mm × 30 mm) was clamped between the grips, then the specimen was loaded at a given rate. A rectangular region (20 mm × 20 mm) was marked in the center region of the sample. By doing so, the deformation mechanism was investigated, as shown in Figure 7.

Figure 7

Illustration of VHB 4910 specimen under tension test [48].

The area of the marked region was recorded during the loading process and compared with the simulation results of the Yeoh model and the modified model, as shown in Figure 8. During the simulation process using the Yeoh model, the material was treated as homogeneous isotropic, the parameters and more details can be found in Reference [47] and in our previous work [49]. The minimum values of l and n were 0.08 mm and 15 for the modified model. The deformation of the sample can be summarized as: when strain ranged from 0 to 4, a homogeneous deformation could be observed; however, when strain ranged from 4 to 12, the deformation was affected by shrinking behavior of the boundary. As can be seen from Figure 8, the comparison between the simulation results of the modified and the experimental data shows good agreement with relative error 5.2% (28% for Yeoh model).

Figure 8

Relationship between the marked area and tensile strain [48].

The stress–strain curves of VHB 4910 specimens under different loading speeds (0.05/s and 0.22/s) are shown in Figure 9. It was observed that VHB 4910 demonstrated typical elastomeric tensile behavior, which means a greater stress could be observed under larger strain or larger strain rate. In a state of small strain ( < 4), both Yeoh model and the modified model were valid for the nonlinear mechanical behavior of VHB 4910. However, in a state of large strain state ( > 4) only the modified model showed a promising result (ignoring the difference near the failure strain, 13).

Figure 9

Relationship between stress and strain under different loading rates [48].

The difference of the two models can be explained as follows: the middle part of the specimen was shrinking during the loading process. This part reduced quickly using the modified model compared with the Yeoh model (especially for larger strain ). Additionally, the equivalent elastic modulus calculated using the modified model (near the boundary region) became much larger compared with the Yeoh model (along the z-direction). Because of this, the material can be regarded as nonhomogeneous, which can further affect the deformation of the material.

3.3. Simulation of Circular Dielectric Elastomer

A circular elastomer membrane with radius R was fixed on a movable frame and coated with flexible electrode on both sides (radius = 7.5 mm) in the reference state. Then, the membrane was pre-stretched to a level = 75 mm in the pre-stretched state, as shown in Figure 10a. In the activated state, the radius of active area changed into by applying an external voltage V on the active area, as shown in Figure 10b.

Figure 10

Working principle of a circular dielectric elastomer actuator [47]: (a) Pre-stretched state; (b) Activated state.

The relationships between radial strain of active area ( − 1) and time under different pre-stretched ratio and applied voltage V are plotted in Figure 11.

Figure 11

Fitting the improved model and Yeoh model to the circular dielectric elastomers: (a) = 3; (b) = 4; (c) = 5.

From Figure 11, some conclusions can be drawn as follows. Subjected to the same voltage, the radial strain is increased with the pre-stretched level. When the pre-stretched strain was held constant, the radial strain increased with increasing voltage. The radial strain changed temporality and quickly when at time t < 300 s, then remained stable over a longer period, t >> 300 s. The relative errors of the Yeoh model were 18% and 21% for test = 4, V = 3500 V and test = 5, V = 3000 V, respectively, as shown in Figure 11b,c. However, the relative errors of the modified model were 4.6% and 6.2%, which are much more acceptable.

As for lower pre-stretched level and voltage, the simulation results of the two models were close, as can be seen from Figure 12a,b. However, at higher pre-stretched level and voltage, the boundary of the activated area extended with time more quickly using the modified model compared with the Yeoh model, as shown in Figure 12c,d.

Figure 12

Simulation results of quarter model of the circular dielectric elastomer at time t = 300 s: (a) Yeoh model = 3 and V = 2000 V; (b) Modified model = 3 and V = 2000 V; (c) Yeoh model = 5 and V = 3000 V; (d) Modified model = 5 and V = 3000 V.

4. The Effect of Internal Particles on the Electromechanical Properties

4.1. Determination of the Material Parameters: Local Density

In this section, the local distribution of polymer chains is investigated, providing a basis for the next section. Particles with radius r are regarded as homogeneous distribution inside the matrix ideally, as the representative volume element (length d) shown in Figure 13a. The advantage of this assumption is that the computational time can be saved rather than considering the real microstructure of the particle reinforced polymer. From Figure 13b–d, the effect of volume fraction v = , segment number n, particle radius-to-segment length ratio , and parameter k on the local density (P) was studied and compared with that of random distribution (Prand). The calculation results show that local density increased with increasing k and decreasing . From Figure 13a,e,f, the local densities of different positions (points A, B and C) were calculated. The calculation results indicate that local density was larger when the point was closer to the boundary. Beyond that, local density increased with segment number n and volume fraction v, as shown in Figure 13.

Figure 13

Density of polymer chain: (a) representative volume element; (b) point A- = 10, k = 1; (c) point A- = 10, k = 10; (d) point A- = 100, k = 1; (e) point B- = 10, k = 1; (f) point C- = 10, k = 1.

4.2. Simulation of the Physical Parameters

The effects of internal particles on material electromechanical properties are mainly discussed in this section. The particles (CCTO) were regarded as homogeneously distributed inside the material. With different particle volume fractions, the reinforcements of the particles are plotted in Figure 14, from which the equivalent elastic modulus can be observed as a function of l, n, and k.

Figure 14

The effect of internal geometry on the equivalent elastic modulus: (a) = 100; k = 1; (b) = 10; k = 1; (c) = 1; k = 1; (d) = 10; k = 2; (e) = 10; k = 10.

As can be seen in Figure 14, the equivalent elastic modulus generally increased with increasing segment number n. With a large particle radius-to-segment length ratio ( = 100), the equivalent elastic modulus increased with the volume fraction, as shown in Figure 14a. Differently, with a small particle radius-to-segment length ratio ( = 10 and = 1), the equivalent elastic modulus increased with the volume fraction over a range of segment numbers and then decreased, as shown in Figure 14b,c. When comparing Figure 14b,d, the reinforcement of cross-linking characteristic can be observed: the equivalent elastic modulus was generally enhanced (about 10%). In Figure 14e, the reinforcement is more obvious than that of k = 2:70% for volume fraction 0.6; 10% for volume fraction 0.05.

In addition, the effect of internal particles on the dielectric property was evaluated. According to the research on CCTO–PDMS composites, the dielectric permittivity for CCTO after sintering was almost 105 under a frequency of 105 Hz; the dielectric permittivity for CCTO after the grounding process was 102–103 under frequency < 102 Hz (chosen as ); for PDMS, the dielectric permittivity almost held constant (chosen as ) over a broad frequency range (10−1–106 Hz). In contrast to the equivalent elastic modulus mentioned above, the dielectric permittivity was mainly affected by volume fraction, although there was a small increase in at n = 106 compared to that at n = 1–10 (about 10%), as shown in Figure 15. The main reason for the difference is the colossal dielectric permittivity of CCTO.

Figure 15

Effect of internal filler on the dielectric permittivity of polymer-based composite.

4.3. Simulation of CCTO-PDMS

The details of the fabrication and the test process of CCTO-PDMS with different volume fractions of CCTO (2.3%, 5.1%, and 8.4%) can be found in the literature [25]. In this section, the parameters = 25, n = 42, and k = 2 were employed in the simulation. For comparison purposes, the calculation results of Lichtenecker’s logarithmic law and Yamada’s model of mixing of composites are added in Figure 16a. Using the same parameters, the electromechanical behavior of CCTO-PDMS was also calculated, as shown in Figure 16b. The Lichtenecker logarithmic law and Yamada’s model can be expressed as [25,50]:

Figure 16

Simulation results of the CaCu3Ti4O12 (CCTO)–PDMS: (a) dielectric permittivity; (b) electromechanical behavior [25].

5. Conclusions

In summary, we have proposed a numerical method to capture the electromechanical behavior of polymer-based composites. The distribution of polymer chains was modified considering several factors: segment number, segment length, external geometries, and internal fillers. Generally, inside the system the chains were regarded as random distribution; however, near any surfaces (inside or outside), the distribution of the chains was restricted. The advantage of the model is that many physical parameters related to local polymer chain can be investigated.

As for VHB 4910 film, it is suggested that the polymer chain is restricted along the thickness direction. The simulation results show that when the deformation increased due to electromechanical coupling, the activated part had the potential to increase the instability and further deform owing to the lower corresponding elastic modulus. More importantly, the electromechanical coupling effect was increased by this phenomenon compared with the isotropic assumption. For a complex system, it seems that the electric field is a function of the surface characteristics, which change during the loading process.

For the CCTO–PDMS composite, the calculation results indicate that the electromechanical properties are greatly affected by internal fillers. With a higher density of CCTO, a higher elastic modulus and dielectric was generally found. Moreover, interaction between the filler and the matrix also contributes to the reinforcement, up to a critical level. The proposed model was successfully used to estimate the experimental data in the literature, and can promote the design and development of composite materials.