Extended Finite Element for Electromechanical Coupling

Extended Finilte Element for Electromechanical

Coupling

Veronique Rochus andDaniel Rixen

Delft

University

of

Technology faculty

3mE

Dpt. of Precision and MicrosystemsEngineering,Engineering Dynamics Mekelweg 2, 2628 CD Delft, The Netherlands

V.Rochus

@tudelft.nl,

D.J.Rixen@tudelft.nl Phone.

+31-(0)

15-278 18 81

Abstract

In MEMS modelling the electro-mechanical coupling takesanimportantplace. Indeedmanydevicesuse electro-static forces as actuator. Thenumerical modelling of this

type ofproblem needs astrong coupling between the

me-chanics and the electrostatic field. In fact when the struc-ture is moving, the electrostatic field around it has to be modifiedin consequence. The first solution isto usefinite element methodtomodel the electrostatic field. Inthiscase

the mesh hastobeupdated depending onthedisplacement of thestructure. Many researches have beenperformedto

deformproperly the electrostatic mesh, but when large dis-placementaretaken intoaccount,theelements become dis-torted. Furthermore, when thepull-in is achieved, the elec-trodes arein contactand thelayer of electrostatic elements istotally squeezed. The second usual solution isto usethe boundary element method tomodel the electrostatic field.

Inthiscase,thereare no moreremeshingproblem, but the

computational time islarger andsingularity problems ap-pearswhen theelectrodes becomein contact.

One solution for this remeshing problem is to use

ex-tended finite elements (X-FEM) which are a new type of elements tailored to simulate problems involving discon-tinuities and moving boundaries. Initially this

methodol-ogy wascreated for crackpropagationproblems [1, 3, 5],

but itsapplicationhas been extendedtoseveral other prob-lems such as elastic problem involving inclusions,

flow-structure interaction and solidification problems. In this

papertheconceptof extended finiteelements is appliedto

develop modelling approaches for the electro-mechanical

coupling. The method will be illustrated here for a

one-dimensionalproblemandimplementationissuesrelativeto

thetwo-dimensionalcase arediscussed.

1. Introduction

This research aims at modelling of electro-mechanical

coupling thatnormallytakesplace in some micro-electro-mechanical systems (MEMS) like micro-resonator or RF

switches. Theproblem canbe described as a conducting mechanicalstructurewithapplied voltage, whichgenerates

asurrounding electrostatic field. The electrostaticfield, in

itsturn, causes an appearanceofelectrostaticforce,applied

to the structure. This type ofproblemis a strongly non-linear problem since theelectric domainchanges with the deformation of the structure. The usual numerical tech-niques to model this type ofelectro-mechanical problem

are the finite element method and the boundary element

method. The mechanical structure is usually simulated

by a finiteelement model and the electrostatic domain is solved by either finite element method or boundary

ele-mentmethod. For both cases someproblemsappearwhen

the structureundergoes large displacements and when the

electrodescomeintocontact. Indeed, the electrostatic finite element mesh has to be modified as the structure moves. Moreover, theelectrostatic meshcanbeseverely deformed

ifthe structure undergoes large displacements.

Further-more, when theelectrodes comeintocontact,theelements between theelectrodes havetobe deleted. Theboundary element method proposes a partial solution to this prob-lem: the electrostatic domain is meshedonlyonthe

bound-ary hence allowing large displacements to the structure.

Howeverthis increases the computation time. Moreover,

when the structure comes into contact, theboundary

ele-ment can no longer be applied since it requires that agap

exists between theelectrodes. Inorderto simplifyand im-prove modelling ofstructuresmoving in an electric field,

we propose to make use of the concept of eXtented

Fi-niteElements. They are a new type ofelements tailored

tosimulateproblems involving discontinuities and moving boundaries.

The basic idea istohaveanelectrostatic mesh covering the entire domain and that doesnotchange while the struc-turepartis moved within thefield. The electro-mechanical problem is consideredas abi-materialproblem where the mechanics and theelectricityarecomputed andcoupledon a single element. Following the variational approach de-velopedin [4], electrostatic forcesmaybe derived and

ap-pliedatthe interface of theelement. The electromechani-calproblemmaythen besolved and the resultscorrespond

verywellto theanalytical solution for a onedimensional problem. A short discussionwill be also given about

im-plementation issues of this techniquesin 2dimensions.

2. eXtended Finite Element

Theory

The extended finiteelement method consistsin discretis-ing the entire electro-mechanical problem with a fixed mesh andinfollowing the interface between twodomains through this mesh. At the interface thephysicalfieldorits gradientare no morecontinuous. Tomodel this disconti-nuity, special shape functions areused toenrich theusual discretisation. For instance, the mechanical displacement

u is enhanced by discontinuous shape functions Mi such

as:

u(x, t)

=

ZNi

(x)

Ui

(t)

+

ZMj

(x,

t)Aj

(t)

whereNi are the standard shape functions and

Mj

are the enriched shape functions taking the discontinuity into

ac-count. Newunknowns

Aj

areintroducedtomodel the

dis-continuity.

There are different ways to createtheseadditional shape functions. Moes [2] proposes to define these functions basedonthe standardshape functionsby the relation:

My(x,t)

=Nj(x)O(x,t)

enrichment0 becomes:

)a

= forO<1 < F

®b=

1

-r

for1<l

<1

Theplotof these functions is illustrated in Figure 2 and (2)

where0(x, t) is defined by:

("-(x,it)

= lVj NN- ,V

i i (3) _a b

where ](x, t) is the levelset fielddescribing the location of the interface, andvi is the value of the level setof the

node i. The level set is described by a surface function

intersecting the problem planatthe interface.

First this methodology will be applied in a

one-dimensional electro-mechanical example and the

mechanical problem will be solved by computing electro-static forces atthe interface. Then this methodology will

beadaptedtothetwodimensionproblem.

3.

Implementation

in

One

Dimension

3.1 Shape Functions in ID

Inthis section theshape functions forapuremechanical

bi-material elementarepresentedinonedimension. Inthat

case the usual shape functions Ni on a reference element

arethefollowing:

0 F

Figure2. Extended finiteelement.

the enrichedshape functionsare:

Mla M2a Mlb M2b (1-1 1r 1-11

fl 1-ll

(6)

The sum of these two shape functions corresponds to the

expression of0(see equation (5)). f Mla+M2a= r

Mlb +M2b= 1 -

(7)

N, = 1-n

N2r= (4)

These shape functions are represented in Figure 1. They

allowstomodel linear behaviourintheelement.

I1I

1

whichcorrespond of the plotinFigure2

3.2 ApplicationtoElectro-mechanicalCoupling

The extended finite elements methodology is now

ap-plied to electro-mechanical coupling in one dimension.

The element is divided intwosub-domains: the mechanical domain called "a"(in grey) whichrepresentsaconducting material and the electrostatic domain "b" (in white)

repre-sentinganon-conductingmedium. The part "b" represents, forinstance, the airinwhich the structuremoves(see Fig-ure3).

0

k

=0

Figure 1: Linear shape functionsin ID. a

U,

b

L U,

Then the enriched shape functions are added to model

the discontinuity. The enriched shape functions,

follow-ing the approach of Moes [2], are given by (2) where 0

is defined by the level setTequal tothe signed distance

between thepointxand thediscontinuitywhich is located

atadistance F from the first node. In the IDmodel, this

Figure 3: Extended finite element.

Inelectro-mechanicalproblems, thephysical unknowns

arethedisplacementu and thevoltage 4. Both fields will

(5)

No I

C =V

be discretised on domaina andb. Hence we assume that there isamechanical field andanelectrostatic fieldinboth

partsand define:

{

Ua Ub

Oa

Ob

N1

U1

+N2 U2+

MlaA

l+M2aA2 N1 U1+N2 U2+MlbA1 +M2bA2

Nl1'l

+

N2P2

+

MlaBi

+

M2aB2

Nl4'1

+N2P42 +MlbB1

+

M2bB2

whereNiarethe standardshapefunctions, and Mia and Mib arethe enrichedshape functions for each domains:

{

N1N2 x L x L M Ma= (1 L) x M2a x Mib- (I- )- (9) M2b LL L/

where1is theposition of the interface.

The unknowns of this problem are U1, U2, (cl and(D2,

thedisplacement and the electric potential at the

extremi-ties of the element. and thenewunknownsA1, A2, B1 and

B2 usedtomodel the discontinuity of the mechanical field and the electricpotential.

3.2.1 Electrostatic Potential

Along the extended element including the conductor

structure and the airgapbetween theelectrodes, the

elec-tric field is discontinuous. Indeed thevoltage iscontanton the conductor(mechanicalparta) and decreases linearlyon the electric domain b.

Apotential difference is applied between the extremities of the element

(Dl

=V and(P2=0)as showninFigure3. The stiffness matrix associatedtothe electrostatic fieldmay

becomputed byintegrationoverbothpartsof the element:

6qT

OO6q5

a60

a6d

I

FLa6O a60

~ ~Ca-jx+~ b dx

¢¢

2

Joax aax

2

Jlax bax

(10) whereCa andEb arethepermittivity of domainaandb, re-spectively. Considering the structure as aperfect

conduc-tor, thevoltageonthisparthastobeconstant. Tokeepthe

voltage constantonthe mechanicalpart,thepermittivityof this domain isimposedtobeaverylargenumbercompared

tothe voidpermittivity. Inthiscase wewill take Ca=1 and

-b=

-0-Solving thepureelectrostaticproblem, the obtained

po-tentialalong the elementisconstantonthe mechanical part and decreaseslinearly between the electrodesasplottedin Figure 4.

3.2.2 Electrostatic Forcesatthe Interface

Tocompute the electrostatic forces appliedatthe inter-face between themechanicalstructureand the electrostatic

domain, an energetic approach has been chosen as

pro-posed bythe author inpaper [4]. This method consists in

determining the electrostatic forcesatthe nodes ofafinite

elementbyanintegrationonits volume. The finite element

formulation is:

fe';ecau DTF (grad6u)dQ

(I

1)

x10

ElementLength [m]

Figure 4: Electrostatic potential inone ex-tendedfiniteelement.

with

Q 0 2aa

F= 8 a ) (12)

where D is the electrostatic displacement. Inone

dimen-sion thisexpression is reducedto:

felec1

2

aJo

(ax

ax

(13)

The samemethod isnow appliedto the extended finite

element. The electrostatic forces iscomputedoneach sub-domain of the elementby:

T 6U I Oa 2 a6adx IL a(b 2 aU

felec6)u=- J(ax)Ca a dx-2 j ( a ) Eb a dx

(14)

The electrostatic potential being constant on the

conduc-tor structure, the firsttermdisappears and theelectrostatic

forcesarecomputed only by the integrationondomain b. 3.2.3 Electro-mechanical Coupling

Now the complete electro-mechanical problem will be

considered. The mechanical stiffness ofan extended

one-dimensional element isobtainedby:

TKuu6

a6u a36u

~L 36u

a36u

6uTKax Ea axdx+ -' b d

2 ax ax 2I ax ax

(15)

whereEa andEb are the Young'smodulus of the domain

a and b, respectively. In the present case, a mechanical

behaviour existsonlyonthe domainaandEbwill beset to

zero.

Theequilibrium position of the electro-mechanical prob-lemmaybe obtainedsolving thesystem:

KUUUI=1 felec

l

Ko54,0

qelec (16)

where thearrayucontains the mechanicaldegrees of free-dom ui andAi and contains (i and Bi. qeiec are the

,

-tzz

chargesonthe nodes. Thissystemofequation is non-linear since theelectrostatic force

felec

isanon-linear function of the electric potential (see equation (14)) and because the position of the interface1 after deformation hastobe taken

into account inthecomputation of theelectrostatic stiffness

Koo

andforce

feiec

b

(D1,D2)

20F

(CI

C)2

15 10 0.5 1 1.5 2 2.5 3 x10

Displacementof theinterface [in]

Figure 5: Displacement of the interface(dots)and

analytical solution(plain line).

InFigure 5 thedisplacement obtained with theextended

finite element method (represented with dots) is compared

totheanalytical solution (plain line). Thesetworesults fit

very well for the stablepartofthecurve.

4.

2D

Electrostatic

Problem

We willnowinvestigate possible shapefunctions for ex-tending the triangular linear elements in twodimensions.

Triangular elements areoften usedinpracticeinautomatic

meshing tools andarethusverycommoninpractical

mod-els.

4.1 Moes' Shape Functions

For a triangular form Moes [2] proposes for 0 the

fol-lowing function:

= ,= iVj Nj(x,y)+ VjiNi(x,y) (17)

i i

wherevi isthe value of the level setof the node i. Letus assumethat,inthe reference coordinatespace,the interface

passesthroughthe nodes (C1,C2)and(D1, D2)where C1 =

0andD1 = 1-D2asshown infigure6. The levelsetmay

be chosenas:

uW4t)=a4 +

bq

+c (8

wherea=(D2 -C2), b =-D1 andc=C2D1 sothatv=0

corresponds to the equation of the interface between the

twodomains. Note that this choice isnotuniquefora,b,c.

Figure6. Extended finite element in 2D.

Figure7: Enrichment hat function 0.

The enrichedshape functions and the linear shape

func-tions are: Mal

(Ni

= I-

4-TL

Ma2,

N1 =

~

Ma3 N3 =qMbl ~~~Mb2 kMb3

-2(b

+

c),(I

-

4-,1)

-2(b+c)rlt -2(b +

c)flq

2(c+a

-

cr)(I- -,)

2(c

+

a0

-ccq)r

2

(c

+

a0,-

cq),q

(19) The hat function 0 modellingthe discontinuity isplotted

infigure 7. "a= ZMai i ®b = ZMbi i (20) We can observe that the line of discontinuity is not

hor-izontal, which seems to indicate that these extended

for-mulationmightnotbe suitabletoimposeaconstant poten-tial (as needed foraconductor inside forinstance). Letus

check if this extended shapefunctions aresuitable for the

Soletusinvestigate ifthe extended fieldabove can be

usedtomodelanelectric potentialconstant on one domain

and varyingontheother. Thepotentialis discretised by:

{ (Pa=Nl(A1+ N2 2+N3(l3+ MlaB1 +M2aB2 + M3aB3

lPb

=Nl4l+

N2'D2

+

N3

(l3+

MlbBl

+

M2bB2

+

M3bB3

(21)

Considering that thepotential isconstant onthe

quadrangu-lardomain a (see Figure 7), we impose

aa;

=0andaOa =0

for all values of4 andrj on the domain a. After

develop-ment weobtain thefollowingconstraints:

Therefore we will introduce two successive changes of

variables correspondingto two isoparametric transforma-tion of thephysicalspace: a first between(X, Y) and ( ,rl), andasecond between(t,rl) and (s, t), aspresentedin Fig-ure 8. Notethat the second transformation is defined

sep-arately for the triangular and the quadrangular part ifthe partitionedelement. (X Y,) ~y tt

01

) (l1 =(D2 and

B1

=B2=B3 <1)

(cl

-

c3l)

(22) 2(b+c)

Onlytwounknowns areleft, asexpected, thepotential(l3 and(l1 = D2.

We will now consider that thepotential is constant on

thetriangular domain b sothat

a0b

=°and

aOb

=0 for all

values of4and rj onthe domainb. Theconditions to have a constantpotentialonthis domainare

B2 =B1 =B3 and B=

cl)1

La- 2 B1

2a

2c

(23)

(X, Y) (-1,1) r ) (XY) (,-3 (0,1) -I11)

In this case, settingthetriangular domain b to a constant

potential (l3 implies arelation between cl1 and(2. This

canbe understoodasfollows.

The conditionB2=B1 =B3 enforces thelinearityof the

discretisationfieldondomainb andas aconsequence,also thelinearity of the fieldonthequadrangular domaina. The relation between cl1 and cl2 implied bythe second set of

constraints in (23) imposes that the values of the electric potentialatthe nodes"1", "2","C" and"D" arecoplanar. So the situation where cl3 is set at a potential V (higher electrode) andcl1 and

cl2

are onthegrounded lower elec-trode, ispossible onlyifthe interface isparalleltothelower electrode.

Thus it is impossibleto impose aconstantpotential

in-side thetriangularpart of the extendedelement anda

lin-earfieldonthequadrangularpart. We musttherefore

con-clude that the shape functions derived from the approach described in [2] not suitable for electro-mechanical

mod-ellinginthevicinity of conductors. Hence inthefollowing section, webuild adifferent extended field tocircumvent thisshortcoming.

4.2 Quadratic Enriched Shape Functions

The underlying idea of this new approach can be best understoodby observing thatifthe extendedelementwas

buildoutoftwofiniteelements(one for the conductor and

onefor the non-conducting part), it would straightforward

to simulate the behaviour of electro-mechanical problem. So,wewilltry to usequadrangularshape functions for the

trapezoidal part andtriangularshape functions for the

tri-angularpart.

(-1,-1) (1,-1) (0,0) (1,0)

Figure 8: Successive transformations for an extended

tri-angularelement.

4.2.1 First change of variables

The first change of variables is identical for both do-mainsaandb andmaybe expressed by therelation:

{ X=N1X1+N2X2 +N3X3

Y=N1Y1+N2Y2+N3Y3

N1 = (1 N2=

N3=1

r'-,)

(24) Thecoordinates of the pointCandDisobtain by

comput-ing the intersection between the levelsetboundary and the edge of the triangle. Inthe second space, theposition of thesepointsare:

{C

'C 0 (XC-X1i) (X3 -X1) T1D (X3-XD) (X3-X2) (XD-X2) (X3-X2) (25) 4.2.2 Second Transformation-Quadrangular Part

The relation between the reference space (s, t) andthe in-termediate space (4,r') is given by the following

isopara-metric transformation: t(s,t)

l

(s,

t) 1Ml

+42M2+tDM3+tCM4

rlMl +h2M2

+TIDM3 +TICM4

The last two shapefunctions will be used to enhanc

solution fieldsince the shape functionsM1 and M2 arn sociatedtothenodes 1 and 2 that already define the an

tude of the shape functionsN1 and N2 in the basic elen

Thetotaldiscretisation of the extended field is thus:

Oa

= N1

(4(S t),

rl

(s,~t))

'D1

+N2

(4(S, t),

rl

(s, t) ) (D2

+N3

(4(S,t),rl(S,t))

(3 +M4

(s,t)BC

+M3

(s,t

Itis easy to verify that the part of the hat function

domaina canbe found by settingBC=BD:

(Da

=

M3 +M4

This is depictedinfigure 9.

(27)

z the

e

as-

npli-nent.

4.2.3 SecondTransformation-Triangular Part

The shape functions for the triangular domain are the

same astheinitialchange ofvariables: G1 =

(I

-s-t) SG2=S

t G3=t

and the

isoparametric

transformationbecomes: 4(s,t) l

(s,

t)

(33)

(34)

3G1 +4CG2+4DG3

3G1

+TicG2

+

IDG3

As for the quadrangular part only the shape function not

~)BD

associatedtothe basic node

(here

G2

and

G3)

will define

(28) the enrichment field. The total enhanced

shape

functions

e in are

thus,

fordomain b

Ob

= Nl

(4(S t),~' (s,~t))

'D1

+N2

(4(S, t),q(s, t))

(D2

+N3

(4(S, t),

rl

(s, t) )

(D3

+G2

(s, t)BC

+G3

(s, t)BD

(29)

(35)

Obviously since the addeddegrees of freedom

BC

and BD

used here and used intheshape functions (28) in domain a areidentical, the continuity of thepotentialfield is guar-antiedwhile the electric field, gradientof thepotential,can

be discontinuous due to the enrichment field. Again we

have that the hatfunction0in the extended theory is

(pb

=

G2

+

G3

(36)

0 0 08

~-~o:2~O-~o0.6

Figure

9: Enrichment of domaina.

Let us now compute the electric stiffness matrix. The

electric stiffness matrix is computed by the relation:

K¢))a= Js BaT(X,y)

aBa(X,Y)dXdY

(30)

where Sa is the surface ofintegrationin the space (X,Y)

and where Ba(X,Y) is the matrix of electric field shape

functions obtained from the derivatives of the potential shape functions. The change of variablesbetween the space

(X, Y)

and the referencespace (s, t)may be obtained by the

relation (24) and (27) which allowus towriteX(s, t), Y(s, t)

and the JacobianJof this transformation. Theelectric stiff-ness matrix maythen be computedonthe reference space

by:

Kta

J.

a

C)

J

Ba(s,t)

det(J)

ds dt

a

(31)

where aNj

Ba(s,

t)

=

as

at DM3 at3 1 at (32)

Theelectric stiffnessmatrix iscomputed by the relation

Ko,b

B

T(X,

y)

bBb(XY)

dXdY (37) whereSb is the surface ofintegration inthe space (X, Y). Using (24) and (34) the electric stiffness matrix may be

computedonthe referencespaceby

K¢¢b = X BbT(,J

TT

(J

1

Bb(s,t)

det(J)dsdt

where aN( Bb(S,t aas at aN3 as aN3 at DG3 as DG3 at (38)

(39)

4.2.4 Simple Verification ofthe Element

Wewill consider theverysimple2D caseof a unit square

domain where avoltage of 1V is imposed on thetop and

where thelower edge is grounded. Half of the domain is

conducting and the domain is modelled with 2triangular extendedelementsasbuiltintheprevious section (See fig-ure 10).

First the interface between conductor and vacuum is

takenparalleltotheelectrodes. The computed electric

po-tential is plotted in figure 11: it is observed that the

V=1V

0.8 06

04-L~~~

_V=O

F'igure

10:

Simple

2Dmodel.

Figure 13: Electrostatic potential computed by

twoeXtended Elements.

1.5

14

05~

0

Figure 11 Electrostatic potential two eXtended Elements.

computed by

Next the electricpotential is also computed for thecase where the interface is not parallel to the extremities as

shown infigures 10. The computed potential is shown in

13. Thepotential is again piecewise linear and clearly the

new shape functions of the eXtended Finite Elements can

properly handle the computation of the electric potential evenifthe interface isoblique.

V=lV

V=OV

Figure 12: Simple 2D model.

5.

Conclusions

The issue of mesh moving is a real challenge when

modelling electro-mechanical devices with finite elements. Thispaperinvestigates the application of the Extended

Fi-nite Elementapproachtomodel the motion ofastructurein

anelectrostatic field. The electro-mechanical forcesare

de-rived from the variationalmethodology proposedin paper

[4].

Whenappliedto asimpleonedimensionalproblem, the eXtended FiniteElement approach finds theexactsolution for the strongly coupled electro-mechanical problem: the

exactelectrostaticpotential along the element is retrieved and, under the action of the electrostatic forcesonthe inter-face, thecorrectrelation between deformation andapplied voltage is found.

For the two dimensional case, we have discussed why the enhancementstrategyproposedin[3] isnotsuitable for modelling the electric field jumpinpracticalproblems. We

propose in this paper different enhancement shape

func-tions thatguarantythat thepotential fieldacross a conduc-tor/vacuum interface canbe properly approximated. Itis expected that thesameapproachcanbeappliedinthree di-mensions. This will be investigate in the future together with the global efficiency of the newelements in solving

theelectromechanicalcouplingofrealmicrosystems.

Acknowledgments

The authorswant toacknowledge thesupportof the

Koi-ter Institute of the Delft University ofTechnologyand of the MicroNed program financed by the ministry of

eco-nomical affairs of The Netherlands. The first author

ac-knowledges the financial support of the Belgian National Fund for Scientific Research.

References

[1] T. BELYTSCHKO AND T. BLACK, Elastic Crack

Growth in Finite Elements with Minimal Remeshing,

International Journal of Numerical Methods in Engi-neering, Vol. 45,No.5(1999),pp.601-620.

[2] N.

MoE~s,

M. CLOIREC, P. CARTRAUD, AND J. F.

REMACLE,Acomputational approachto handle

com-plexmicrostructure geometries,ComputerMethods in

Applied Mechanics and Engineering, Vol. 192 (2003),

[3] N. MOES, J. DOLBOW, AND T. BELYTSCHKO, A

Finite Element Method for Crack Growth without

Remeshing, International Journal ofNumerical

Meth-odsinEngineering,Vol. 46 (1999),pp. 131-150. [4] V. ROCHUS, D. J. RIXEN, AND J.-C. GOLINVAL,

Monolithical Modeling ofElectro-Mechanical

Cou-pling in Micro-Structures, International Journal for Numerical Methods in Engineering, Vol. 65, No. 4

(2006),pp.461-493.

[5]

N. SUKUMAR, T. BELYTSCHKO, C. PARIMI,

N. MOIES, AND U. SHUJI, Modeling Holes and

In-clusions by LevelSets in the ExtendedFinite Element Method,ComputerMethodsinApplied Mechanics and