Analysis of electromechanical systems based on the absolute nodal coordinate formulation

The absolute nodal coordinate formulation (ANCF) approach has been successfully used to analyze bodies undergoing large deformations in multibody dynamics applications. In this study, the ANCF is extended to the analysis of coupled electromechanical systems. To this end, the electrostatic equations are solved by means of conventional plane finite elements, and the ANCF is used to describe the geometrically nonlinear elastic deformation of a thin beam. Bidirectional coupling between electrostatic and elastic domains was introduced using an iterative staggering algorithm. The results illustrate that the ANCF approach can be applied to electromechanical problems when objects are discretized using beam and plate elements. Two numerical examples of microbeams subject to an electrostatic field are used to validate the proposed solution strategy and to reveal characteristic features of fully coupled electromechanical solutions accounting for finite strain theory.


Introduction
Elements of micro-electromechanical systems (MEMS) typically experience a combination of large deformations and large reference motions. Microsensors, actuators, accelerometers, and gyroscopes are examples [15]. Some of these devices make use of the piezoelectric effect as well as pneumatic and thermal actuation. Quite often in MEMS applications, mechanical movement is caused by electrostatic forces. As a result, coupled electromechanical simulations are essential for the successful design and testing of such devices [3,31].
Although the general theory of coupled nonlinear electromechanical analysis is well established (see, for example, [7,26]), there are applications where the theory or its combination with specific numerical algorithms must be updated. One such area is the analysis of thin plate-or beam-like structures. These kinds of structures can often be found in micro-electromechanical systems ( [2,6,16,17,27]).
One possible approach to the analysis of electromechanical systems that undergo large deformations and large reference motions is to apply the absolute nodal coordinate formulation (ANCF), which was originally introduced for multibody dynamics analysis [33]. The formulation is a nonlinear finite element approach that can be used for a wide variety of applications; however, it has been mostly used for beam-and platelike structures [9]. Some of the most recent developments include enhancements in numerical performance [14,18,20,34]. Applications of the ANCF include soft machines [13], biological tissues [29,30], cables and wires [12,19], pneumatic elements [37], dielectric elastomer actuators [24], railroad track structures [33], and rotating shafts [4]. As these applications imply, the ANCF is especially useful when analyzing beam-or plate-like structures that experience large deformations and large rotations.
As shown in the literature, the ANCF has already been successfully used for the coupled electromechanical analysis of some micro-electromechanical systems. ANCF beam elements within the framework of the Euler-Bernoulli beam theory assumption were used to simulate the performance of surface-bonded piezoelectric actuators [10]. Later, [28] successfully extended the ANCF approach to simulate piezoelastic laminated plates. Other attempts to apply the ANCF to MEMS simulations have been reported by [23,38]. However, in all the above studies, application of the ANCF to electromechanical analysis was limited to systems that exploit the piezoelectric effect. Another typical case for MEMS applications, a charged beam-like structure in the presence of an external electric field (which includes calculation of this field outside of the structure under consideration), has not to date been addressed in the ANCF literature.
In practice, the introduction of general electrostatic loads in mechanical simulations performed using the ANCF approach can be accomplished by employing an iterative staggering coupling algorithm [21,35]. The practical importance of the subject is that successful electromechanical coupling with ANCF elements paves the way to high-accuracy simulations of thin structures experiencing large deformations in the presence of an electric field.
The objective of this paper is to introduce a novel computational procedure for solving coupled electromechanical problems involving large deformations in thin beam-like structures. The approach makes use of the absolute nodal coordinate formulation, which has not yet been applied to general electromechanical analysis. Compared to conventional coupled finite element analysis with solid elements, the advantages of the proposed approach are brought about by making the beam assumption and using the ANCF in the mechanical part of the analysis. At the same time, the procedure includes a direct solution of the electrostatic field equations accounting for the geometrical changes that occur as the structure deforms mechanically. This combination of features makes the approach suitable for analyzing various MEMS that include thin structures subjected to large deformation.

Electrostatic formulation
When there is no magnetic field, Maxwell's equations reduce to Faraday's law with no induction term and Gauss law. They can be expressed as (SI units are assumed): where E is the electric field (field strength), D is the electric displacement field or electric induction, ρ is the volumetric charge density, and ∇ is the customary notation for the Nabla operator. Vectors E and D are connected by constitutive relations and in the linear case are proportional with a coefficient called permittivity of the material ε as follows: The solution of Eq. (1.1) is usually sought in terms of the scalar electric potential ϕ, i.e., E = −∇ϕ, which automatically satisfies the first equation in (1.1) and leads to expressing Eq. (1.2) as: For a homogeneous material (i.e., ε = f (x, y)), Eq. (4) can be written in the form of Poisson's equation, using the symbol for the Laplacian operator: A solution of Eq. (5) can be found by applying the finite element method [39]. In a two-dimensional case, for the quadrilateral elements, linear shape functions in local bi-normalized coordinates −1 ξ, η 1 are as follows: They lead to the local matrix of gradients: The Jacobian matrix J, which relates global physical and local bi-normalized coordinates, enables the transition to derivatives with respect to global coordinates: Under the assumption that the considered domain is represented by a homogeneous material with permittivity ε, the permittivity matrix P e for each element can be calculated using numerical integration as follows: After discretization and use of the variational approach, Eq. (5) can be expressed as In Eq. (11), V is the global vector of unknown values of the electric potential at the nodes of the finite element mesh. The correspondence between the positions of unknowns in the global vector V and in the local vector V e (the vector consisting of electric potential at 4 nodes of the element e) can be expressed using the incidence matrix a e : V e = a e V .
Matrix P in Eq. (11) denotes the global matrix of the element permittivity matrices P e (10): where N e is the number of elements in the finite element mesh. Finally, the right-hand side in Eq. (11) is represented by L = L V + L p , which is the global vector responsible for the contribution of point charges L p and volumetric charge density present in the region L p . The latter one is calculated as follows: In Eq. (14), l e can be obtained using shape functions matrix N. It can be expressed as: Note that point charges L p have been introduced in the right-hand side of (11) to provide an easy way to account for the presence of thin structures (beams or shells) in the electrostatic domain when the distribution of charges is known.
Once the solution of (11) is found, it is possible to determine electrostatic forces acting on electrically charged structures. The electrostatic force acting on a point charge q in the presence of on electric field E can be expressed as follows: The electrostriction term is not considered in (16), because the electrostrictive strains are small for most non-ferroelastic materials. In particular, for the silicon, which is the most popular material in MEMS, the authors of [36] report an electrostriction coefficient of 1.8 × 10 −22 m 2 /V 2 . For the strongest electric field considered in our research, 5.7 × 10 6 V/m, that results in electrostrictive strains less than 6 × 10 −9 , it is negligible in comparison with the elastic strains caused by the Coulomb forces.
In terms of the finite element discretization described above, the force acting on a point charge q, placed at node n belonging to element e, can be therefore calculated as: where B| n is the gradient matrix evaluated at the coordinates of node n. For nodes that belong to several adjacent elements, the final value of the electrostatic force can be calculated by averaging the results obtained by applying (17) to each of the elements. Once calculated, the force determined is then used as a nodal force in the structural finite element analysis, which is performed by use of the absolute nodal coordinate formulation (ANCF). Equations (16) and (17) for electrostatic forces require that electrical charges are known, which is not always the case. It is typical in electromechanical devices for the voltage applied to a conducting structure to be controlled without knowing a priori what is the distribution of charges. In that case, the Maxwell stress tensor [11] can be used to find distributed electrostatic forces. The Maxwell stress tensor is a second-order tensor, and in dyadic notation, its part related to the electrostatic field can be expressed as: where E ⊗ E is the tensor product of these vectors, E is the magnitude of the electric field vector E, and I is the second-order identity tensor. The vector of electrostatic traction (negative pressure) acting on a surface of an electrically charged body, characterized by the external normal vector n, can be calculated using the Maxwell stress tensor at this point as: In the case of a conductor with a prescribed value of the electric potential on its surface, which is often the case in electromechanical devices, the electric field E is perpendicular to that surface and co-linear with the normal vector n. In that case, substitution of (18) into (19) results in a simple expression for the magnitude of electrostatic traction (negative pressure) p, Using finite element discretization, for a point with local coordinates (ξ, η) on the surface of an element e, the value of the electrostatic traction in matrix notation can be calculated as follows: The electrostatic pressure defined by (21) for each finite element makes it possible to calculate an equivalent distribution of nodal forces, which are then applied in the structural finite element analysis in the same way as the forces defined by (17).

Absolute nodal coordinate formulation
In the absolute nodal coordinate formulation, the position field, r, and the displacement field, u h , of an arbitrary point within an ANCF beam element in the current (deformed) configuration in the Cartesian coordinate system can be interpolated via the isoparametric property using equal-order interpolations as follows: where N m is a shape function matrix, q is a vector of nodal coordinates, u is a vector of nodal displacements, x = {x, y, z} are the physical coordinates, and ξ = {ξ, η, ζ } are the bi-normalized local coordinates. As implied by its name, the absolute nodal coordinate formulation uses absolute positions and their partial derivatives with respect to local coordinates (i.e., slope vectors) as nodal coordinates. For a node i of the three-node three-dimensional element illustrated in Fig. 1, these coordinates are: where the following notations for slope vectors are used: These slope vectors, known also as extensible director vectors in FE-literature, define the orientation and the deformation of the beam cross section. This is a trademark of the ANCF. Similar approaches have also been studied in finite element literature [5,32]. The kinematics of a beam element based on the ANCF is depicted in Fig. 1. As the Figure shows, the position field can be described using the initial position R and displacement field u h vectors as The shape function matrix for a d-dimensional element based on n base functions can be written as follows: where I is a d × d identity matrix and N 1 . . . N n are element-specific shape functions. Considering the relationship in Eq. (25), the deformation gradient F can be written as where I is an identity matrix and ∇u h is the material displacement gradient. The Green-Lagrange strain tensor can be written as A vector of elemental internal forces can be found by where S is the second Piola-Kirchhoff stress tensor. It can be expressed for a linear elastic isotropic St. Venant-Kirchhoff material as where λ and G are Lamé parameters. In this study, a two-dimensional linear two-node ANCF element with selective reduced integration scheme [25] is used, because it is compatible with the four-node bilinear FE plane element in the subject electromechanical analysis. The shape functions can be written as follows: where parameters l x and l y are the length and height of the element in the local physical coordinates x = {x, y}.
The local natural (bi-normalized) coordinates ξ = {ξ, η} are used with numerical integration over the volume for practical reasons yielding ξ = 2x l x , η = 2y l y . The approach proposed in [8,25] to eliminate Poisson locking is used here. The material coefficients matrix based on the plane stress assumption is divided into material coefficients matrix D 0 where the Poisson effect is not taken into account, and material coefficients matrix D v , where E is Young's modulus and ν is Poisson's ratio. Due to the plane stress assumption, the stress-strain constitutive relation is written via Voigt notations as follows: where ε = [E 11 E 22 2E 12 ] T is a strain vector. A vector of elemental internal forces for the selective reduced integration scheme can be written as where J e = det ∂ R ∂ξ is the Jacobian determinant, and stress tensors S 0 and S v consist of components of the stress vectors σ 0 and σ v . Due to the highly nonlinear nature of the internal forces in Eq. (34), the elemental stiffness matrix is found by applying a forward finite-difference scheme. The global system of mechanical equations in discretized form can be written as: where K is the global stiffness matrix, U is the global vector of nodal displacements, and F is the global vector of external forces, which includes both internal forces calculated according to (34) and nodal electrostatic forces either calculated according to (17) or found from the electrostatic pressure (21), depending on the problem under consideration. On the basis of the above, the electromechanical problem is solved according to the following procedure: 1. Solve the electrostatic FE equations (11) for the unknown nodal potential V . 2. Obtain gradients and calculate electrostatic forces according to (17) or (21). 3. Solve the mechanical FE equations (35) for the nodal displacements U. Note that the right-hand part of (35) incorporates electrostatic forces calculated in step 2. 4. Update the FE mesh used for the electrostatic solution in accordance with the obtained displacements U. Since the displacements U are only calculated for nodes of the structure presented in the mechanical analysis, an appropriate smoothing technique must be used to ensure an acceptable shape for all elements used in the electrostatic analysis after the coordinates update. 5. Return to step 1 of the algorithm with the updated FE mesh.
The loop is broken when the difference in results of two consequent structural analyses is below a certain value (convergence criterion). This method of coupling is usually referred to as staggered algorithm and has been proven to provide reliable results in numerous cases [21,35].

Numerical examples
This Section introduces two numerical examples aimed to test the proposed approach. Both of them represent relatively simple electromechanical problems, which are solved by applying the coupled technique where the ANCF element is loaded with electrostatic forces using a staggering algorithm. The implementation of the algorithm was accomplished via MATLAB.
The first example is the bending of an electrically charged cantilever beam with a permanent charge distribution under the influence of an external electric field. This problem is similar to the problems that are often used to illustrate electromechanical coupling [1,3,22] and allows direct calculation of the electrostatic nodal forces according to (17). Figure 2 depicts the analyzed system. A 200 µm-long beam, with 15µm width and 4µm thickness, Young's modulus 169 GPa and Poisson's ratio 0.3, carries the electrical charge 1.05 × 10 −11 C/µm and is placed in the external electrical field 0.1 V/µm. The idea is to calculate the deformed shape of the microbeam subjected to the action of electrostatic forces by means of the coupled electromechanical analysis described above.
As a first step, a finite element solution of the electrostatic problem is obtained using the approach described in the previous Section. The fact that the dimensions of the cross section of the beam are very small in comparison with the length of the beam makes it difficult to represent the beam's real geometry in the electrostatic domain. Instead, the beam is represented by a linear set of nodes in the two-dimensional electrostatic mesh. The finite element mesh of the electrostatic problem is presented in Fig. 3. The geometry of the beam and corresponding 21 nodes, where concentrated 10 −10 C charges are applied, is highlighted in the Figure.  Solving the finite element equation (11) for the system determines the distribution of electric potential in the domain under consideration. This is presented in Fig. 4, which reveals that the point charges that simulate the electrically charged beam disturb the linear distribution of the electric potential. For comparison, the same electrostatic problem was solved using the commercial ANSYS software, which resulted in an almost identical distribution of electric potential. The use of quadratic elements by ANSYS explains the slight difference in the geometry of isolines close to the point charges (Fig. 5).
Knowing the obtained distribution of the electric potential in the area around the charged beam makes it possible to calculate the electrostatic forces that act at the nodes of the beam according to (17). These forces were then used as input for the structural analysis, which was performed as described in Sect. 3. After the analysis, the displacements obtained at the nodes of the beam were used to update the electrostatic mesh, and both analyses were repeated in a loop, i.e., via a staggering algorithm. The loop was repeated until the difference in displacements obtained for two consequent structural analyses was below 10 −10 µm. In the considered example, convergence was achieved after 19 iterations. The obtained deformed shape of the beam and corresponding mesh of the electrostatic domain from the last iteration of the staggering algorithm are presented in Fig. 6. To avoid distortion of elements due to large deformation of the beam, a simple mesh morphing algorithm was implemented. In this algorithm, all electrostatic mesh node positions were updated with increments in coordinates proportional to displacements of the nodes belonging to the beam. The proportionality coefficients depend on the distance to the nodes. Figure 7 illustrates the solution of the electrostatic problem for the last iteration of the staggering algorithm, which clearly reflects the deformation of the beam.
To verify the solution strategy, a series of problems was solved with different voltages applied at the top boundary of the electrostatic domain. The problem was solved with fully coupled or uncoupled solutions and geometrically linear or nonlinear approaches using the ANCF. "Fully coupled" refers to the coupling approach that makes use of the staggering iterative algorithm. "Uncoupled" refers to the approach in which the electrostatic problem is solved only once for the undeformed configuration. The forces determined are then used for the structural analysis. "Linear" and "nonlinear" refer to the kinematics of the continuum that is used in the ANCF finite element simulations. In the nonlinear approach, the ANCF solution was obtained iteratively in the framework of finite strain theory, using the Green-Lagrange strain tensor in the constitutive relations of otherwise linear elastic material, and updating the stiffness matrix at each iteration. To reveal the importance of this approach, a linear solution was also found based on infinitesimal strain theory. Obtaining the linear solution does not require iteration within each step of the coupling algorithm. The comparison of the four solutions is presented in Fig. 8.
As it is seen from Fig. 8, the uncoupled solutions provide somewhat unrealistic results even for an applied field of moderate strength. The uncoupled approach results in deflections that are almost 2 times lower than those obtained via the fully coupled approach. Figure 8 also reveals the greater importance of using finite strain theory for a fully coupled analysis. The difference between linear and nonlinear solutions, though small for the uncoupled case, becomes considerable when a fully coupled analysis is performed for large values of applied voltage. Among the four presented solutions, only the result of uncoupled linear ANCF analysis is in fact linear. The other solution based on the infinitesimal theory, though linear in the elastic domain, is nonlinear overall because of coupling with the electrostatic domain.
The second numerical example considers voltage being applied to a conducting structure. The applied voltage results in a redistribution of charges, which may change as the structure deforms and the configuration of the electrostatic field changes. Such problems are typically considered in MEMS literature, and the particular problem considered in this paper is about an 80 µm-long silicon cantilever microbeam with 10 µm width and 0.5 µm thickness, which is positioned 0.7 µm above the ground plane [1,22]. The following mechanical properties were assumed for the silicon beam: Young's modulus 169 GPa and Poisson's ratio 0.3. A voltage is applied to the beam, which results in the redistribution of charges, generates an electrical field in the gap between the beam and the grounding plane and produces electrostatic traction that results in the beam bending. The value of applied voltage is increased until convergence of the coupled algorithm cannot be achieved. This means that a solution with a positive gap between the beam and the ground cannot be found, and the beam  The resulting minimum gap between the beam and the ground plane (initial gap-maximum deflection) as a function of applied voltage-comparison with results reported in [22] will in reality contact with the ground plane. That situation is usually referred to as "pull-in effect," and the corresponding value of electric potential is called "pull-in voltage." Results of the simulation in comparison with the results reported in [22] for the same problem are presented in Fig. 9. (Numerical values for the results by [22] were measured from the figure in the published article.) A good correlation is observed between the results obtained with the proposed technique for the coupled analysis and those obtained in [22], which explores the meshless finite cloud and boundary cloud methods (Fig. 9). For the considered problem, the authors of [22] report the pull-in voltage being 2.35 V. In the simulation performed using the approach introduced here, 2.35 V was the largest voltage for which an equilibrium was found, which implies that the pull-in voltage is insignificantly higher (less than 2.36 V). Results obtained for the same problem using the finite element and boundary element methods [1] yielded a pull-in voltage of 2.39 V. The conclusion is that all three solution strategies mentioned provide numerical results that correlate well.

Conclusions
This work introduces a finite element procedure for coupled electromechanical analysis based on the absolute nodal coordinate formulation (ANCF). A staggering algorithm was used to couple two-dimensional electrostatic analysis with mechanical analysis for beam-type structures. The procedure was implemented in MATLAB. Two numerical examples of microbeams subject to an electrostatic field were analyzed using the proposed solution strategy. Comparison with published results obtained for the same problem by other numerical techniques reveals good correlation and validates the proposed approach. Results obtained for various analysis options highlight the importance of using fully coupled solution strategies for such problems even if the applied field is of moderate strength. The difference in the results produced by geometrically linear and nonlinear ANCF beam elements strongly depends on the intensity of the deformation, and for the conditions when electrostatic loads are relatively small, a linear approach may be used.
The implemented procedure and results of its testing on model problems illustrate the applicability of the ANCF approach in coupled electromechanical finite element analysis. Because of the intrinsic properties of the ANCF, the proposed approach is able to offer good efficiency for problems where beam-type structures are subjected to large deformation in an electrical field.