On the influence of free space in topology optimization of electro-active polymers

This study investigates the impact of the surrounding free space on the topology optimization (TO) of electro-active polymers (EAPs). It is well understood that, under the application of an electric field, the deformation of an EAP is not solely determined by the field distribution within the body, but also by the distribution in the free space surrounding it. This is particularly true for electronic EAP, which are emerging as leading candidates for developing artificial muscles. Our study specifically focuses on understanding the influence of the free space in the context of density-based TO. We model the free space as an extended void region around the design domain. Our numerical experiments focus on EAP actuators and take into account their geometrical nonlinear behavior. The results show that incorporating the surrounding free space has a significant impact on the performance of the optimized EAPs with low electric permittivity. This makes it essential to consider in real-world applications.


Introduction
Electro-active polymers (EAP) are a class of smart materials that react to electric stimulation with bulk deformation. EAPs are used in soft robots (Lee et al. 2017;Gu et al. 2021), artificial muscles (Bar-Cohen 2002;Carpi et al. 2011), and in a broad range of industries (Bashir and Rajendran 2018).
Numerous examples of finite-strain-coupled numerical models can be found in the literature to simulate EAPs (Yang and Batra 1995;Dorfmann and Ogden 2005;Vu et al. 2007;Zwecker et al. 2011;Ask et al. 2012;Hossain et al. 2012;Skatulla et al. 2012;Büschel et al. 2013;Gil and Ortigosa 2016). In all of the aforementioned works on EAPs, only their material bodies were considered and discretized using the finite element method. This is a reasonable approach for simulating EAPs with higher permittivity, e.g. piezoelectric polymers.
However, in the case of electronic EAPs, where the electric permittivity is approximately one order of magnitude greater than that of vacuum, the electric energy stored in the surrounding free space also has a significant influence on the behavior of EAP. This influence was studied by Steinmann (2010, 2012) and Steinmann (2011), by considering the free space through a coupled BEM-FEM method. In Pelteret et al. (2016), a mixed variational formulation for quasi-incompressible electro-active or magneto-active polymers, which accounts for the influence of the surrounding free space was studied. In Steinmann and Vu (2017) the computational challenges in the simulation of EAPs were presented and the authors concluded that to build a complete picture of what happens inside an EAP, what happens outside deserves due attention.
In the recent decades, topology optimization (TO) has been increasingly utilized to generate designs based on mathematical optimization techniques that go beyond human intuition (Sigmund and Maute 2013;Deaton and Grandhi 2014). TO has been used to optimize electromechanical actuators for MEMS application Qian and Sigmund (2013), as well as dielectric elastomer actuators by Wang et al. (2019) and Chen et al. (2020). In the recent work of Ortigosa et al. (2021), TO was applied to optimize EAPs. However, in all the aforementioned TO problems, the electric field in the surrounding free space was disregarded. This approach is reasonable considering the materials used and the focus on optimizing thin film EAPs with an electric field applied perpendicular to the film. Nevertheless, there is still a lack of research on electronic EAPs where the electric energy stored in the surrounding free space also plays a significant role. In this study, we specifically consider this class of material. Additionally, unlike the previous works that concentrated on thin films, our investigation of 2D structures does not involve any dimension that is significantly smaller than the others. These distinctions allow us to explore the influence of free space on the optimized structure.
The layout of this paper is as follows: Sect. 2 introduces the essence of nonlinear electro-elasticity with free space. Section 3 describes TO for the electro-elastic problem. Finally, Sect. 4 illustrates numerical examples where topology optimization of EAP is carried out with the consideration of free space. Finally, Sect. 5 provides some concluding remarks.

Nonlinear electro-elasticity with freespace
This section presents a brief introduction to the formulation and computation of geometrically nonlinear electro-elasticity immersed in the surrounding free space.

Kinematics and state equations
Consider an EAP B 0 immersed in the surrounding free space (vacuum) Ω 0 in their material configuration at time t 0 . Subsequently, due to deformation, the EAP and the free space will occupy a spatial configuration B t and Ω t at time t ( t > t 0 ), as shown in Fig. 1. The transformation of the EAP into the spatial configuration is defined by a deformation map . The deformation map is associated with the deformation gradient F , which is defined as follows: Here, ∇ X is the gradient operator with respect to material coordinates X . The deformation gradient F is associated with the Jacobian J and the right Cauchy-Green strain tensor C , which are defined as follows: We can characterize the quasi-static elastic problem through the balance of linear momentum, which can be written as follows: (1) F = ∇ X .
Here, the operator Div represents divergence with respect to material coordinates X . The jump, denoted by [[•]] , is defined as the difference in the value of (•) when transitioning from inside the material to the free space outside The total Piola stress is represented by the symbol P , while f 0 denotes the mechanical body force per unit undeformed volume in B 0 . The term t 0 represents the applied mechanical traction force per unit undeformed area on Γ t 0 ⊂ Γ 0 . The prescribed displacement on Γ 0 ⊂ Γ 0 , is denoted by . Further, Γ t 0 ∪ Γ 0 = Γ 0 and Γ t 0 ∩ Γ 0 = � . We can characterize the electrical problem through the electrostatic Maxwell equations and restrict ourselves to materials without free currents and free electric charges, which leads to: where D is the electric displacement vector in Lagrangian form and 0 is the free surface charge density on Γ 0 ⊂ Γ 0 . Furthermore, Faraday's law results in where E is the electric field vector in Lagrangian form, is the scalar electric potential and is the prescribed electric potential on Γ 0 ⊂ Γ 0 , with Γ 0 ∪ Γ 0 = Γ 0 and Γ 0 ∩ Γ 0 = �.

Constitutive modeling
Next, to solve the state problem described by Eqs. (3)-(5), a constitutive relation is required. To this end, we define an appropriate total energy density functional = (F, E) in terms of the unknown deformation map and electric potential . The total energy density function for the bulk material, denoted as b (F, E) , is defined as: where and are the two Lamé coefficients of the classical neo-Hookean material law and n dim is the dimension of the space. The material constants c 1 and c 2 are responsible for electromechanical coupling and r is the relative electric permittivity of the bulk material.
To model the influence of the free space, we define an appropriate energy functional in the free space as: To ensure that we can interpolate the energy functional of the bulk material and the free space in TO, we must consider pseudoelastic properties in the functional above. However, to minimize the contribution of these properties, we introduce a dimensionless coefficient of the order of 10 −15 , denoted by . In the equation above, C 0 represents the fourth-order linear elasticity tensor, while 0 = 8.854 × 10 −12 F m is the electric permittivity of free space.
Furthermore, Ortigosa et al. (2021) noted that the electromechanical component of the energy functional in Eq. 7 compromises numerical stability due to the high non-convexity of the functional with respect to the deformation gradient tensor F . Therefore, we adopt the stable version of the energy functional proposed in Ortigosa et al. (2021). This version eliminates the dependence of the electroelastic component on the deformation gradient, resulting in the following definition: We acknowledge that the simplification in Eq. 8 does not accurately account for the Maxwell stress. However, we conducted a comparison of the electric energy density defined in Eqs. 7 and 8 for our specific problem setting in Sect. 2.4. Our comparison revealed that the difference between the two formulations is not significant in our problem setup. Therefore, it is reasonable to utilize the numerically stable version of the energy function as defined in Eq. 8. (8) Having defined the energy density functions we can now define the total Piola stress P and the referential dielectric displacement D as Furthermore, the second derivatives of the energy density functional yields the constitutive tangent tensor, namely, the fourth-order elasticity tensor C , the third-order piezoelectric tensor P , and the second-order dielectric tensor , defined respectively as

Variational formulation in nonlinear electro-elasticity
Here we discuss the variational approach to derive the weak form of the governing equations, (3) and (4). To this end, we define the total potential energy functional as follows: The stationary point min max Π ⇒ Π = 0 defines the equilibrium solution of the system. To determine the stationary point, we take the variation of the energy functional with respect to variations of the displacement field, resulting in: Similarly, we take the variation of the energy functional with respect to variations of the electric potential, resulting in (9) P = F and D = − E . (10) Here, P b and P f represent the total Piola stress obtained from the energy density functions of the bulk and free space, respectively, as defined in Eq. 9. Similarly, D b and D f represent the referential dielectric displacement of the bulk and free space as defined in Eq. 9. In order to solve for the stationary condition in Eq. 12 and 13, we utilize a Newton-Raphson scheme where we linearize Π with respect to the incremental displacement field Δ and the incremental electric potential Δ , resulting in the following system of equations: Here, the fourth-order elasticity tensor for the bulk and free space is denoted by C b and C f , respectively. Similarly, the third-order piezoelectric tensor for the bulk and free space is denoted by P b and P f , and the second-order dielectric tensor for the bulk and free space is denoted by b and f , respectively. The next step is to solve the linearized system, Π( ) + Π( )

Numerical example
In the first example, we validate our formulation using the benchmark problem presented in Vu and Steinmann (2010). We consider a square plate with dimensions of 60 µm × 60 µm, which contains a square hole of dimensions 30 µm × 30 µm. The free space surrounding the plate has dimensions of 120 µm × 120 µm. The lower and upper edges of the plate have prescribed electric potentials of −500V and 500V, respectively. The setup with the inclusion of free space is illustrated in Fig. 2. The constitutive parameters for the energy density function in Eq. 6 are taken from Vu and Steinmann (2010), with = 0.05 MPa, = 0.06 MPa, c 1 = 0.2 0 , c 2 = 2 0 , and r = 5 0 .
To illustrate the significance of the consideration of free space in the simulation, we carried out the simulation in two setups. In the first setup, the bulk material is immersed in free space, and in the second setup, we consider only the bulk material. In Fig. 3, we present the electric potential distribution in the deformed configuration for both setups. The magnitude of the maximum displacement in the setup with the inclusion of free space is 3.1 µm, whereas without free space, it is 2.4 µm, which is 22% less. To get a closer look at the difference, in Fig. 4, we plot the displacement along the top face. This plot clearly shows that considering free space is important, and the results obtained are also in agreement with the results in Vu and Steinmann (2010).
In the second example, we demonstrate the complexities in nonlinear coupled electro-mechanics. We consider a linear guiding actuator setup, as shown in Fig. 5. During the simulation, we consider different sizes of free space truncation, namely, 4 µm, 16 µm, 64 µm, 256 µm, and 512 µm, as illustrated in Fig. 6a. After solving the problems with different sizes of free space truncation, we plot the electric potential along the diagonal from bottom left to top right in Fig. 6b. It is important to note the potential distribution on the right side. The relationship between the size of free space truncation and the potential at the top right of the bulk material is not monotonous. This complex behavior motivated our subsequent study on the influence of free space size on topology optimization. Furthermore, we conducted a comparison between the electric energy densities in free space as defined by Eqs. 8 and 9. To evaluate the differences, we introduced a normalized difference metric given by the equation: For this analysis, we chose a free space truncation of 256 µm. The visualization of the normalized difference in electric energy density is depicted in Fig. 7. Notably, the dissimilarity is predominantly observed near the boundary where the potential is prescribed, with the energy density difference amounting to only 0.1% . Based on this observation, it is justified to employ the stable version of the energy function defined in Eq. 9.

Topology optimization
The goal of TO is to obtain an optimal distribution of EAP material to maximize the displacement at an output port, while satisfying a given volume constraint. Mathematically, the optimization problem is defined as follows: where is the variable pseudo-density (design variable), l is an indicator vector with values of 1 on the degrees of freedom where the displacement has to be maximized, V is the restriction on the allowed volume of the material, e is (17) the pseudo-density corresponding to the e-th finite element, and N e is the number of finite elements (design variables).
Since the pseudo-density varies continuously from void to solid (i.e., from [0, 1]), we interpolate the energy density function according to the well-known SIMP method (Bendsoe and Sigmund 2003). In our case, we interpolate the total energy density functions of the solid (Eq. 6) and free space (Eq. 8) to define the interpolated energy density function SIMP as follows: where p is the penalization parameter, set to 5 in all our numerical studies.
Furthermore, to prevent the appearance of checkerboard patterns and to overcome the dependency on the mesh size in the final design, we have adopted a density filter (Bourdin 2001;Bruns and Tortorelli 2001). The filtered pseudodensity ( ̃ ) is defined as follows: where i is a circular neighborhood centered at the center of element i which contains all the elements whose center are within a radius R.
Due to the use of the radius filter, some intermediate densities may still exist, which can be reduced using a smooth Heaviside projection filter (Wang et al. 2011). Here, controls the sharpness of the Heaviside projection, and is the threshold density at which the transition takes place. The interpolated energy density function in Eq. 18, which utilizes the projected density ( ), is given by: Furthermore, the volume constraint in Eq. 17 must also consider the projected pseudo-density .
In addition, when the design is geometrically symmetric, we apply an additional projection-based symmetric constraint based on the method proposed by Vatanabe et al. (2016). This constraint helps to maintain the symmetry of the optimized solution, which can be affected by finite precision errors in floating-point computations.

Sensitivity analysis and design update
The sensitivity analysis is performed using an adjoint method. The steps to perform adjoint sensitivity analysis involving electro-elastic problems are detailed in Ortigosa et al. (2021). Based on the computed sensitivities, the design is updated using an optimality criteria method (OCM), as described in Bendsoe and Sigmund (2003).

Results
In this section, we present the results of our study on the influence of free space on the optimized design. We provide two benchmark examples: (a) linear guiding actuator and (b) gripper actuator. We increase the size of the free space truncation considered in the computation and observe its impact on the optimized design.
The material parameters used to define the energy functional in Eqs. 6 and 8 are common for all examples and are as follows: = 0.1 MPa, = 0.05 MPa, 0 = 8.854 × 10 −12 F m , c 1 = 0.2 0 , c 2 = 2 0 , r = 5 0 , and = 10 −15 . The optimization parameters are also common for all examples and are listed in Table 1.
In our implementation of the OCM (Bendsoe and Sigmund 2003), we set the lower Lagrange value as 0, and initialize the upper Lagrange value as 1. In subsequent iterations, the upper Lagrange value is set as 10 times the Lagrange value from the previous iteration. In the smooth Heaviside projection filter (Eq. 21), we initialize = 1 and keep it unchanged until iteration 100. After that, we update i = 1.04 × i−1 until we reach the value of indicated in Table 1.

Linear guiding actuator
The first example we consider is a linear guiding actuator. When an electric potential difference of 20 V is applied, the actuator maximizes the displacement at the output port u out in the horizontal direction against a reaction force of f out = 1 N applied by a body being pushed. The problem We discretize the truncated free space with a coarser mesh, as illustrated in Fig. 8. The use of a coarser mesh improves the numerical performance. Furthermore, the finite elements of the free space are not included in the optimization design space, and thus, do not require any special treatment in the topology optimization routine.
To investigate the impact of free space on the optimized design, we conducted TO for the linear guiding actuator while considering different sizes of free space truncation in the computation. The sizes considered were 0.5L = 30 µm, L = 60 mum, 2L = 120 µm, 4L = 240 µm, and 6L = 360 µm, as shown in Fig. 9g.
In Figs. 9a-f, we present the results of TO for different free space truncations. It is observed that the consideration of free space significantly affects the optimized structure. This difference is clearly visualized in Fig. 9h, where we use a shape generation algorithm  to extract the outline of the structure. The optimized structures converge as we increase the free space truncation, with the difference in the optimized structure being minimal for free space sizes of 4 times (240 µm) and 6 times (360 µm) the bulk dimension. This convergence is also evident in the evolution of the objective plot in Fig. 10. Finally, we note that the difference in the objective value for free space size of 0 and 6L (360 µm) is 9.98%.
In Fig. 11, we present the potential distribution along the diagonal of the design domain to further examine how the consideration of different sizes of free space influences the potential distribution inside the design domain (bulk). Here, we observe a convergence in potential distribution for a free space size of 4 times ( 240 m) and 6 times (360 µm) bulk dimension.
Next, we consider the same problem setup as shown in Fig. 5, but on a larger scale, to investigate whether the influence of free space remains consistent for larger domain sizes. To achieve this, we set L = 60 mm and maintain an electric potential difference of 20 kV. Consequently, the truncated free space sizes are now 0.5L = 30 mm, L = 60 mm, 2L = 120 mm, 4L = 240 mm, and 6L = 360 mm, as depicted in Fig. 12g.
The TO results for different free space truncations are presented in Figs. 12a-12f. It is evident that there are no significant differences in the geometry of the designs as the domain size increases. This can be clearly observed by comparing the extracted shape in Fig. 12h with the shape extracted from a smaller domain in Fig. 9h. Furthermore, this similarity is also apparent in the evolution of the objective shown in Fig. 13. Based on these findings, we can conclude that the influence of free space is independent of the domain size.

Gripper actuator
The next example we studied is a gripper actuator. When actuated with an electric potential of 20V, the actuator arms at the output port, u out , move towards each other in the vertical direction against a reaction force f out = 1 N applied by a body being gripped. The problem setup is illustrated in Fig. 14 and L = 60 µm in our study.
We conducted a similar study as the previous example to investigate the influence of free space on the optimized structure. We considered different sizes of free space truncation in the computation and presented the results of TO in Fig. 15a-f. It can be observed that the variations in the topologies converge as we increase the size of free space truncation. This convergence is best visualized in Fig. 15h, where we presented the outline of the structures. The convergence of the design is also reflected in the evolution of the objective plot in Fig. 16. The differences in the design are due to the variation of the potential distribution inside the bulk when different sizes of free space are considered. To further examine the variation of the potential distribution inside the bulk, we plotted the potential along the diagonal of the design domain in Fig. 17 for different truncation of free space. An important observation is the convergence in the variation as we increase the size of truncation of free space. In our case, we observed the convergence of potential distribution for free space sizes of 4 times ( 240 m) and 6 times (360 µm) of bulk dimension. Finally, the difference in the objective without considering free space and with a free space size of 6 times (360 µm) that of bulk dimension is 24.77%.

Conclusion
In our study, we investigated the impact of the surrounding free space on the structural optimization of EAPs (EAPs). We modeled the free space as the void region in densitybased TO. We increased the size of the free space considered for computation and observed the results of TO. Although there may not be significant visual differences when considering the free space, there are differences in the performance. Considering that these materials will be used as actuators, for example in biomedical devices, these differences might play an important role. Thus, we propose that the consideration of free space has an impact on the optimized structure. To accurately capture the influence of free space, we recommend considering a free space that is at least 6 times the size of the bulk material under consideration. Furthermore, we investigated to assess whether the influence of free space remains consistent as the domain size increases. Our observations revealed that the influence of free space is indeed independent of the domain size.
The consideration of free space also plays an important role in magneto-active polymers. Furthermore, the influence of free space on shape optimization and stress constraints has not been studied yet, and this could be explored using a sequential optimization approach Stankiewicz et al. 2022). This constitutes the direction for future works.