## 1. Introduction

The electric field generated by the electrode can be applied in various engineering fields. One of the examples is the dielectrophoresis(DEP) phenomenon, which makes it easy to handle small particles that have traditionally been difficult to handle with mechanical forces, by forming a suitable electric field according to the purpose(Pohl *et al*., 1978). This means that the particles could move under a non-uniform electric field in a specified direction. Studies dealing with small particles using the DEP have widely been carried out as follows: particle separation, DNA analysis, and so on. To further extend these studies, it is important to obtain an appropriate electric field for the purpose. In general, the influence of electric field is also influenced by the material properties of the particles or medium to be treated, the shape of electrode is also significantly influenced as well(Park *et al*., 2005). Therefore, this paper aims to study the electrode with the optimal shape to produce the appropriate electric field according to the purpose.

To achieve this goal, Yoon(Yoon *et al*., 2010) used a topology optimization technique to obtain the non-uniform electric field of high intensity. However, topology optimization has disadvantages in that the boundary of the obtained optimal shape is discrete, the distinction of the structure boundary is not clear, and the optimal shape with low productivity is obtained. There is also a paper on the optimal shape of electrodes using a genetic algorithm that requires a high performance computation device and long computation time, which shows a limitation of the method. Therefore, we employ a shape design optimization method.

In this study, we study a gradient-based shape design optimization method using finite element method. The easiest choice of design variable in the shape design optimization is the nodal coordinates of the finite element model. In this case, however, it is difficult to obtain an unrealistic model with irregular boundary due to too many design variables and to maintain an adequate finite element mesh in the optimization process. Therefore, to resolve this issue, the design variables are parameterized by b-spline functions. Also, in the shape design optimization, boundary variations are performed at each optimization iteration. The nodal coordinates in the domain must be updated to prevent mesh distortion problem. The re-meshing has the disadvantage that the optimization convergence may be deteriorated due to the sudden change of the mesh structure and the loss of original nodal connectivity. To overcome the aforementioned issues, we employ a mesh regularization scheme(Choi *et al*., 2015). The method minimizes a Dirichlet energy functional to obtain an effective mapping between the parametric and the physical domains. The variance of the Jacobian is minimized to obtain the mesh uniformity.

This paper proposes a continuum-based design sensitivity analysis method developed for electrostatic problem. By applying the techniques that are used to obtain smooth boundaries and lower numerical instability, we verify our proposed method through optimal electrode design and optimization history in numerical examples. Finally, we draw some conclusions derived from this research and present possible topics that can be studied in the future.

## 2. Shape deisgn optimization

We describe the governing equations and finite element formulations for the electrostatic problems to obtain the optimal shape of electrode for the desired performance. A continuum-based shape design sensitivity expressions are derived for the general performance functional using both the direct differentiation method and the adjoint method.

### 2.1 Variational formulation of electrostatic problem

Consider the electrostatic problem shown in Fig. 1.

The general governing equation of electrostatic problem is

In addition to the governing equation (1), the dirichlet boundary conditions are imposed as

The trial solution space *Y* is defined as

and the space *Y* for the virtual electric potential as

Using the virtual electric potential *V* that satisfies the homogeneous boundary conditions, the weak form of Eq. (1) is written as

A bilinear electric energy form is defined as

and a linear load form as

Since the objective function can include second derivative of the response, we use quadratic shape functions for approximating the geometry and the response.

### 2.2. Material derivative

Consider the variation of domain from an original domain *Ω* to a perturbed domain *Ω _{τ}* as shown in Fig. 2.

Suppose that only one parameter *τ* defines a transformation *T*. The mapping $T:\overrightarrow{x}\to \overrightarrow{{x}_{\tau}},\overrightarrow{x}\in \Omega $ is given by

and

A design velocity field $\overrightarrow{S}$ that is equivalent to a mapping rate can be defined as

The point-wise material derivative of response *V* at $\overrightarrow{x}\in \Omega $ is expressed as

where *V*′ and ∇*V* are the partial derivative and gradient of *V* , respectively.

Consider the general performance functional in domain and boundary integral forms

and

The first order variations with respect to the shape design parameter *τ* are derived as

and

where $\kappa =\nabla \hspace{0.17em}\xb7\hspace{0.17em}\overrightarrow{n}$ is used. For a detailed derivation of shape sensitivity expressions, the interested readers may refer to the reference(Haug *et al*., 1986)

### 2.3 Shape sensitivity analysis-direct differentiation method

Taking the first order variation of the bilinear electric energy form and linear load form which respect to shape design parameter *τ*, we have followings

where ${a}^{\prime}\left(V,\overline{V}\right)$ and ${\mathcal{l}}^{\prime}\left(\overline{V}\right)$ denote the explicit variation terms with the dependence of their arguments on the shape design parameter suppressed. For brevity of the problem, the external loads are assumed independent of shape variations. Using the linearity of energy and load forms and the fact that ${V}^{\prime}=\dot{V}-\nabla V\hspace{0.17em}\xb7\hspace{0.17em}\overrightarrow{S}$, Eqs. (16) and (17) are rewritten as

where

The load terms are assumed to be independent of shape variations. By solving equation, the response sensitivity $\dot{V}$ can be obtained. However, in general, not all the response sensitivities are required to obtain the sensitivity of the performance measure. Thus, the direct differentiation method is somewhat inefficient in terms of computation costs. To improve such problem, the following adjoint variable method can be used to efficiently compute the design sensitivity of a specified performance function.

### 2.4 Shape sensitivity analysis-adjoint variable method

Consider a general objective function that may be written in integral form as

where *V*∈*H*^{2}(*Ω*) and function *g* is continuously differentiable with respect to its arguments. By using the material derivative formulas from Eq. (14), the variation of the functional in Eq. (22) can be obtained as

Using the relation in Eq. (11), the partial derivatives in Eq. (23) can be rewritten in terms of $\dot{V}$ as:

Note that $\dot{V},\nabla \dot{V}$, and $\nabla \left(\nabla \dot{V}\right)$ depend on velocity field $\overrightarrow{S}$.

An adjoint equation is introduced by replacing $\dot{V}\in \overrightarrow{Y}$ in Eq. (24) by the virtual electric voltage $\lambda \in \overrightarrow{Y}$ and by equating the sum of terms involving *λ* with the energy bilinear form, yielding the adjoint equation for the adjoint variable *λ* as

To take advantage of the adjoint equation, evaluate Eq. (25) at $\overline{\lambda}\in \dot{V}\in \overline{Y}$ to obtain the following,

Similarly, the design sensitivity equation (18) may be evaluated at $\overline{V}\in \lambda $, since both belong to *Y* , to obtain

Recalling that energy bilinear form *a*(∙, ∙) is symmetric in its arguments, we can conclude that the left sides of Eqs. (26) and (27) are equal. Thus, the right sides of both equations yield the following useful relation

By substituting Eq. (28) into Eq. (24), the expression of *ψ*′ in terms of *V*, *λ*, and $\overrightarrow{S}$ is obtained as

Note that the evaluation of the design sensitivity formula in Eq. (29) requires solving the variational Eq. (5) for *V* . Similarly, the variational adjoint Eq. (25) must be solved for adjoint variable *λ*. Using finite element analysis, solving for *λ* is efficient if the boundary-value problem for *V* has already been solved, since all that is required is adapting the solution to the same set of finite element equations with a different right side(adjoint load).

## 3. Shape design optimization techniques in finite element analysis

We introduce techniques used to obtain smooth boundaries and lower numerical instability in finite element shape optimization. First, we describe the shape-design parameterization where design variables are parameterized using B-spline functions. Second, we describe how the mesh regularization scheme in this method is applied in finite element framework.

### 3.1 Design variable parameterization

In finite element analysis, the nodal coordinates on boundary are usually taken as shape variables. However, this independent nodal movement approach makes very large set of design variables since all nodes on the design boundary must be adopted as design variables. This means that the design variable number on the fine mesh increases exponentially. Also the conventional method has a tendency to generate unrealistic wiggly designs due to large design space. This unrealistic design not only lowers the manufacturability of the resulting optimal shape, but also does not meet geometrical requirements on the regularity of the boundaries.

To resolve such problems in taking nodal coordinates on the boundary as shape variables, we use the parameterization technique that is not to take all nodal coordinates on boundaries as design variables, but to take only some representative nodal coordinates as design variables. Although various functions can be used as parametric functions, this paper adopts the B-spline function as a parametric function(Braibant and Fluery *et al*., 1984). The B-spline function is easy to express complex geometry with a small set of design variables and is flexible enough to satisfy the regularity constraint at the boundary. Additional control points can be introduced without increasing the degree of the curve. Therefore, the parameterization of design variables using B-spline functions avoids unrealistic design. Also, the H-refinement property of the b-spline curves makes it possible to easily calculate the interpolation values between design variables.

### 3.2 Mesh regularization scheme

In shape design optimization, since the structural boundary varies at every iteration, maintaining high mesh quality is an important issue. Thus, it is still an inherent challenge to update inner node coordinates after boundary variations. Generally re-meshing scheme is used to obtain high mesh quality if the boundary variation is significantly large and complicated. However, re-meshing scheme should not be frequently used, since the large variation of mesh structure may lead to sudden variations of the objective function or violation of the constraints which prevents smooth convergence to an optimal shape. Furthermore, when the re-meshing scheme is used, the connectivity between elements and nodes changes. Therefore, when performing shape optimization, various methods of maintaining high mesh quality without re-meshing scheme have been studied. In this paper, we use a mesh regularization scheme that minimizes the following Dirichlet energy functional with convexity properties.

Since the Eq. (37) is convex, the minimization of equation is equivalent to satisfying Laplace’s equation. If it satisfies Laplace’s equations, mapping $F:X\left({x}_{1},{x}_{2}\right)\to \tilde{X}\left({u}_{1},{u}_{2}\right)$ can be achieved and also bijective mapping can be obtained by the RKC theorem.

## 4. Numerical examples

In this chapter, we demonstrate the effectiveness of the methodologies presented in this paper using various numerical examples. Section 4.1 uses the sensitivity formulation derived from chapter 2. The obtained expressions are verified by comparing with the finite difference ones. In section 4.2, shape design optimization is performed on 2D electrodes using the shape sensitivities derived from chapter 2 and the techniques described in chapter 3.

### 4.1 Shape design sensitivity verification

The derived shape sensitivity is compared with the finite difference one to verify the derived analytical expressions. Since the shape sensitivity provides the search direction, it has a significant role in gradient based optimization. The basic shape of electrode is thin parallel plates. Therefore, in this paper, we try to model and interpret the most basic 2D parallel plate structure as the following figure.

${\Gamma}_{D}^{+}$ and ${\Gamma}_{D}^{-}$ are each considered anode and cathode. In anode boundary ${\Gamma}_{D}^{+}$ the prescribed electric potential value is 10, and in cathode boundary ${\Gamma}_{D}^{-}$ the prescribed electric potential value is 0. The parallel electrode geometry model is discretized with 2025 elements and 8281 nodes. To obtain the necessary design velocity field for the shape DSA, the shape variation in Fig. 3(a) is considered. The black arrows represent velocity field which is comprised of all vertical directions to the surface on Dirichlet boundary. As using the design variable parameterization method, we don’t need to take all nodes on Dirichlet boundary as design variables to express the design velocity field. We take only 32 nodes to represent design velocity field. The design velocity field is easily obtained from the perturbed node positions.

where $\delta \overrightarrow{{x}_{I}}$ s the difference of node positions between the original and the perturbed models. So, the design velocity field is also interpolated by the same general 9-node quadratic shape function which is used geometry field expression. The amount of design perturbation is 0.01% and a linear velocity field is assumed. In Table 1, the shape sensitivity of electric voltage from equation is compared with the finite difference results. The objective function is the electric voltage on some selected points shown in Fig. 3(b). As shown in Table 1, a good agreement of shape sensitivity is observed.

Now that we have obtained very accurate and efficient sensitivity as shown Table 1, we are ready to perform shape design optimization for the electrostatic problem.

### 4.2 Shape design optimization for 2D electrodes

Some numerical examples for the shape optimization are demonstrated to show the applicability and effectiveness of the proposed method. In all examples, the permittivity of medium is assumed a unit value. If the only one kind of medium exists, ϵ can be considered a constant value which can be eliminated in governing equation. And also, all kind of load terms in equation are assumed equal to zero. We only consider electrode shape optimization where electric charge density *q* and electric flux $\overrightarrow{E}\hspace{0.17em}\xb7\hspace{0.17em}\overrightarrow{n}$ on surface don’t exist. Consider an electrostatic problem where we have only an anode of 10V and a cathode of 0V. Fig. 4

For first example, the objective is to find an optimal design for a maximum DEP force. The DEP force is one of the useful techniques for dealing with the movement of fine particles. Thus, we make the following objective function since the DEP force has a gradient of the square of electric field as its variable

where *Ω _{c}* is an element located in the center of the whole domain. Thus, the shape optimization problem for maximum DEP force at center point is stated as

The shape sensitivity of the volume can be written as

Consider the following 4-electrode square model.

The black arrows represent design velocity field and only 22 nodes on Dirichlet boundary are taken as design variables for design parameterization method. The initial model is discretized with 3,721 elements and 900 nodes. The B-spline function order for design variable parameterization is 8 and the constraint volume is 50%. Fig. 5 shows optimal shape and other results with initial model.

As shown in Fig. 5, it can be seen that the force of the DEP force tends to increase due to the concave shape of the central part. This is similar to the smooth quadratic curve obtained by assuming a polynomial function for the same example by Y. Huang(Huang *et al*., 1991). The research analytically obtained the shape of the electrode that maximizes the DEP force in the paper. In fact, this shape has been experimentally verified to maximize the DEP force. Therefore, it can be seen that the simulation results in this study agree with the known experimental results.

To proceed with the optimal design in a more complex geometry, consider the case where the initial geometry is curved. In general, it is known that the shape smoothly inward curved is close to the optimal shape for maximizing DEP Force. Therefore, we investigate the influence of the shape of the electrode on the size of the DEP force by performing the optimization in the commonly known curved geometry. The initial geometry with proper curvature, which does not distort the mesh distribution, is made to be Fig. 6(a) and optimization is performed shown in Fig. 6(b). The design velocity field is the same as in the previous example and the order of B-spline function for design variable parameterization is 9 and the constraint volume is 50%.

As can be seen in Fig. 6, an electrode shape that is more inward than a smooth curve which is obtained in previous example. From this result, it can be seen that as the curve has more concave curvature inward, the magnitude of the DEP force increases. Therefore, if a higher magnitude of force is required to use the DEP force, it would be better to form the electric field by designing the center of the electrode more concave.

## 6. Conclusions

Using the finite element method and the efficient adjoint DSA method, a shape sensitivity analysis is derived for electrostatic problems in a steady state. In the FEM-based approach, it is possible to vary the boundary of a structure by directly placing the position of the nodal points for the design variables. In the FEM-based shape optimization, the shape of the boundary is wiggled due to many design variables, so the mesh quality is degraded at every optimization iteration. To solve this problem, design variables were parameterized by B-spline function and minimized Dirichlet energy function for every optimization iteration to maintain high mesh quality. Since the sensitivity of an analytically derived electrostatic problem is applied to an optimization algorithm, we can obtain a smooth convergence history. Through shape optimization research, we have obtained an optimal electrode shape that produces an electric field distribution. The achievement of this research will help to produce electrodes for the purpose of future research.