Electrostatic model of dielectric elastomer generator based on finite element

When dielectric elastomer materials are used for power generation, bias voltage is applied at both ends of dielectric elastomer film, and there are equal amounts of heterogeneous charges on both sides of the film, so Maxwell electrostatic force is always coupled in the process of power generation. In order to investigate the distribution of Maxwell stress in dielectric elastomer material under electric field, the electrostatic model of dielectric elastomer generator is established in COMSOL finite element simulation software environment in this paper. The distribution of electrostatic force is studied from two aspects of theoretical derivation and simulation, and the magnitude and direction of electrostatic force are determined. The simulation results show that the Maxwell electrostatic force can be equivalent to the tensile force along the film plane and the extrusion force perpendicular to the plane, and they are the same.

Dielectric elastomer (DE) is a new type of smart material that could be used as an actuator to convert electrical energy into mechanical energy. In recent years, it has been studied in the field of power generation for energy recovery [1][2][3][4] . DE has some advantages, such as large strain, light weight and high flexibility [5][6][7][8] . Compared with other power generation materials, DE has higher energy conversion efficiency and specific energy density 9,10 . Therefore, DE has good application prospects and application fields, such as wind generators, wave energy generators, body-worn generators, etc. 11,12 .
Dielectric Elastomer Generator (DEG) refers to a "sandwich" three-layer variable capacitor that uses a dielectric elastomer as a carrier, two layers of upper and lower flexible electrodes, and an electrode wire. Its power generation principle is variable capacitor power generation 13 . Under the action of external mechanical force, the capacitance of DE increases when DE is stretched. At this time, under the action of bias power supply, the initial charge is increased. When DE shrinks, the anisotropic charges in the upper and lower surface electrodes are pushed away due to the increase of thickness, and the isotropic charges are squeezed closer due to the decrease of area, thus increasing the charge voltage. That is, the capacitance decreases, the charge remains unchanged, and the voltage increases, so that the mechanical energy can be converted into electrical energy 14 .
Dielectric elastomer, as an intelligent energy exchange material, can be used not only in bionic driving materials, but also in the field of power generation. When used in bionic driving materials, Maxwell electrostatic force is the essential driving source of driving materials. Applying voltage to the upper and lower electrodes of the film, the upper and lower electrodes accumulate equal charges with opposite polarities, which results in Maxwell electrostatic force driving film deformation. When applied in the field of power generation, bias voltage is applied at both ends of the film, and the same amount of heterogeneous charges are accumulated on both sides of the film, so Maxwell electrostatic force is always coupled in the process of power generation, which affects the electromechanical characteristics of DEG [15][16][17][18] .
In order to investigate the distribution of Maxwell stress in dielectric elastomer materials under electric field during power generation 19,20 , the electrostatic model of dielectric elastomer materials was established in COM-SOL finite element simulation software environment based on theoretical deduction. Based on the simulation results, the electromechanical characteristics of DEG are studied.

Theoretical Model
Electromagnetic field essentially generates volume force in the medium. The resultant force of electromagnetic field acting on the medium can be calculated by the method of volume force density, or by the surface force produced by Maxwell stress tensor 21 . The Maxwell stress tensor can be expressed as: where E is the electric field strength and D is the potential shift. And E is the gradient of voltage, which can be expressed as: where V is the voltage. Potential shift D can be expressed by electric field strength E.
where ε 0 is the dielectric constant in vacuum and ε r is the relative dielectric constant. Thus, formula (1) can be simplified as follows: The surface force produced by Maxwell stress tensor in DE film under electric field is given below, as shown in Fig. 1. Let O be a point on a surface of thin film and establish coordinate system XYZ through O point 22 . The normal unit vector of O-crossing point is n , and the electric field E and Z axes are in the same direction. In this coordinate system, the three components of electric field E are: Then Maxwell stress tensor can be expressed in matrix form: The electric field force acting on the unit area of the film is as follows: If n and X, Y and Z axes are in the same direction in turn, the electric field force components in the three main directions can be obtained as follows: where P ex , P ey and P ez components are the projections of P e along the x-axis, y-axis, and z-axis, respectively. According to the projection result, P ex component and P ey component are positive, and P ez component is negative.
Next, from the point of view of energy conservation, the force acting on DE film under electric field is analyzed.
The DE film coated with electrodes is analogized to a deformable parallel plate capacitor. Therefore, when the capacitor is charged, an equal amount of dissimilar charge Q is accumulated at both ends of the positive and negative electrodes. At this time, the capacitance of the film can be expressed as: www.nature.com/scientificreports/ Among them, U denotes the voltage at both ends of the capacitor, S denotes the area of the capacitor plates, and d denotes the distance between the capacitor plates.
According to the formula of capacitive electric energy: W = 1 2 CU 2 , when the external power supply charges the DE film, the electric energy is: After applying voltage, the thickness and area of the film are deformed due to the existence of Maxwell electrostatic force, so the increment of W e is as follows: The first term in the formula represents the energy change caused by the change of charge when the power supply charges the thin film, and the second term represents the energy differential which is converted into mechanical energy during the deformation of the thin film. In the charging process, the work done by charge and external force is equal to the change of electric energy on DE film.
where σ EA is the equivalent stress of Maxwell stress in the plane direction and σ EZ is the equivalent stress of Maxwell stress in the thickness direction.
Comparisons (11) and (12) are available: It can be simplified as follows: The above two methods show that the Maxwell stress has σ EA in the plane direction and σ EZ in the thickness direction, and both of them are equal in size and are 1 2 ε 0 ε r E 2 .

Simulation Study
The flow chart for numerical simulation process of DEG electrostatic model based on finite element is shown in Fig. 2. DE material will undergo stress and deformation under electric field, and there is a coupling between force field and electric field. From the above theoretical analysis, it is known that the Maxwell stress is generated when the electric field is applied. COMSOL Multiphysics finite element software is famous for its strong multiphysical field coupling characteristics. It can easily couple various physical fields and analyze their coupling characteristics 23 . COMSOL Multiphysics 5.0 was used in the study. In the electrostatic module of COMSOL Multiphysics, the electric field and force field are coupled based on Maxwell stress tensor. The formulas are as follows: www.nature.com/scientificreports/ In order to study the Maxwell stress distribution of DE material under electric field, an electrostatic model of DE material was established in COMSOL environment. In order to study the convenience, a rotationally symmetrical circular thin film sample is selected for analysis. The schematic diagram is shown in Fig. 3. The electrodes of the same size are coated on the upper and lower sides of the circular DE matrix. The DC voltage is applied to the upper and lower electrodes. The whole material is in an infinite air region, thus forming the whole electrostatic model.
As we know from the above figure, the whole model is rotationally axisymmetric. In the finite element modeling, the electrostatic model can be reduced to a two-dimensional rotationally symmetric model for easy solution and analysis. The schematic diagram is shown in Fig. 4. As can be seen from the figure, after applying voltage, the upper and lower electrodes accumulate the same amount of heterogeneous charges. The heterogeneous charges  www.nature.com/scientificreports/ gathered at the upper and lower electrodes attract each other and produce force F z , while the same charges on the same side of the electrode repel each other and produce force F r . Among them, r d and z d are the radius and thickness of DE matrix, r a and z a are the radius and thickness of coated electrode respectively, and the bias voltage is V. The electrostatic module of COMSOL Multiphysics is selected for two-dimensional rotational symmetry modeling. Firstly, the geometric model is established. The radius r d of DE matrix is 14 mm, the thickness z d of DE is 60 um, the radius r a of electrode is 7 mm, and the thickness z a is 20 um. All materials are contained in the circular air region, as shown in Fig. 5. The relative dielectric constant of DE film is 4.7 and that of air is 1. Physical field constraints are applied to the model, 3000 V DC voltage is added to the upper electrode, and the lower electrode is grounded to form potential difference. The application of voltage causes the charge to accumulate on the surface of the upper and lower electrodes, and the surface charge density is ρ . The divergence of the electric field strength is equal to the charge density divided by the dielectric constant. Substituting charge density = dielectric permittivity of medium × electric field strength. So Maxwell stress tensor formula can be expressed by charge density as follows: where ρ s is the charge density.
In order to improve the simulation accuracy, the DE matrix and electrodes are divided into regular quadrilateral meshes, and the boundary is refined. The triangular meshes are used in the external air domain, as shown in Fig. 6.
The steady-state solution of the model is carried out and the post-processing analysis is carried out.

Results Analysis
As mentioned above, the charge is concentrated on the surface of the electrode. Taking the upper electrode as an example, the surface of the electrode is analyzed, as shown in Fig. 7. The simulation results show that the charge density on the outer surface is 3.93 × 10 −11 , the charge density on the side is 3.93 × 10 −11 , and the charge density on the inner surface is 0.00208. It can be seen that the charge density on the outer surface is two to three orders of magnitude different from that on the side and the inner surface. According to formula (17), the surface stress tensor of the outer surface is almost negligible, and the stress mainly concentrates on the inner surface and the side. The z-direction stress tensor of the inner surface is 52,222 N/m 2 , while the r-direction stress tensor is 3.93 × 10 −11 N/m 2 . It can be seen that the Maxwell stress tensor on the inner surface is extruded inward from the vertical film surface, while the r-direction stress is almost zero, which can be neglected. The stress of the lower electrode is the same.
The simulation results show that the surface stress tensor in the r direction of the upper electrode side is 25,882 N/m 2 , while the surface stress tensor in the lower electrode side is the same, so the surface stress tensor It can be concluded that the compressive stress p z of the whole circular film is 52,222 N/m 2 perpendicular to the plane direction of the film, and the tensile stress p r of the plane direction is 51,764 N/m 2 . The stress calculated by formula (17) is 52,018 N/m 2 . Therefore, the following relationship can be obtained by removing the error of simulation solution:

Conclusions
In this paper, DEG electrostatic model is established, and the distribution of Maxwell stress is determined by theoretical deduction and simulation analysis, and a unified result is obtained. The following conclusions are drawn: 1. According to the theory of electromagnetics and conservation of energy, it is deduced that the film is subjected to the compressive stress in the plane direction and the tensile stress in the plane direction, both of which are 1 2 ε 0 ε r E 2 .
(18) p r = p z = 1 2 ε 0 ε r E 2  www.nature.com/scientificreports/ 2. The electrostatic finite element model is established in COMSOL electrostatic module. By applying physical field constraints, the charge density and surface stress tensor on the electrode surface are obtained. Similarly, the compressive stress and the tensile stress in the plane direction of the vertical film are obtained, both of which are 1 2 ε 0 ε r E 2 .