Modeling Electromagnetics:
Electrostatics

Sebastián García, Software Developer, WRI.

(RowDefault[, (Spacer1), , (Spacer1), ])  #WolframTechConf

Agenda

Workflow for electrostatic models

Discontinuities

Multiphysics actuator

Electrostatics

Electrostatic equation:

∇·(-(Subscript[ϵ, 0])∇V(X)+P(X))=(Subscript[ρ, v])

PDEComponent and boundary conditions

ElectrostaticPDEComponent

ElectricPotentialCondition

ElectricFluxDensityValue

ElectricSymmetryValue

ElectrostaticPDEComponent

({{V,       , electric potential,       , [V]}, {X = {x, y, ...},       , spatial independent variables,       , [m]}, {P,       , "Polarization",       , [C/(m^2)]}, {(Subscript[ϵ, r]),       , "RelativePermittivity",       , [-]}, {(Subscript[P, r]),      , "RemanentPolarization",       , [C/(m^2)]}, {d,       , "Thickness",       , [m]}, {A,       , "CrossSectionalArea",       , [(m^2)]}, {(Subscript[ϵ, 0]),       , "VacuumPermittivity",       , [F/m]}, {(Subscript[ρ, v]),       , "VolumeChargeDensity",       , [C/(m^3)]}})

ElectrostaticPDEComponent[{V[x,y,z],{x,y,z}},<|"RelativePermittivity"->(Subscript[ϵ, r]),"VolumeChargeDensity"->(Subscript[ρ, v])|>]

Electrostatic monograph

Electrostatics monograph:

Overview example (3D capacitor)

Electrostatic equation

Sources

Secondary Quantities

Boundary Conditions

Multiple materials

What happens at the interface when fields are computed?

Discontinuities

To compute the fields D and E:

E=-∇V

D=(Subscript[ϵ, 0])(Subscript[ϵ, r])E=-(Subscript[ϵ, 0])(Subscript[ϵ, r])∇V

What if the relative permittivity is inhomogeneous in the model?

Cross-section of parallel  plate capacitor

Parallel plate capacitor

LC=2;
wdC= 0.5;
d0C = 0.1;
thC=0.4*d0C;
regionMarkersC=<|"dielectric"->1,"air"->2|>;
meshC=ToElementMesh[ToBoundaryMesh[],];

varsC={V[x,y],{x,y}};
parsC=<|"RelativePermittivity"->Piecewise[{{7,ElementMarker==regionMarkersC["dielectric"]}},1]|>;

Cross-section of parallel  plate capacitor

Parallel plate capacitor

LC=2;
wdC= 0.5;
d0C = 0.1;
thC=0.4*d0C;
regionMarkersC=<|"dielectric"->1,"air"->2|>;
meshC=ToElementMesh[ToBoundaryMesh[],];

varsC={V[x,y],{x,y}};
parsC=<|"RelativePermittivity"->Piecewise[{{7,ElementMarker==regionMarkersC["dielectric"]}},1]|>;

bcs={ElectricPotentialCondition[ElementMarker==2,varsC,parsC,<|"ElectricPotential"->0|>],ElectricPotentialCondition[ElementMarker==3,varsC,parsC,<|"ElectricPotential"->1*^4|>]};

Cross-section of parallel  plate capacitor

Parallel plate capacitor

LC=2;
wdC= 0.5;
d0C = 0.1;
thC=0.4*d0C;
regionMarkersC=<|"dielectric"->1,"air"->2|>;
meshC=ToElementMesh[ToBoundaryMesh[],];

varsC={V[x,y],{x,y}};
parsC=<|"RelativePermittivity"->Piecewise[{{7,ElementMarker==regionMarkersC["dielectric"]}},1]|>;

bcs={ElectricPotentialCondition[ElementMarker==2,varsC,parsC,<|"ElectricPotential"->0|>],ElectricPotentialCondition[ElementMarker==3,varsC,parsC,<|"ElectricPotential"->1*^4|>]};

VfunC=NDSolveValue[{ElectrostaticPDEComponent[varsC,parsC]==0,bcs},V,{x,y}∈meshC]

DiscontinuousInterpolatingFunction

VfunCDIF = DiscontinuousInterpolatingFunction[VfunC]

DiscontinuousInterpolatingFunction

VfunCDIF = DiscontinuousInterpolatingFunction[VfunC]

E=-∇V

EfieldC=-Grad[VfunCDIF[x,y],{x,y}];

EvaluateOnElementMesh

D=(Subscript[ϵ, 0])(Subscript[ϵ, r])E

Dfield=Through[EvaluateOnElementMesh[{x,y},8.854*^-12*parsC["RelativePermittivity"]*EfieldC,meshC][x,y]];

Dfield/.{x->0.25,y->0}

EvaluateOnElementMesh

D=(Subscript[ϵ, 0])(Subscript[ϵ, r])E

Dfield=Through[EvaluateOnElementMesh[{x,y},8.854*^-12*parsC["RelativePermittivity"]*EfieldC,meshC][x,y]];

Dfield/.{x->0.25,y->0}

Dfield=Through[EvaluateOnElementMesh[{x,y},8.854*^-12*parsC["RelativePermittivity"]*EfieldC,meshC,"MarkerPriority"->{2,1}][x,y]];
Dfield/.{x->0.25,y->0}

Modeling of electrostatically
actuated MEMS

Sumant, P . S ., N . R . Aluru, and A . C . Cangellaris . "A Methodology for Fast Finite Element Modeling of Electrostatically Actuated MEMS ." International Journal for Numerical Methods in Engineering . Vol 77, Number 13, 2009, 1789 - 1808.

A micro-electro-mechanical (MEM) device typically consist of thin, movable beams or electrodes suspended over a fixed electrode.

Mesh

bmesh = ToBoundaryMesh[];
mesh=ToElementMesh[bmesh,"RegionHoles"->{{0,3*^-6},{4*^-5,-1*^-6}},MaxCellMeasure->{"Length"->5*^-6}];
Show[mesh["Wireframe"],PlotRange->{{-0.85*^-4,0.85*^-4},{-1.7*^-5,1.3*^-5}}]

Mesh

bmesh = ToBoundaryMesh[];
mesh=ToElementMesh[bmesh,"RegionHoles"->{{0,3*^-6},{4*^-5,-1*^-6}},MaxCellMeasure->{"Length"->5*^-6}];
Show[mesh["Wireframe"],PlotRange->{{-0.85*^-4,0.85*^-4},{-1.7*^-5,1.3*^-5}}]

Variables and parameters

vars={V[x,y],{x,y}};
pars=<||>;

Mesh

bmesh = ToBoundaryMesh[];
mesh=ToElementMesh[bmesh,"RegionHoles"->{{0,3*^-6},{4*^-5,-1*^-6}},MaxCellMeasure->{"Length"->5*^-6}];
Show[mesh["Wireframe"],PlotRange->{{-0.85*^-4,0.85*^-4},{-1.7*^-5,1.3*^-5}}]

Variables and parameters

vars={V[x,y],{x,y}};
pars=<||>;

PDEComponent

op=ElectrostaticPDEComponent[vars,pars]

Boundary conditions

bmesh["Wireframe"["MeshElementMarkerStyle"->Blue,"MeshElementStyle"->{Blue,Green,Red,Yellow}]]

Boundary conditions

bmesh["Wireframe"["MeshElementMarkerStyle"->Blue,"MeshElementStyle"->{Blue,Green,Red,Yellow}]]

(Subscript[Γ, p])={ElectricPotentialCondition[ElementMarker==1||ElementMarker==2,vars,pars,<|"ElectricPotential"->0|>],ElectricPotentialCondition[ElementMarker==3,vars,pars,<|"ElectricPotential"->20|>]};

Solve

Vfun=NDSolveValue[{op==0,(Subscript[Γ, p])},V,{x,y}∈mesh];

Legended[DensityPlot[Vfun[x,y],{x,y}∈mesh,],BarLegend[]]

Solid mechanics coupling

Solve electrostatic model

Compute pressure

P=(1/(2ϵ))((|D|)^2)

Solid mechanics coupling

Solve electrostatic model

Compute pressure

P=(1/(2ϵ))((|D|)^2)

Compute deformation

Update electric flux density

|(Subscript[D, def])(x)|=|(Subscript[D, 0])(x)|(G/(G-v(X)))

Repeat  2-4 until convergence

Geometry and mesh

bmeshSM=ToBoundaryMesh["Coordinates"->{{-7.5*^-5,2*^-6},{7.5*^-5,2*^-6},{7.5*^-5,4*^-6},{-7.5*^-5,4*^-6}},"BoundaryElements"->{LineElement[{{1,2},{2,3},{3,4},{4,1}}]}];
meshSM=ToElementMesh[bmeshSM,MaxCellMeasure->{"Length"->1*^-6}];
bmeshSM["Wireframe"]

Geometry and mesh

bmeshSM=ToBoundaryMesh["Coordinates"->{{-7.5*^-5,2*^-6},{7.5*^-5,2*^-6},{7.5*^-5,4*^-6},{-7.5*^-5,4*^-6}},"BoundaryElements"->{LineElement[{{1,2},{2,3},{3,4},{4,1}}]}];
meshSM=ToElementMesh[bmeshSM,MaxCellMeasure->{"Length"->1*^-6}];
bmeshSM["Wireframe"]

Variables and parameters

varsSM={{u[x,y],v[x,y]},{x,y}};
pars["SolidMechanicsModelForm"]="PlaneStress";
pars["YoungModulus"]=170*^9;
pars["PoissonRatio"]=0.34;
pars["Thickness"]=1;

Boundary conditions

(Subscript[Γ, wall])=SolidFixedCondition[ElementMarker==4 ||ElementMarker==6||ElementMarker== 8,varsSM,pars];

(Subscript[Γ, load])=SolidBoundaryLoadValue[ElementMarker==1|| ElementMarker==2||ElementMarker==3,varsSM,pars,<|"Pressure"->{Indexed[BoundaryUnitNormal[x,y],1],Indexed[BoundaryUnitNormal[x,y],2]}*electrostaticPressure[Vfun,0,{x,y}]|>];

PDE

pdeSM={SolidMechanicsPDEComponent[varsSM,pars]==(Subscript[Γ, load]),(Subscript[Γ, wall])};

PDE

pdeSM={SolidMechanicsPDEComponent[varsSM,pars]==(Subscript[Γ, load]),(Subscript[Γ, wall])};

Solve

displacement=NDSolveValue[pdeSM,{u[x,y],v[x,y]},{x,y}∈meshSM];

Iteration

maxdisp=Max[Abs[displacement[[2]][[0]]["ValuesOnGrid"]]];
olddisp=0;

While[Abs[(maxdisp-olddisp)/maxdisp]>1*^-9,
(Subscript[Γ, load])=SolidBoundaryLoadValue[ElementMarker==1|| ElementMarker==2||ElementMarker==3,varsSM,pars,<|"Pressure"->{Indexed[BoundaryUnitNormal[x,y],1],Indexed[BoundaryUnitNormal[x,y],2]}*electrostaticPressure[Vfun,displacement,{x,y}]|>];
displacement=NDSolveValue[{SolidMechanicsPDEComponent[varsSM,pars]==(Subscript[Γ, load]),(Subscript[Γ, wall])},{u[x,y],v[x,y]},{x,y}∈meshSM];
olddisp=maxdisp;
maxdisp=Max[Abs[displacement[[2]][[0]]["ValuesOnGrid"]]];
]

Post - processing

deformedMesh=ElementMeshDeformation[meshSM,displacement,"ScalingFactor"->10];
strain=SolidMechanicsStrain[varsSM,pars,displacement];
stress=SolidMechanicsStress[varsSM,pars,strain];
vonMisesStress=PDEModels`VonMisesStress[varsSM,pars,stress];
deformedvonMisesStress=ElementMeshInterpolation[deformedMesh,vonMisesStress[[0]]["ValuesOnGrid"]];
Show[bmesh["Wireframe"],DensityPlot[deformedvonMisesStress[x,y],{x,y}∈deformedMesh,PlotRange->All,PlotLegends->BarLegend[Automatic,LegendLabel->"[Pa]"]],AspectRatio->(1/4),PlotRange->{{-0.000085`,0.000085`},{-0.000017`,0.000013`}},ImageSize->Large,PlotLabel->Style["Von Mises strees",Bold,15]]

Future work:

Electric currents

Magnetic Fields

Electromagnetic Fields

Questions?

Needs["NDSolve`FEM`"]

electrostaticPressure[potential_InterpolatingFunction,displacement_,{x_,y_}]:=
Module[{Efield,absEfield,eps0,absD0,absD,vdisp,G,ePressure},
Efield=-Grad[potential[x,y],{x,y}];
absEfield=Sqrt[Efield[[1]]^2+Efield[[2]]^2];
eps0=1/((4 Pi*10^-7)*(299792458)^2);
absD0=eps0*absEfield;
absD=absD0;
If[displacement=!=0,
vdisp=Abs[displacement[[2]]];
G=2*^-6;
absD=(absD0*G)/(G-vdisp);
];
ePressure=(absD^2)/(2*eps0)
]