Computing electrostatic force in pde toolbox

14 views (last 30 days)
dileesh pv on 5 Jan 2021
I want to compute the electrostatic force on a conductor after solving the Poisson equation in the domain. The code to solve the Poisson equation is given below (This is based on the code given in a pde example: Finite Element Analysis of Electrostatically Actuated MEMS Device)
h = 1;
g0 = h;
b = 5*h;
bg = b;
LX =5*bg;
LY = 2*g0+h;
hmax = 0.2;
x1 = [0;(LX-bg)/2;(LX-b)/2;LX/2;(LX+b)/2;(LX+bg)/2;LX]; % X-coord of init mesh
y1 = [0;g0;g0+h/2;g0+h;LY]; % Y-coord of init mesh
rect_domain = [2,6,x1(1), x1(2), x1(6) x1(7), x1(7), x1(1),y1(1),y1(1),y1(1),y1(1),y1(5),y1(5)]'; % electrostatic domain
rect_movable = [3,4,x1(3) x1(5), x1(5), x1(3),y1(2), y1(2),y1(4),y1(4)]';% Conductor region
rect_movable = [rect_movable;zeros(length(rect_domain)-length(rect_movable),1)]; % padding zeros to match with the size of rect_domain
gd = [rect_domain,rect_movable];
ns = char('rect_domain','rect_movable');
ns = ns';
sf = 'rect_domain-(rect_movable)';
dl = decsg(gd,sf,ns);
model = createpde;
geometryFromEdges(model,dl);
pdegplot(model,'EdgeLabels','on','FaceLabels','on')
xlabel('x-coordinate, meters')
ylabel('y-coordinate, meters')
axis([0,LX,0, LY])
axis tight
V0 = 0;
V1 = 10;
applyBoundaryCondition(model,'dirichlet','Edge',[6],'u',V0); % ground electrode
applyBoundaryCondition(model,'dirichlet','Edge',[3,4,8,9],'u',V1); % Conductor
specifyCoefficients(model,'m',0,'d',0,'c',1,'a',0,'f',0);
% Generate mesh
generateMesh(model,'Hmax',hmax);
pdeplot(model)
xlabel('x-coordinate, meters')
ylabel('y-coordinate, meters')
axis([0,LX,0, LY])
axis tight
% Poisson solver
results = solvepde(model);
u = results.NodalSolution;
figure
pdeplot(model,'XYData',results.NodalSolution,...
'ColorMap','jet');
title('Electric Potential');
xlabel('x-coordinate, meters')
ylabel('y-coordinate, meters')
axis([0,LX,0, LY])
In the example cited( Finite Element Analysis of Electrostatically Actuated MEMS Device ), They have used a function handle to compute the electrostatic pressure as shown below;
pressureFcn = @(location,state) - ...
CalculateElectrostaticPressure(results,[],location);
The function used in the function handle has been given as;
function ePressure = ...
CalculateElectrostaticPressure(elecResults,structResults,location)
% Function to compute electrostatic pressure.
% on the movable electrode.
% Inputs:
% elecResults: Electrostatic FEA results
% structResults: Mechanical FEA results (optional)
% location: The x,y coordinate where pressure is obtained
%
% Output:
% ePressure: Electrostatic pressure at location
%
% location.x : The x-coordinate of the points
% location.y : The y-coordinate of the points
xq = location.x;
yq = location.y;
% Compute the magnitude of the electric field from the potential
% difference.
% The permittivity of vacuum is 8.854*10^-12 farad/meter.
epsilon0 = 8.854e-12;
% Compute the magnitude of the electric flux density.
absD0 = epsilon0*absE;
absD = absD0;
% If structResults (deformation) is available,
% update the charge density along the movable electrode.
if ~isempty(structResults)
% Displacement of the movable electrode at position x
intrpDisp = interpolateDisplacement(structResults,xq,yq);
vdisp = abs(intrpDisp.uy);
G = 2e-6; % Gap 2 micron
absD = absD0.*G./(G-vdisp);
end
% Compute the electrostatic pressure.
ePressure = absD.^2/(2*epsilon0);
end
Here, I have two queries;
• What are the arguments (location,state) passed into the function CalculateElectrostaticPressure. I have to compute the electrostatic pressure on the surface of the conductor) - I have little or no clue about how this function handle works. It works perfectly fine in the example!
• How to integrate this pressure to obtain the total force on the conductor surface (edges E3,E4,E8 & E9)?