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.

Graphical Abstract Figure
Graphical Abstract Figure
Close modal

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, Ω0. The motion of the body is described by the smooth mapping φ(X,t) of points X in Ω0 to their position x=φ(X,t) in the current domain Ωt. The deformation gradient F=Xφ locally describes the deformation and J=det(F).

The hyper-elastic material response will be governed by the right Cauchy–Green tensor which is defined as
C=FTF
(1)

2.1 Balance Relations.

In the electro-quasi-static case, the boundary Ω0, with the outward unit normal N to Ω0, is split such that Ω0=Ω0h,ϕΩ0g,ϕ and Ω0h,ϕΩ0g,ϕ=. The Maxwell–Faraday equation for quasi-electrostatics reads
X×E=0
(2)
which implies that the electric field E is conservative and consequently is obtained from the gradient of a scalar electric potential ϕ
E=Xϕ,inΩ0
(3)
ϕ=ϕg,onΩ0g,ϕ
(4)
where the boundary condition in Eq. (4) represents the prescribed electric potential on the boundary part Ω0g,ϕ to Ω0.
In absence of free volume charge density in the dielectric material, Gauss’s law in terms of the electric displacement D is formulated as
XD=0,inΩ0
(5)
DN=DN,onΩ0h,ϕ
(6)
where the boundary condition in Eq. (6) represents the prescribed electric charge on the boundary part Ω0h,ϕ.
The boundary for the balance of linear momentum for quasi-static electro-mechanics is similarly split such that Ω0=Ω0h,φΩ0g,φ and Ω0h,φΩ0g,φ=. This yields
XT+ρ0f=0,inΩ0
(7)
TN=TN,onΩ0h,φ
(8)
φ=φg,onΩ0g,φ
(9)
where f is the mechanical body load per unit volume and T the total first Piola-type stress tensor. The boundary condition in Eq. (8) represents the prescribed traction on Ω0h,φ and Eq. (9) the prescribed displacement on Ω0g,φ.

2.2 Constitutive Model.

The dielectric EAP model used in this work is described in detail by Ask et al. [2729] where also viscoelasticity is included. As demonstrated by Dorfmann and Ogden [26], an augmented, free energy function Ω(F,E) can be introduced such that
Ω=Ωm+Ωmel+Ωel
(10)
where the augmented free energy function is assumed to be separated into a purely mechanical Ωm, an electro-mechanical Ωmel, and a purely electric Ωel part. The purely mechanical response is modeled by neo-Hookean hyper-elasticity
Ωm=12K(J1)2+12G(C¯:I3)
(11)
where K and G are the bulk and shear modulus, respectively, and C¯=J2/3C. The electro-mechanical coupling part Ωmel is taken as
Ωmel=12ε0εrJEC1E
(12)
where ε08.8541012 F/m denotes the vacuum permittivity and εr the constant relative permittivity. The purely electric part Ωel is taken as
Ωel=ceEE
(13)
From Eq. (10), the total Piola stress tensor T and material electric displacements D can be obtained as
D=ΩE=ε0εrJEC12ceE
(14)
and
T=ΩF=K(J1)JFT+GJ2/3[F13[C:I]FT]+ε0εrJ[(EF1)(EC1)12(EC1E)FT]
(15)

2.3 Variational Forms and Finite Element Discretization.

The variational form of the electro-quasi-static problem in Eq. (5) in the material configuration Ω0 is
δWelt=δWexteltδWintelt=0
(16)
for all δϕ that vanishes on Ω0g,ϕ. In Eq. (16)
δWextelt=Ω0h,ϕδϕDNdA
(17)
δWintelt=Ω0XδϕDdV
(18)
Starting from the electro-mechanical balance of linear momentum in Eq. (7) and assuming no mechanical body forces, the variational form in the material configuration Ω0 is
δWmec=δWextmecδWintmec=0
(19)
for all δφ that vanishes on Ω0g,φ. In Eq. (19)
δWextmec=Ω0h,φδφTNdA
(20)
δWintmec=Ω0Xδφ:TdV
(21)

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

The electric potential only requires one additional degrees-of-freedom (DOF) in addition to the displacement degrees-of-freedom. The same shape functions are used to interpolate the displacements φ and potential ϕ
φφh=Nφaφϕϕh=Nϕaϕ
(22)
where Nφ and Nϕ are the global shape functions, aφ and aϕ are corresponding nodal values. Based on the Galerkin approach, the virtual displacements and electric potential are approximated as
δφNφδaφ,δϕNϕδaϕ
(23)
where δaφ and δaϕ are the virtual element nodal values.
Inserting the finite element approximations in Eqs. (22) and (23) into Eqs. (16) and (19) and making use of the arbitrariness of δaφ and δaϕ render the residual equations in matrix format
rϕ=Ω0h,ϕNϕ,TDNdAΩ0Bϕ,TDdV=0rφ=Ω0h,φNφ,TTNdAΩ0Bφ,TTdV=0
(24)
where Bφ and Bϕ holds derivatives of the global shape functions for displacements and potential.
To find a solution to the residuals rϕ=0 and rφ=0, the Newton method is adopted which is based on linearizations of Eqs. (16) and (19)
0=δWelt+dδWelt0=δWmec+dδWmec
(25)
Assuming dead loading condition, the linearizations are found as
dδWelt=dδWintelt=Ω0XδϕDeltXdϕdVΩ0XδϕDmix:XdφdVdδWmec=dδWintmec=Ω0Xδφ:Dmix, TXdϕdV+Ω0Xδφ:Dmec:XdφdV
(26)
where the tangents are defined as
Delt=2ΩEE,Dmix=2ΩEFDmix, T=2ΩFE,Dmec=2ΩFF
(27)
Note that Delt is second-order, Dmix third-order, and Dmec fourth-order tensors.
The finite element discretization of Eqs. (25) and (26) provides the linearized system of equations can now be summarized in matrix format as
KkΔak=rk,[KφφKφϕKϕφKϕϕ][ΔaφΔaϕ]=[rφrϕ]
(28)
which provides the Newton update
ak+1ak+Δak
(29)
that continues until the residual rk is sufficiently small.
The stiffness matrices in matrix format are provided by
Kϕϕ=Ω0Bϕ,TDeltBϕdVKφϕT=Kϕφ=Ω0Bϕ,TDmixBφdVKφφ=Ω0Bφ,TDmecBφdV
(30)
where Dmec, Dmix, and Delt are the Voigt representation of the material tangents in Eq. (27).

3 Topology Optimization Formulation

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

3.1 Filter, Projection, and Material Interpolation.

The optimization design field is regularized by smoothing the piecewise uniform densities ρ such that a continuous ρ~ is obtained. The smoothing is performed using the partial differential equation (PDE) filter by Lazarov and Sigmund [30]
{l02ΔXρ~+ρ~=ρ,in ΩDDXρ~N=0,on ΩDD
(31)
where the filter length scale parameter l0 controls the level of smoothing.
The PDE filter in Eq. (31) is discretized using the same finite element discretization as in the state problem, which render
Kfρ~=Tfρ
(32)
where Kf and Tf can be calculated a priori.
To reduce the amount of gray elements, the smooth Heaviside projection introduced by Wang et al. [31]
ρ¯=H(ρ~)=tanh(βη)+tanh(β(ρ~η))tanh(βη)+tanh(β(1η))
(33)
is used. In Eq. (33), β and η control the sharpness and the projection threshold, respectively. The projection in Eq. (33) was applied on the filtered density field ρ~ in every Gauss point.
Solid isotropic material with penalization (SIMP), see Bendsøe and Sigmund [32], is used to interpolate the terms in Eq. (10) as
ΩSIMP=(δm+(1δm)ρ¯pm)Ωm+(δmel+(1δmel)ρ¯pmel)Ωmel+(δel+(1δel)ρ¯pel)Ωel
(34)
where δm, δmel, and δel are small positive factors that prevent the system matrices from being singular and pm, pmel, and pel are the penalization exponents.

3.2 Optimization Formulation.

The considered optimization problem is formulated as
(TO){minρg0s.t{gi0,i=1,,nconstr0ρe1,e[1,,nelm]
(35)
where the objective function in all examples is
g0=lTa
(36)
where l is a constant vector with nonzero entries in the DOF which should be optimized, a is the solution vector at the terminal load step, and ρe is the element density in element e.
To limit the volume V of the solid region, a volume constraint is enforced
g1=VαVDD10
(37)
where VDD is the volume of the design domain ΩDD, α the allowed volume fraction, and where
V=ΩDDρ¯(ρ)dV
(38)
Remaining constraint functions are described in Sec. 3.3.

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.

Fig. 1
Illustrative example of the VTM where arrows indicate heat flow. Top figure shows two disconnected solid regions. One that is connected to the Dirichlet part of the boundary and will therefore have a lower temperature and one isolated which therefore will have a high temperature. In the bottom figure, the two regions are connected to each other as well as to the Dirichlet part of the boundary. Consequently, the whole solid region will have a low temperature.
Fig. 1
Illustrative example of the VTM where arrows indicate heat flow. Top figure shows two disconnected solid regions. One that is connected to the Dirichlet part of the boundary and will therefore have a lower temperature and one isolated which therefore will have a high temperature. In the bottom figure, the two regions are connected to each other as well as to the Dirichlet part of the boundary. Consequently, the whole solid region will have a low temperature.
Close modal

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, ΩeliTg. The rest of the boundary, Ω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.

The virtual thermal conduction problem
{X(k(ρ¯)XT(X))+Q(ρ¯)=0 for XΩeliTT(X)=0 for XΩeliTgXT(X)N=0 for XΩeliTh
(39)
is solved as a 2D problem in the undeformed configuration in each electrode ΩeliT and will give the temperature distribution T(X). Connectivity is now enforced by constraining the maximum allowed temperature as
g2=TpλTref0
(40)
where λ is a dimensionless factor, the maximum temperature Tp in the solid domain is approximated using the smooth p-norm of the temperatures, i.e.,
Tp=(ΩeliT(ρ¯T)pdA)1/p
(41)
The scaling factor ρ¯ in Eq. (41) removes high virtual temperatures in void regions. Tref is a reference peak temperature calculated with Eq. (41) when ρ=1 in the whole electrode.
The SIMP formulation [32] is used to interpolate both the thermal conductivity coefficient k and the heat supply Q in Eq. (39) as
k(ρ¯)=(ϵ+ρ¯pk(1ϵ))k0,Q(ρ¯)=ρ¯pQQ0
(42)
where k0=1 W/m/K and Q0=1 W/m3 are arbitrary since they are formulated in relation to a reference temperature. The value of ϵ should be small but sufficiently large to ensure that the stiffness matrix does not become singular. In this work, ϵ=105 was used.
The virtual temperature problem is solved using conventional FEM, i.e.,
Kc(ρ¯)ac=fc(ρ¯)
(43)

3.4 Sensitivity Analysis.

The design updates are obtained using the gradient-based method of moving asymptotes (MMA), see Svanberg [33], which uses gradients of the objective and constraint functions with respect to the design variables ρ. Because a PDE filter is used, the chain rule is applied to the objective and constraint functions such that
giρ=giρ~ρ~ρ
(44)
where ρ~ρ is obtained from the filter in Eq. (31).
The sensitivities of the volume constraint g1=VαVDD1 with respect to the filtered design variables are calculated using
Vρ~=ΩDDρ¯ρ~dV=ΩDDρ¯ρ~ρ~ρ~dV=ΩDDHNTdV
(45)
where Eq. (33) was used.
The sensitivity of the objective function g0(a)=lTa is calculated with the adjoint method wherein the objective function is augmented with λTr(a,ρ~)=0, where λ is the adjoint and rT=[rφ,Trϕ,T] the residual. The sensitivities with respect to the filtered densities ρ~ is now computed as
g~0(a,ρ~)ρ~=g~0ρ~+g~0aaρ~=λTrρ~+(lTλTK)aρ~
(46)
To annihilate the implicit sensitivity aρ~, λ is assigned to fulfill
Kλ=l
(47)
where K is the stiffness matrix in the last load step. Note that a homogeneous Dirichlet boundary condition is enforced in the DOF where displacement or potential is prescribed.
Once Eq. (47) is solved, the sensitivity can be computed as
g~0(a,ρ~)ρ~=λTrρ~
(48)
To compute the sensitivity of the connectivity constraint g2 in Eq. (40), the adjoint approach is used again but now replacing λTr with λcTrc where rc=Kcacf. Using identical approach as in Eqs. (46)(48) results in
g~2ρ~=g2ρ~λcTrcρ~
(49)
where
Kcλc=(g2ac)T
(50)
Note that the sensitivities of the connectivity constraint will only be nonzero for the nodes corresponding to the electrode.

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 109 to reduce the influence of the void region, while avoiding a singular stiffness matrix.

Table 1

Material properties for the EAP material

K (MPa)G (MPa)ce (NV-2)εr()
0.60.114.7
K (MPa)G (MPa)ce (NV-2)εr()
0.60.114.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.

Table 2

Penalization exponents

pmpmelpel
363
pmpmelpel
363

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 Δϕ 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 α=40% is used. All finite element models are based on 3D solid elements.

Fig. 2
(a) Illustration and geometry of the cylinder, which is discretized using 3D solid elements. Symmetry with respect to x=0 and z = 0 is assumed. To prevent rigid body motion, the structure is properly constrained. (b) The solution domain for the virtual temperature problem. The boundary conditions used for the VTM problem in both electrodes are homogeneous Dirichlet on the left side and homogeneous Neumann on the remaining sides, but with different radius r.
Fig. 2
(a) Illustration and geometry of the cylinder, which is discretized using 3D solid elements. Symmetry with respect to x=0 and z = 0 is assumed. To prevent rigid body motion, the structure is properly constrained. (b) The solution domain for the virtual temperature problem. The boundary conditions used for the VTM problem in both electrodes are homogeneous Dirichlet on the left side and homogeneous Neumann on the remaining sides, but with different radius r.
Close modal
Table 3

Geometrical parameters for cylinder shell in Fig. 2 

Case123
r (mm)404040
lz (mm)4080120
t1 (mm)111
t2 (mm)111
Δϕ (kV)151515
Elements4×40×1324×80×1324×120×132
Case123
r (mm)404040
lz (mm)4080120
t1 (mm)111
t2 (mm)111
Δϕ (kV)151515
Elements4×40×1324×80×1324×120×132

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,

l={1nyDOF on top edge1nyDOF on bottom edge0else
(51)
where n is the total number of nodes on the both edges and the sign accounts for in which direction the displacement should be minimized.

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 β was initiated to β=1 and then increased by 1 every 20th iteration to a maximum of 15. The parameter η=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 |g0kg0k1g0k|<105, and until all constraints were fulfilled.

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

Optimization algorithm without connectivity

Algorithm 1

Set initial values on all parameters

Set ρ for the initial design

While|g0kg0k1|>105|g0k| and continuation terminated do

  If applicable, update parameters β and pk,

  Apply filter to ρ to obtain ρ~ from Eq. (32)

  Project filtered densities ρ~ to obtain ρ¯ from Eq. (33)

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

  Calculate constraints gi, i=1,,nconstr

  Compute sensitivities using the adjoint method

  Solve MMA to update element densities ρ using MMA

end While

4.1 No Connectivity Constraint.

In the first example, only a volume constraint is introduced, i.e.,
(TO)1{minρg0=uout=lTas.t{g1=VαVDD100ρe1,e[1,,nelm]
(52)
where s.t. means subjected to.

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

Fig. 3
Projected densities ρ¯ in the design domain for the final optimized structures described in Table 3 without connectivity constraints. As a result, two discontinuous, finger-like islands are formed on top and bottom sides: (a) lz=40 mm, g0=−7.46 mm, (b) lz=80 mm, g0=−7.56 mm, (c) lz=120 mm, g0=−7.51 mm, and (d) convergence of objective function.
Fig. 3
Projected densities ρ¯ in the design domain for the final optimized structures described in Table 3 without connectivity constraints. As a result, two discontinuous, finger-like islands are formed on top and bottom sides: (a) lz=40 mm, g0=−7.46 mm, (b) lz=80 mm, g0=−7.56 mm, (c) lz=120 mm, g0=−7.51 mm, and (d) convergence of objective function.
Close modal

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 α=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.

Fig. 4
Deformed state for the whole optimized structure for case 1 in Table 3 compared with a structure with straight electrode parts fulfilling the volume constraint. Element densities with ρ<0.0001 have been removed and the electric potential ϕ is plotted on the remaining elements. Note that the optimized structure has about 3% larger absolute value of the objective function g0 as compared to the straight case: (a) optimized structure g0=−7.46 mm and (b) straight structure g0=−7.26 mm.
Fig. 4
Deformed state for the whole optimized structure for case 1 in Table 3 compared with a structure with straight electrode parts fulfilling the volume constraint. Element densities with ρ<0.0001 have been removed and the electric potential ϕ is plotted on the remaining elements. Note that the optimized structure has about 3% larger absolute value of the objective function g0 as compared to the straight case: (a) optimized structure g0=−7.46 mm and (b) straight structure g0=−7.26 mm.
Close modal

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.

Adding the connectivity constraints described in Eq. (40) for both electrodes, see Fig. 2 for the virtual temperature solution domain, yields two extra constraints, g2 and g3 for the inner and outer electrodes respectively, i.e.,
(TO)2{minρg0=uout=lTas.t{g1=VαVDD10g2=TpTrefλmax0g3=TpTrefλmax00ρe1,e[1,,nelm]
(53)

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 λmax=10.

When enforcing the connectivity constraint in Eq. (40), the optimization prefers in general thin connections which could cause problems in the design updates if too much material are removed from them, potentially causing a reduction in connectivity. To combat this and stabilize the optimization procedure, an additional constraint is enforced in each electrode
λminTrefTp0
(54)
where a lower bound on the temperature Tp is enforced. In the numerical examples, this constraint is inactive in the final designs, meaning it does not constrain them, but it stabilizes the solution procedure. The lower bound started at λmin=0 to allow the constraint to be inactive initially and then increased at iteration 40 to λmin=λmax0.01, close to the upper bound.

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.

Fig. 5
Projected densities ρ¯ in the deformed design domain for the final optimized structures described in Table 3 with connectivity constraints. A comparison between the Fig. 5 designs with existing solid phase connection to the Fig. 3 designs without connections shows that enforcing connectivity decrease the performance by about 6–9%. Similar finger-like structures are also still formed: (a) lz=40 mm, g0=−6.85 mm, (b) lz=80 mm, g0=−6.90 mm, (c) lz=120 mm, g0=−7.07 mm, and (d) convergence of objective function.
Fig. 5
Projected densities ρ¯ in the deformed design domain for the final optimized structures described in Table 3 with connectivity constraints. A comparison between the Fig. 5 designs with existing solid phase connection to the Fig. 3 designs without connections shows that enforcing connectivity decrease the performance by about 6–9%. Similar finger-like structures are also still formed: (a) lz=40 mm, g0=−6.85 mm, (b) lz=80 mm, g0=−6.90 mm, (c) lz=120 mm, g0=−7.07 mm, and (d) convergence of objective function.
Close modal

Since the upper limit for the virtual temperature constraint, λ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 λmax{6,10,14} is shown. In this comparison, λmin is initiated to 0 and then changed to λmin=λmax0.01 at iteration 40. In all cases, a clear connection is formed between the upper and lower parts. The material distribution is similar for λmax=6 and λmax=10 while λmax=14 gives a noticeable difference. The fingers previously shown in Figs. 3 and 5 are again present.

Fig. 6
The projected densities ρ¯ in the deformed design domain for different choices for λmax and where λmin=λmax−0.01 for case 1 in Table 3 at the last iteration: (a) λmax=6, g0=−6.68 mm, (b) λmax=10, g0=−6.85 mm, (c) λmax=14, g0=−6.91 mm, and (d) convergence objective function
Fig. 6
The projected densities ρ¯ in the deformed design domain for different choices for λmax and where λmin=λmax−0.01 for case 1 in Table 3 at the last iteration: (a) λmax=6, g0=−6.68 mm, (b) λmax=10, g0=−6.85 mm, (c) λmax=14, g0=−6.91 mm, and (d) convergence objective function
Close modal

This numerical experiment shows that the choice of λ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 λmax, in general, results in lower, more negative g0, meaning larger displacements. This is expected. Higher λ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/m3)

V =

volume (m3)

B =

derivatives of shape functions

C =

Cauchy–Green tensor

ce =

electric parameter (N/V2)

D =

electric displacement (C/m2)

DN =

prescribed electric charge on boundary (C/m2)

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/m3)

Tp =

approximate maximum temperature (K)

Tref =

reference temperature (K)

VDD =

volume in design domain (m3)

TN =

prescribed traction

Greek Symbols

α =

allowed volume fraction

β =

sharpness parameter

ε0 =

vacuum permittivity (F/m)

εr =

relative permittivity

η =

threshold parameter

λ =

adjoint vector

λc =

adjoint vector for VTM

λ =

factor in VTM constraint

λmax =

maximum factor in VTM constraint

λmin =

minimum factor in VTM constraint

ρ =

element density

ρ~ =

filtered element density

ρ¯ =

projected element density

ϕ =

scalar electric potential (V)

ϕg =

prescribed scalar electric potential (V)

φ =

deformation mapping (m)

φg =

prescribed displacement (m)

Ω =

augmented free energy function

Ωm =

mechanical free energy function

Ωmel =

electro-mechanical free energy function

Ωel =

electric free energy function

Ω0 =

material configuration

Ω0 =

boundary to material configuration

Ωt =

spatial configuration

Ωt =

boundary to spatial configuration

Superscripts and Subscripts

m =

mechanical terms

t =

spatial configuration

el =

electric terms

mel =

electro-mechanical terms

0 =

material configuration, objective function

ϕ =

potential component

φ =

displacement component

Acronyms and Abbreviations

EAP =

electroactive polymer

FEM =

finite element method

MMA =

method of moving asymptotes

SIMP =

solid isotropic material with penalization

TO =

topology optimization

VTM =

virtual temperature method

References

1.
Bar-Cohen
,
Y.
, ed.,
2004
,
Electroactive Polymer (EAP) Actuators as Artificial Muscles: Reality, Potential, and Challenges
, 2nd ed.,
SPIE Press
,
Bellingham, WA
.
2.
Granlund
,
G.
,
Wallin
,
M.
,
Günther-Hanssen
,
O.
,
Tortorelli
,
D.
, and
Watts
,
S.
,
2024
, “
Topology Optimization of Compliant Mechanisms Under Transient Thermal Conditions
,”
Comput. Methods Appl. Mech. Eng.
,
418
, p.
116478
.
3.
Zhu
,
X.
,
Zhao
,
C.
,
Wang
,
X.
,
Zhou
,
Y.
,
Hu
,
P.
, and
Ma
,
Z.-D.
,
2019
, “
Temperature-Constrained Topology Optimization of Thermo-Mechanical Coupled Problems
,”
Eng. Optim.
,
51
(
10
), pp.
1687
1709
.
4.
Xu
,
X.
,
Wu
,
Y.
,
Zuo
,
L.
, and
Chen
,
S.
,
2021
, “
Topology Optimization of Multimaterial Thermoelectric Structures
,”
ASME J. Mech. Des.
,
143
(
1
), p.
011705
.
5.
Xing
,
J.
,
Luo
,
Y.
,
Deng
,
Y.
,
Wu
,
S.
, and
Gai
,
Y.
,
2023
, “
Topology Optimization Design of Deformable Flexible Thermoelectric Devices for Voltage Enhancement
,”
Eng. Optim.
,
55
(
10
), pp.
1686
1703
.
6.
Lee
,
S.
, and
Tovar
,
A.
,
2013
, “
Topology Optimization of Piezoelectric Energy Harvesting Skin Using Hybrid Cellular Automata
,”
ASME J. Mech. Des.
,
135
(
3
), p.
031001
.
7.
Ortigosa
,
R.
,
Martínez-Frutos
,
J.
,
Ruiz
,
D.
,
Donoso
,
A.
, and
Bellido
,
J. C.
,
2021
, “
Density-Based Topology Optimisation Considering Nonlinear Electromechanics
,”
Struct. Multidiscipl. Optim.
,
64
(
1
), pp.
257
280
.
8.
Ortigosa
,
R.
, and
Martínez-Frutos
,
J.
,
2021
, “
Multi-Resolution Methods for the Topology Optimization of Nonlinear Electro-Active Polymers at Large Strains
,”
Comput. Mech.
,
68
(
2
), pp.
271
293
.
9.
Martínez-Frutos
,
J.
,
Ortigosa
,
R.
, and
Gil
,
A. J.
,
2021
, “
In-Silico Design of Electrode Meso-Architecture for Shape Morphing Dielectric Elastomers
,”
J. Mech. Phys. Solids
,
157
, p.
104594
.
10.
Ortigosa
,
R.
,
Martínez-Frutos
,
J.
, and
Gil
,
A. J.
,
2023
, “
Programming Shape-Morphing Electroactive Polymers Through Multi-Material Topology Optimisation
,”
Appl. Math. Model.
,
118
, pp.
346
369
.
11.
Bortot
,
E.
,
Amir
,
O.
, and
Shmuel
,
G.
,
2018
, “
Topology Optimization of Dielectric Elastomers for Wide Tunable Band Gaps
,”
Int. J. Solids Struct.
,
143
, pp.
262
273
.
12.
Sharma
,
A. K.
,
Kosta
,
M.
,
Shmuel
,
G.
, and
Amir
,
O.
,
2022
, “
Gradient-Based Topology Optimization of Soft Dielectrics as Tunable Phononic Crystals
,”
Compos. Struct.
,
280
, p.
114846
.
13.
Dev
,
C.
,
Stankiewicz
,
G.
, and
Steinmann
,
P.
,
2023
, “
On the Influence of Free Space in Topology Optimization of Electro-Active Polymers
,”
Struct. Multidiscipl. Optim.
,
66
(
8
), p.
187
.
14.
Liu
,
S.
,
Li
,
Q.
,
Chen
,
W.
,
Tong
,
L.
, and
Cheng
,
G.
,
2015
, “
An Identification Method for Enclosed Voids Restriction in Manufacturability Design for Additive Manufacturing Structures
,”
Front. Mech. Eng.
,
10
(
2
), pp.
126
137
.
15.
Li
,
Q.
,
Chen
,
W.
,
Liu
,
S.
, and
Tong
,
L.
,
2016
, “
Structural Topology Optimization Considering Connectivity Constraint
,”
Struct. Multidiscipl. Optim.
,
54
(
4
), pp.
971
984
.
16.
Swartz
,
K. E.
,
Tortorelli
,
D. A.
,
White
,
D. A.
, and
James
,
K. A.
,
2022
, “
Manufacturing and Stiffness Constraints for Topology Optimized Periodic Structures
,”
Struct. Multidiscipl. Optim.
,
65
(
4
), p.
129
.
17.
Donoso
,
A.
, and
Guest
,
J. K.
,
2019
, “
Topology Optimization of Piezo Modal Transducers Considering Electrode Connectivity Constraints
,”
Comput. Methods Appl. Mech. Eng.
,
356
, pp.
101
115
.
18.
Luo
,
Y.
,
Sigmund
,
O.
,
Li
,
Q.
, and
Liu
,
S.
,
2020
, “
Additive Manufacturing Oriented Topology Optimization of Structures With Self-Supported Enclosed Voids
,”
Comput. Methods Appl. Mech. Eng.
,
372
, p.
113385
.
19.
Luo
,
Y.
,
Sigmund
,
O.
,
Li
,
Q.
, and
Liu
,
S.
,
2022
, “
Topology Optimization of Structures With Infill-Supported Enclosed Voids for Additive Manufacturing
,”
Addit. Manuf.
,
55
, p.
102795
.
20.
Donoso
,
A.
,
Aranda
,
E.
, and
Ruiz
,
D.
,
2022
, “
A New Approach Based on Spectral Graph Theory to Avoiding Enclosed Holes in Topology Optimization
,”
Comput. Methods Appl. Mech. Eng.
,
393
, p.
114769
.
21.
Donoso
,
A.
,
Aranda
,
E.
, and
Ruiz
,
D.
,
2023
, “
A New Method for Designing Piezo Transducers With Connected Two-Phase Electrode
,”
Comput. Struct.
,
275
, p.
106936
.
22.
Donoso
,
A.
,
Aranda
,
E.
, and
Ruiz
,
D.
,
2023
, “
A Continuous Model for Connectivity Constraints in Topology Optimization
,”
Struct. Multidiscipl. Optim.
,
66
(
4
), p.
71
.
23.
Jackson
,
J. D.
,
1962
,
Classical Electrodynamics
,
Wiley
,
New York
.
24.
Eringen
,
A. C.
, and
Maugin
,
G. A.
,
1990
,
Electrodynamics of Continua. 1: Foundations and Solid Media
,
Springer
,
New York Heidelberg
.
25.
Kovetz
,
A.
,
2000
,
Electromagnetic Theory
,
Oxford Science Publications, Oxford University Press
,
Oxford; New York
.
26.
Dorfmann
,
A.
, and
Ogden
,
R. W.
,
2005
, “
Nonlinear Electroelasticity
,”
Acta Mech.
,
174
(
3
), pp.
167
183
.
27.
Ask
,
A.
,
Menzel
,
A.
, and
Ristinmaa
,
M.
,
2012
, “
Phenomenological Modeling of Viscous Electrostrictive Polymers
,”
Int. J. Non-Linear Mech.
,
47
(
2
), pp.
156
165
.
28.
Ask
,
A.
,
Menzel
,
A.
, and
Ristinmaa
,
M.
,
2012
, “
Electrostriction in Electro-Viscoelastic Polymers
,”
Mech. Mater.
,
50
, pp.
9
21
.
29.
Ask
,
A.
,
Menzel
,
A.
, and
Ristinmaa
,
M.
,
2015
, “
Modelling of Viscoelastic Dielectric Elastomers With Deformation Dependent Electric Properties
,”
Procedia IUTAM
,
12
, pp.
134
144
.
30.
Lazarov
,
B. S.
, and
Sigmund
,
O.
,
2011
, “
Filters in Topology Optimization Based on Helmholtz-Type Differential Equations
,”
Int. J. Numer. Methods Eng.
,
86
(
6
), pp.
765
781
.
31.
Wang
,
F.
,
Lazarov
,
B. S.
, and
Sigmund
,
O.
,
2011
, “
On Projection Methods, Convergence and Robust Formulations in Topology Optimization
,”
Struct. Multidiscipl. Optim.
,
43
(
6
), pp.
767
784
.
32.
Bendsøe
,
M. P.
, and
Sigmund
,
O.
,
1999
, “
Material Interpolation Schemes in Topology Optimization
,”
Archive Appl. Mech.
,
69
(
9
), pp.
635
654
.
33.
Svanberg
,
K.
,
1987
, “
The Method of Moving Asymptotes—A New Method for Structural Optimization
,”
Int. J. Numer. Methods Eng.
,
24
(
2
), pp.
359
373
.