Journal Search Engine
Search Advanced Search Adode Reader(link)
Download PDF Export Citaion korean bibliography PMC previewer
ISSN : 1229-3059(Print)
ISSN : 2287-2302(Online)
Journal of the Computational Structural Engineering Institute of Korea
Vol.32 No.4 pp.223-231


Shape Design Optimization of Electrode for Maximal Dielectrophoresis Forces

Hong-Yeon Jeong1, Seonho Cho1†
1Dept. of Naval Architecture and Ocean Engineering, Seoul National Univ., Seoul, 08826, Korea
Corresponding author: Tel: +82-2-880-7322; E-mail:
May 30, 2019 June 5, 2019 June 10, 2019


A continuum-based design sensitivity analysis(DSA) method is developed for electrostatic problems. To consider high order objective functions, we use 9-node finite element basis functions for analysis and DSA methods. As the design variables are parameterized with B-spline functions, smooth boundary variations are naturally obtained. To solve mesh entanglement problems during the optimization process, a mesh regularization scheme is employed. By minimizing the Dirichlet energy functional, mesh uniformity can be automatically achieved. In numerical examples for maximizing dielectrophoresis forces, the numerical results are compared with well-known electrode geometries and the obtained characteristics are discussed.

최대 유전영동력을 위한 전극의 형상 최적설계

정 홍 연1, 조 선 호1†
1서울대학교 조선해양공학과


정전기 문제에 대한 연속체 기반 설계 민감도 해석(DSA) 방법을 해석적으로 유도하였다. 고차 항을 포함한 목적 함수를 고려하기 위해 해석 및 DSA 방법을 위해 9 노드 유한요소법 기반 함수를 형상 함수로 사용하였다. 최적화 과정에서의 설계 변수를 B- 스플라인 함수로 매개 변수화하여 비현실적인 형상이 아닌 부드러운 경계를 가진 최적 형상을 얻을 수 있었다. 유한요소법을 이용한 최적화 과정에서 일반적으로 발생하는 메쉬 얽힘 문제를 해결하기 위해 메쉬 균일화 기법을 사용하였 다. 이 기법은 디리쉴릿 에너지 범함수를 최소화함으로써 메쉬 균일성을 자동으로 얻을 수 있게 한다. 몇 가지 수치 예제들 을 통해 DEP 힘을 최대화하기 위한 평행판의 최적 형상을 얻어낸다. 이를 기존에 실험적으로 검증된 평행판의 최적 형상과 비교하여 그 특성을 논의하였다.

    Ministry of SMEs and Startups

    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

    2 V = q

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

    V = V 0 o n Γ D

    The trial solution space Y is defined as

    Y = { V H 1 ( Ω ) : V = V 0 o n Γ D }

    and the space Y for the virtual electric potential as

    Y ¯ = { V ¯ H 1 ( Ω ) : V = 0 o n Γ D }

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

    Ω V ¯ · V d Ω = Ω q V ¯ d Ω Γ N ( V ¯ E · n ) d Γ N for a l l V ¯ Y ¯

    A bilinear electric energy form is defined as

    a ( V , V ¯ ) = Ω V · V ¯ d Ω

    and a linear load form as

    l ( V ¯ ) = Ω q V ¯ d Ω Γ N ( V ¯ E · n ) d Γ N

    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 : x x τ , x Ω is given by

    x τ = T ( x , τ )


    Ω τ T ( Ω , τ )

    A design velocity field S that is equivalent to a mapping rate can be defined as

    S ( x τ , τ ) d x τ d τ = d T ( x , τ ) d τ = T ( x , τ ) d τ

    The point-wise material derivative of response V at x Ω is expressed as

    V ˙ d d τ V τ ( x + τ S ( x ) ) | τ = 0 = V + V · S

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

    Consider the general performance functional in domain and boundary integral forms

    ψ 1 = Ω f ( x ) d Ω


    ψ 2 = Γ g ( x ) d Γ

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

    ψ 1 = d d τ Ω τ f τ ( x τ ) d Ω τ | τ = 0 = Ω d d τ { f τ ( x τ ) | J | } | τ = 0 d Ω = Ω { f ( x ) + f ( x ) · S ( x ) + f ( x ) · ( S ( x ) ) } d Ω = Ω { f ( x ) + · ( f ( x ) S ( x ) ) } d Ω


    ψ 2 = d d τ Γ τ g τ ( x τ ) d Γ τ | τ = 0 = Γ d d τ { g τ ( x τ ) | J | J T n } | τ = 0 d Γ = Γ { g ( x ) + ( g ( x ) · n + κ g ( x ) ) S ( x ) · n } d Γ

    where κ = · 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

    [ a ( V , V ¯ ) ] = a ( V , V ¯ ) + a ( V , V ¯ ) + a ( V , V ¯ )

    [ l ( V ¯ ) ] = l ( V ¯ ) + l ( V ¯ )

    where a ( V , V ¯ ) and l ( V ¯ ) 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 = V ˙ V · S , Eqs. (16) and (17) are rewritten as

    a ( V ˙ , V ¯ ) = l S ( V ¯ ) a S ( V , V ¯ ) , V ¯ Y ¯


    a ( V ˙ , V ¯ ) = Ω V ˙ · V ¯ d Ω

    l S ( V ¯ ) = l ( V ¯ ) l ( V ¯ · S )

    a ( V ˙ , V ¯ ) = a ( V , V ¯ ) a ( V · S , V ) a ( V , V ¯ · S )

    The load terms are assumed to be independent of shape variations. By solving equation, the response sensitivity 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

    ψ = Ω g ( V , V , ( V ) ) d Ω

    where VH2(Ω) 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

    ψ = Ω [ g , V V + g , V · V + g , ( V ) · · ( V ) + · ( g S ) ] d Ω

    Using the relation in Eq. (11), the partial derivatives in Eq. (23) can be rewritten in terms of V ˙ as:

    ψ = Ω [ g , V V ˙ + g , V · V ˙ + g , ( V ) · · ( V ˙ ) g , V V · S g , V · ( V · S ) g , ( V ) · · ( ( V · S ) ) ] d Ω + · ( g S ) ] d Ω

    Note that V ˙ , V ˙ , and ( V ˙ ) depend on velocity field S .

    An adjoint equation is introduced by replacing V ˙ Y in Eq. (24) by the virtual electric voltage λ Y and by equating the sum of terms involving λ with the energy bilinear form, yielding the adjoint equation for the adjoint variable λ as

    a ( λ , λ ¯ ) = Ω [ g , V λ ¯ + g , V · λ ¯ + g , ( V ) · · ( λ ¯ ) ] d Ω , λ ¯ Y ¯

    To take advantage of the adjoint equation, evaluate Eq. (25) at λ ¯ V ˙ Y ¯ to obtain the following,

    a ( λ , V ˙ ) = Ω [ g , V V ˙ + g , V · V ˙ + g , ( V ) · · ( V ˙ ) ] d Ω

    Similarly, the design sensitivity equation (18) may be evaluated at V ¯ λ , since both belong to Y , to obtain

    a ( V ˙ , λ ) = l S ( λ ) a S ( V , λ )

    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

    Ω [ g , V V ˙ + g , V · V ˙ + g , ( V ) · · ( V ˙ ) ] d Ω = l S ( λ ) a S ( V , λ )

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

    ψ = S ( λ ) a S ( V , λ ) Ω [ g , V V · S + g , V · ( V · S ) + g , ( V ) · · ( ( V · S ) ) + ( g S ) d Ω

    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.

    F D = 1 2 X ( u 1 2 + u 2 2 ) d X = 1 2 X ( u 1 , x 1 2 + u 1 , x 2 2 + u 2 , x 1 2 + u 2 , x 2 2 ) d X

    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 ( x 1 , x 2 ) X ˜ ( u 1 , u 2 ) 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.

    Γ D + and Γ D are each considered anode and cathode. In anode boundary Γ D + the prescribed electric potential value is 10, and in cathode boundary Γ 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.

    S ( x ) = N I δ x I

    where δ 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 E · 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

    ψ 1 = 1 2 Ω c { ( ( V ) ) · ( V ) } · { ( ( V ) ) · ( V ) } d Ω c

    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

    M a x i m i z e ψ 1 when Ω d Ω Ω c o n s t r a i n t d Ω .

    The shape sensitivity of the volume can be written as

    d d τ Ω τ d Ω τ | τ = 0 = Ω · S d Ω

    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.


    This work was supported by the Technology Development Program(S2591643) funded by the Ministry of SMEs and Startups(MSS, Korea).



    Electrostatic problem


    Variation of domain


    Parallel electrode geometry


    Four electrodes


    Shape optimization result for 4 electrodes


    Shape optimization for initially curved example


    Comparison of design sensitivity with finite difference method


    1. Braibant, V. , Fluery, C. (1984) Shape Optimal Design using B-Splines, Comput. Methods Appl. Mech. & Eng., 44(3), pp.247~267.
    2. Choi, M.-J. , Cho, S. (2015) A Mesh Regularization Scheme to Update Internal Control Points for Isogeometric Shape Design Optimization, Comput. Methods Appl. Mech. & Eng., 285(1), pp.694~713.
    3. Haug, E.J. , Choi, K.K. , Komkov, V. (1986) Design Sensitivity Analysis of Structural Systems, Academic Press Inc., New York.
    4. Lee, S.-W. , Cho, S. (2013) Isogeometric Shape Design Sensitivity Analysis of Mindlin Plates, J.Comput. Struct. Eng. Inst. Korea, 26(4), pp.255~262.
    5. Park, J. , Kim, B. , Choi, S.K. , Hong, S. , Lee, S.H. , Lee, K.I. (2005) An Efficient Cell Separation System using 3D-Asymmetric Microelectrodes, Lab on a Chip, 5(11) pp.1264~1270.
    6. Pohl, H.A. (1978) Dielectrophoresis: The Behavior of Neutral Matter in Non-uniform Electric Fields, Cambridge University Press, Cambridge, UK.
    7. Yoon, G.H. , Park, J. (2010) Topological Design of Electrode Shapes for Dielectrophoresis based Devices, J. Electrost., 68(6) pp.475~486.