## Abstract

Electroactive polymers (EAPs) deform when subject to an electric field, which is generated by two or more electrodes. To ensure proper function of the EAP, these electrodes are connected to a source and they are therefore required to be continuous such that no isolated islands exist. Increasing an EAP’s performance using topology optimization while ensuring electrode connectivity is the goal of this work. A topology optimization formulation is introduced where electrode connectivity is ensured using the virtual temperature method. Numerical experiments demonstrate that this is an efficient method to guarantee connectivity.

## 1 Introduction

Electroactive polymer (EAP) is a class of materials that, when electrically stimulated, reacts with mechanical deformation. Modeling and design of EAP is an active research area and it has great potential for a number of applications, e.g., artificial muscles and soft robotics [1]. A typical EAP is a dielectric actuator constructed by placing a thin layer of EAP material between two compliant electrodes. When a potential difference is applied to the electrodes, the EAP material will contract in the direction of the electric field and elongate in the in-plane directions.

The deformation in a dielectric polymer plate with electrodes on both sides is mainly due to two physical effects. For large potential differences, the two electrodes will attract each other due to Coulomb forces and thus contract the plate. This is known as the Maxwell effect. Too high electric fields may cause dielectric breakdown, which results in that the dielectric material looses its insulating properties and becomes conducting. Another effect that is prevalent at moderate electric fields is called (inherent) electrostriction. Dielectric polymers have a quadratic relationship between the deformation and the electric field. Reversing the electric field will thus yield the same deformation. This is in contrast to piezoelectric materials, which instead have a linear relationship where the deformation depends on the sign of the electric field [1].

Topology optimization (TO) of multi-physics phenomena has been performed for several different types of couplings, such as thermo-mechanical in, e.g., Granlund et al. [2] and Zhu et al. [3], or thermo-electric in, e.g., Xu et al. [4] and Xing at al. [5]. TO has also been applied to electro-mechanical coupling, like piezoelectric materials in, e.g., Lee and Tovar [6]. However, the research on EAPs is scarce. Ortigosa et al. [7] used TO to optimize the active EAP layer. Based on this model, Ortigosa and Martínez-Frutos [8] developed a multi-resolution model where multiple density values were used in each element. In Ref. [9], the electrode placement was optimized, however not considering electrode connectivity. As a result, some of the optimal electrode layouts were discontinuous. Recently, Ortigosa et al. [10] combined a multi-material approach for the active layer and shape-morphing TO. For the specific case of dielectric elastomers, Bortot et al. [11] used TO to design tunable band gaps using a genetic algorithm. Recently, this problem was solved using gradient-based optimization, see Sharma et al. [12]. In all of the mentioned cases, the influence of the electric field in the surrounding free space was not accounted for. This was very recently investigated by Dev et al. [13].

In certain TO applications, it is crucial to make structures continuous. The virtual temperature method (VTM) by Liu et al. [14], addresses this in the context of additive manufacturing, where a continuous void region is necessary to be able to remove excess support or unused powder material. The method is also denoted as the virtual scalar field method, see Li et al. [15]. The method is based on solving an auxiliary, fictitious thermal problem, where regions with low density have high conductivity and high heat supply, whereas regions with high density have low conductivity and low heat supply. Homogeneous Dirichlet condition is enforced on part of the boundary while the rest is subjected to a homogeneous Neumann condition. Consequently, low density regions in connection with the boundary with the Dirichlet boundary condition will have low temperatures while regions not connected will have a high temperature. By limiting the maximum fictitious temperature, connectivity can be promoted. A modified version of the VTM was used in this work, where the behavior of the solid and void regions has been swapped. A similar approach was used in, e.g., Swartz et al. [16].

Manufacturing constraints in periodic structures were enforced by VTM in Swartz et al. [16] and electrode connectivity for piezo modal transducers was ensured in Donoso and Guest [17]. Luo et al. [18,19] extended VTM to also incorporate a nonlinear heat source term which results in uniform temperature in enclosed voids, regardless of void sizes and locations.

An alternative to VTM is to use spectral graph theory on the discretized density distribution, which Donoso et al. [20,21] used to ensure electrode connectivity in piezo transducers. Recently, the same authors also developed a similar method for a continuous density distribution in Ref. [22].

In EAP structures, electrodes must be connected to a source and it is therefore necessary to ensure that the electrodes are continuous. Combining EAP and TO that ensures electrode connectivity is the main goal of the present work. The shape of the electrodes follows the shape of the active layer. Optimizing the shape of the active layer will therefore implicitly dictate the shape of the electrodes.

The structure of this paper is as follows. In Sec. 2, the EAP model is described. In Sec. 3, the TO formulation is developed. Of special interest to this work is Sec. 3.3, which describes a modified version of the previously mentioned VTM. Results from developed methods are shown in Sec. 4. Lastly, Sec. 5 provides conclusions.

## 2 Continuum Model for Electroactive Polymers

The basic relations that govern electro-mechanically coupled problems will be introduced, and restricted to the electrostatic case. For a more detailed background information, the reader is referred to Jackson [23], Eringen and Maugin [24], Kovetz [25], and Dorfmann and Ogden [26].

The considered body is a collection of material particles labeled by their coordinates $X$ in the material configuration, $\Omega 0$. The motion of the body is described by the smooth mapping $\phi (X,t)$ of points $X$ in $\Omega 0$ to their position $x=\phi (X,t)$ in the current domain $\Omega t$. The deformation gradient $F=\u2207X\phi $ locally describes the deformation and $J=det(F)$.

### 2.1 Balance Relations.

### 2.2 Constitutive Model.

### 2.3 Variational Forms and Finite Element Discretization.

The variational forms are discretized using a standard nonlinear finite element formulation.

## 3 Topology Optimization Formulation

A piecewise uniform density field $\rho $ is introduced in the design domain $\Omega DD$, which is a subset of the undeformed domain $\Omega 0$. Each discretized finite element $e$ in $\Omega DD$ is associated with a uniform element density $\rho e\u2208[0,1]$. The element densities are considered to be the design variables for the optimization and are collected in the vector $\rho $.

### 3.1 Filter, Projection, and Material Interpolation.

### 3.2 Optimization Formulation.

### 3.3 Connectivity Constraints.

For proper functionality, the electrodes must be connected. To ensure this connectivity, the electrodes should be continuous. If not, i.e., when islands are formed, the connectivity condition requires extra threading which is unfavorable from a manufacturing point of view. To illustrate this, consider the two structures in Fig. 1, where the top consists of two disconnected islands whereas the bottom structure is connected. As creation of islands should be prevented, a connectivity constraint is required.

The VTM [14], also known as virtual scalar field method [15] is here used to enforce connection of solid material between the electrodes. In this method, an auxiliary fictitious thermal problem is solved in the electrodes. Solid regions are assigned high heat supply and high conductivity whereas void regions have low heat supply and low conductivity. Homogeneous Dirichlet boundary condition is enforced on the part of the boundary where the electrode is externally connected, $\u2202\Omega eliTg$. The rest of the boundary, $\u2202\Omega eliTh$, is subject to homogeneous Neumann conditions. Solving this auxiliary problem for Fig. 1(a) configuration will result in a low temperature for the region connected to the Dirichlet part of the boundary due to high heat transport out of the domain. The right solid region in Fig. 1(a) is isolated and heat can not leave and will therefore have a high temperature.

In Fig. 1(b) configuration, the two regions are connected which enables conduction from right to left and as a result the temperature will be lower. By limiting the maximum temperature, connectivity can thus be ensured.

^{3}are arbitrary since they are formulated in relation to a reference temperature. The value of $\u03f5$ should be small but sufficiently large to ensure that the stiffness matrix does not become singular. In this work, $\u03f5=10\u22125$ was used.

### 3.4 Sensitivity Analysis.

## 4 Numerical Examples

The material properties for the solid EAP material are given in Table 1 and the void region is modeled like solid properties scaled by a factor of $10\u22129$ to reduce the influence of the void region, while avoiding a singular stiffness matrix.

$K$ (MPa) | $G$ (MPa) | $ce$ (NV^{-2}) | $\epsilon r(\u2212)$ |
---|---|---|---|

0.6 | 0.1 | 1 | 4.7 |

$K$ (MPa) | $G$ (MPa) | $ce$ (NV^{-2}) | $\epsilon r(\u2212)$ |
---|---|---|---|

0.6 | 0.1 | 1 | 4.7 |

Inspired by Ortigosa et al. [7], the electro-mechanical coupling term is penalized more than the purely electrical and mechanical terms to ensure low electro-mechanical coupling in the void and intermediate regions. The penalization exponents used herein are given in Table 2.

Figure 2 shows a quarter of a cylinder consisting of two layers. The outer part is the active layer and the inner part the passive host layer. The active layer is considered as the design domain while the host layer is fixed with only solid material. The electrodes are situated on both sides of the active layer, between which the potential difference $\Delta \varphi $ is applied. The sources for the electrodes are situated along the bottom edges. It is also assumed that $DN=0$. Symmetry is enforced on the $x=0$ and $z=0$ planes. The objective is to maximize the distance in $y$-direction between the top and bottom edges. Multiple cylinders with different lengths $lz$ were used and the other geometrical parameters can be seen in Table 3. In all cases, a volume fraction $\alpha =40%$ is used. All finite element models are based on 3D solid elements.

Case | 1 | 2 | 3 |
---|---|---|---|

$r$ (mm) | 40 | 40 | 40 |

$lz$ (mm) | 40 | 80 | 120 |

$t1$ (mm) | 1 | 1 | 1 |

$t2$ (mm) | 1 | 1 | 1 |

$\Delta \varphi $ (kV) | 15 | 15 | 15 |

Elements | $4\xd740\xd7132$ | $4\xd780\xd7132$ | $4\xd7120\xd7132$ |

Case | 1 | 2 | 3 |
---|---|---|---|

$r$ (mm) | 40 | 40 | 40 |

$lz$ (mm) | 40 | 80 | 120 |

$t1$ (mm) | 1 | 1 | 1 |

$t2$ (mm) | 1 | 1 | 1 |

$\Delta \varphi $ (kV) | 15 | 15 | 15 |

Elements | $4\xd740\xd7132$ | $4\xd780\xd7132$ | $4\xd7120\xd7132$ |

The objective is to maximize the separation between the top and bottom edges, as indicated in Fig. 2 by $uout$ and the thick lines and arrows. The vector $l$ is defined such that only displacements in the $y$-direction on the top and bottom edges are nonzero,

The filter and projection parameters were chosen based on numerical experiments and determined to be appropriate. As an example, a too low value on the filter parameter $l0$ could cause issues by creating too thin connections. Therefore, it was set to be approximately four element sizes. The thresholding parameter $\beta $ was initiated to $\beta =1$ and then increased by 1 every 20th iteration to a maximum of 15. The parameter $\eta =0.5$ was fixed.

Standard parameters for the MMA were used, except for the first 300 iterations where instead $GHINIT=0.1$, $GHDECR=0.9$, and $GHINCR=1.1$ were used. The objective function $g0$ was scaled by a factor $103$ and the constraint function $g1$ by 10. After the continuation ends, convergence was checked and the design updates continues until $|g0k\u2212g0k\u22121g0k|<10\u22125$, and until all constraints were fulfilled.

Algorithm 1 presents in pseudo-code the main steps for implementing the presented methodology.

### Optimization algorithm without connectivity

Set initial values on all parameters

Set $\rho $ for the initial design

**While**$|g0k\u2212g0k\u22121|>10\u22125|g0k|$ and continuation terminated do

If applicable, update parameters $\beta $ and $pk$,

Apply filter to $\rho $ to obtain $\rho ~$ from Eq. (32)

Project filtered densities $\rho ~$ to obtain $\rho \xaf$ from Eq. (33)

Solve state problem in Eq. (24) using Eq. (28)

Calculate constraints $gi$, $i=1,\u2026,nconstr$

Compute sensitivities using the adjoint method

Solve MMA to update element densities $\rho $ using MMA

**end While**

### 4.1 No Connectivity Constraint.

The optimized designs in Fig. 3, consist of two symmetrical but separate islands. Clearly, the optimizer favors discontinuous designs.

Another interesting result is that all designs have finger-like structure. Figure 4 compares case 1 in Table 3 with a similar structure, also with $\alpha =40%$ where the electrodes are straight and symmetrically placed on top and bottom sides. The optimized structure with ”fingers” has about $3%$ better performance than the straight electrode, verifying that the finger-like structure indeed is favorable for the performance.

### 4.2 Connectivity Constraint.

The previous result clearly shows that discontinuous electrodes improve the performance for the described problem, however, as the source for the electrodes are located at the bottom of the cylinder, isolated islands are obtained.

The penalization exponent for the heat supply in the virtual temperature problem was fixed to $pQ=2$. The penalization exponent for the thermal conductivity started at $pk=2$ and increased by 20% every 20th iteration until a maximum $pk=10$ is reached. In the p-norm for the fictitious temperature in Eq. (41), $p=10$ was used according to Donoso et al. [17]. The upper bound for the constraint was fixed at $\lambda max=10$.

The optimal designs for the geometries in Table 3 are presented in Fig. 5 where it is clear that at least one solid material connection exists between top and bottom sides. The connectivity renders 6–9% decrease in performance and similar to the previous example, finger-like structures are formed.

Since the upper limit for the virtual temperature constraint, $\lambda max=10$, is arbitrary other choices could be suitable. In Fig. 6, a comparison of case 1 in Table 3 optimized according to Eq. (53) for different $\lambda max\u2208{6,10,14}$ is shown. In this comparison, $\lambda min$ is initiated to 0 and then changed to $\lambda min=\lambda max\u22120.01$ at iteration 40. In all cases, a clear connection is formed between the upper and lower parts. The material distribution is similar for $\lambda max=6$ and $\lambda max=10$ while $\lambda max=14$ gives a noticeable difference. The fingers previously shown in Figs. 3 and 5 are again present.

This numerical experiment shows that the choice of $\lambda max$ can be made moderately arbitrarily but for all cases, having the constraint present generates a connected structure. By inspection of $g0$ at the last iteration and the convergence plot in Fig. 6(d) it can be seen that higher $\lambda max$, in general, results in lower, more negative $g0$, meaning larger displacements. This is expected. Higher $\lambda max$ allows for higher temperatures and therefore thinner connections, which in turn means that more material can be placed on the top and bottom parts. This is obviously more beneficial, as can be seen in the case with no connection in Fig. 3 where all of the material is placed on the top and bottom parts.

## 5 Conclusions

This work has presented a density-based topology optimization framework for EAPs ensuring continuously connected electrodes. The developed method for connectivity, based on the virtual temperature method, has been shown to provide a simple and efficient procedure for ensuring electrodes that are connected to a given source. Enforcing connectivity decreased the performance by 6–9%, which is a tradeoff to obtain a functioning structure.

Numerical proof-of-concept examples were presented to assess the need for connectivity constraints and to show the presented methodology. Only constraining the maximum virtual temperature in the virtual temperature method, was found to be insufficient due to unstable behavior in the optimization iterations. To stabilize the optimization, a lower bound for the temperature $Tp$ was introduced. Although one should note that this lower bound constraint was not active in the final designs and used solely for stabilizing the optimization process.

The value of the upper bound for the virtual temperature was shown to not greatly affect the connectivity, but as expected higher values resulted in better performance. Also, multiple connections were created between the electrode parts due to only a quarter of the full structure being used in the optimization. For a physical electrode, it is sufficient with only one connection between the different electrode parts. This could be realized by modeling the full structure.

A natural continuation based on this work would be to directly optimize the placement of the electrodes, possibly by introducing the electrodes as an additional material to the EAP. This would provide more design freedom for the optimization and thereby improved performance.

## Acknowledgment

The Swedish Research Council, Grant No. 2021-03851, and the Swedish Energy agency, Grant No. 48344-1, are gratefully acknowledged. The computations were enabled by resources provided by the Swedish National Infrastructure for Computing (SNIC) at Lunarc, partially funded by the Swedish Research Council through Grant No. 2018-05973. Finally, we would also like to thank Professor Krister Svanberg for his implementation of MMA in fortran.

## Conflict of Interest

On behalf of all authors, the corresponding author states that there is no conflict of interest.

## Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

## Nomenclature

- $k$ =
thermal conductivity (W/m/K)

- $p$ =
penalization exponent, parameter for p-norm

- $r$ =
cylinder radius (m)

- $t$ =
time (s)

=*a*solution vector

=*f*mechanical body load

=*l*constant vector for objective function

=*r*residual

=*x*coordinates in the spatial configuration (m)

- $D$ =
material tangent

- $G$ =
shear modulus (Pa)

- $J$ =
Jacobian determinant of deformation gradient

- $K$ =
bulk modulus (Pa)

- $Q$ =
heat supply (W/m

^{3})- $V$ =
volume (m

^{3})=*B*derivatives of shape functions

=*C*Cauchy–Green tensor

- $ce$ =
electric parameter (N/V

^{2})=*D*electric displacement (C/m

^{2})- $DN$ =
prescribed electric charge on boundary (C/m

^{2})=*E*electric field (V/m)

=*F*deformation gradient

- $T$ =
temperature (K)

=*I*identity tensor

=*K*stiffness matrix

=*N*outward unit normal in material configuration, shape functions

=*T*total first Piola-type stress tensor

=*X*coordinates in the material configuration (m)

- $g0$ =
objective function (m)

- $gi$ =
constraint function

- $k0$ =
unit thermal conductivity (W/m/K)

- $lz$ =
cylinder length (m)

- $l0$ =
length scale in PDE filter (m)

- $t1$ =
thickness (m)

- $t2$ =
thickness (m)

- $rc$ =
residual for VTM

- $Kf$ =
stiffness matrix for PDE filter

- $Kc$ =
stiffness matrix for VTM

- $Q0$ =
unit heat supply (W/m

^{3})- $Tp$ =
approximate maximum temperature (K)

- $Tref$ =
reference temperature (K)

- $VDD$ =
volume in design domain (m

^{3})- $TN$ =
prescribed traction

### Greek Symbols

- $\alpha $ =
allowed volume fraction

- $\beta $ =
sharpness parameter

- $\epsilon 0$ =
vacuum permittivity (F/m)

- $\epsilon r$ =
relative permittivity

- $\eta $ =
threshold parameter

- $\lambda $ =
adjoint vector

- $\lambda c$ =
adjoint vector for VTM

- $\lambda $ =
factor in VTM constraint

- $\lambda max$ =
maximum factor in VTM constraint

- $\lambda min$ =
minimum factor in VTM constraint

- $\rho $ =
element density

- $\rho ~$ =
filtered element density

- $\rho \xaf$ =
projected element density

- $\varphi $ =
scalar electric potential (V)

- $\varphi g$ =
prescribed scalar electric potential (V)

- $\phi $ =
deformation mapping (m)

- $\phi g$ =
prescribed displacement (m)

- $\Omega $ =
augmented free energy function

- $\Omega m$ =
mechanical free energy function

- $\Omega mel$ =
electro-mechanical free energy function

- $\Omega el$ =
electric free energy function

- $\Omega 0$ =
material configuration

- $\u2202\Omega 0$ =
boundary to material configuration

- $\Omega t$ =
spatial configuration

- $\u2202\Omega t$ =
boundary to spatial configuration