| computeMForcesfromV(nodes, edges, triangles, potential, geomtype, xyvalues)
|
function perfStats = computeMForcesfromV(nodes, edges, triangles, potential, geomtype, xyvalues)
%% Function to compute forces from electrostatic analysis output from PDE
%% toolbox and do the mechanical analysis to get the deformation
% Compute the magnitude of electric field from potential difference
% E = -grad(potential)
[ex,ey] = pdegrad(nodes,triangles,potential);
absE_tri_midpoints = sqrt(ex.^2+ey.^2);
% Interpolate from triangle data midpoints to node data.
absE = pdeprtni(nodes,triangles,absE_tri_midpoints);
% Compute electrostatic pressure
% permittivity of vacuum (8.854*10^-12 farad/meter)
epsilon0 = 8.854e-12;
electricFluxDensity = epsilon0*absE;
ePressure = CalculateElectrostaticPressure(electricFluxDensity);
% Interpolate the data to compute the electrostatic pressure on a uniform
% grid
xmin=min(nodes(1,triangles)); xmax=max(nodes(1,triangles));
ymin=min(nodes(2,triangles)); ymax=max(nodes(2,triangles));
nt=size(triangles,2);
%nxy=ceil(sqrt(nt/2))+1;
nxy=2*ceil(sqrt(nt/2))+1;
x=linspace(xmin,xmax,nxy);
y=linspace(ymin,ymax,nxy);
[xgrid, ygrid] = meshgrid(x,y);
% Electrostatic Pressure on a uniform grid
ePressure_grid = tri2grid(nodes,triangles,ePressure,x,y);
% Potential on a uniform grid
potential_grid = tri2grid(nodes,triangles,potential,x,y);
% Absolute value of Electric field on a uniform grid
absE_grid = tri2grid(nodes,triangles,absE,x,y);
% Electric flux denisty at the grid points on a uniform grid
electricFluxDensity_grid = tri2grid(nodes,triangles,electricFluxDensity,x,y);
zdata = zeros(size(ePressure_grid));
% Plot the data on the uniform grid
figure
subplot(2,2,1);
surf(x,y,zdata,potential_grid);
curAxis = axis;
axis(curAxis(1:4));
title('Potential Difference')
subplot(2,2,2)
surf(x,y,zdata,absE_grid);
axis(curAxis(1:4));
title('Electric Field')
subplot(2,2,3)
surf(x,y,zdata,ePressure_grid);
axis(curAxis(1:4));
title('Electrostatic Pressure')
% Iteration 1
% Calculate loads and boundary conditions for mechanical analysis
[xCoords, yCoords, pForce, BC, nEdges, nodeIdEdges] = CalculateInitialLoadsBCs(geomtype, xyvalues, ...
x, y, xgrid, ygrid, ...
ePressure_grid);
% Perform mechanical analysis
[perfStats output feaPoints feaTriangles] = mems_cantilever_spmode2(xCoords, yCoords, pForce, BC);
% Plot the results
xDisp = output(1:2:end);
yDisp = output(2:2:end);
% Iteration 2
% Calculate subsequent iteration loads and boundary conditions for
% mechanical analysis
[xCoords, yCoords, pForce, BC] = CalculateLoadsBCs(geomtype, feaPoints, feaTriangles, ...
yDisp, electricFluxDensity_grid, xyvalues, ...
nEdges, nodeIdEdges, xgrid, ygrid, 2);
% Perform 2nd iteration of mechanical analysis
[perfStats output feaPoints2 feaTriangles] = mems_cantilever_spmode2(xCoords, yCoords, pForce, BC);
yDisp2 = output(2:2:end);
% Iteration 3
% Calculate subsequent iteration loads and boundary conditions for
% mechanical analysis
[xCoords, yCoords, pForce, BC] = CalculateLoadsBCs(geomtype, feaPoints2, feaTriangles, ...
yDisp2, electricFluxDensity_grid, xyvalues, ...
nEdges, nodeIdEdges, xgrid, ygrid, 3);
% Save the electrostatic forces, x y coordinate and boundary conditions for
% the next iteration to a mat file
save('inputdata.mat', 'xCoords', 'yCoords', 'pForce', 'BC');
% Perform 2nd iteration of mechanical analysis
[perfStats output feaPoints3 feaTriangles] = mems_cantilever_spmode2(xCoords', yCoords', pForce, BC, 0.01);
yDisp3 = output(2:2:end);
xfea = feaPoints(:,1);
yfea = feaPoints(:,2);
dt = DelaunayTri(xfea,yfea);
bdId = convexHull(dt);
xbd1 = xfea(bdId);
ybd1 = yfea(bdId);
yDisp_Bd1 = yDisp(bdId);
xfea = feaPoints2(:,1);
yfea = feaPoints2(:,2);
dt = DelaunayTri(xfea,yfea);
bdId = convexHull(dt);
xbd2 = xfea(bdId);
ybd2 = yfea(bdId);
yDisp_Bd2 = yDisp2(bdId);
xfea = feaPoints3(:,1);
yfea = feaPoints3(:,2);
dt = DelaunayTri(xfea,yfea);
bdId = convexHull(dt);
xbd3 = xfea(bdId);
ybd3 = yfea(bdId);
yDisp_Bd3 = yDisp3(bdId);
figure
plot(xyvalues(1:4,1), xyvalues(1:4,2));
hold on
plot(xbd1,ybd1+yDisp_Bd1, 'Color', 'r');
plot(xbd2,ybd2+yDisp_Bd2, 'Color', 'g');
plot(xbd3,ybd3+yDisp_Bd3, 'Color', 'm');
legend({
'Undeformed'
'Iteration 1'
'Iteration 2'
'Iteration 3'
},'Location','NorthWest');
hold off
figure
plot([max(abs(yDisp))*1e6 max(abs(yDisp2))*1e6 max(abs(yDisp3))*1e6]);
end
%%
function ePressure = CalculateElectrostaticPressure(electricFluxDensity)
% Function to compute the electrostatic pressure
% Electrostatic pressure is given by
% P = absD0^2/2*epsilon
% where: D0 = epsilon*E
% epsilon = relative coefficient of dielectricity * epsilon0
% assuming, relative coefficient of dielectricity = 1
% permittivity of vacuum (8.854*10^-12 farad/meter)
epsilon0 = 8.854e-12;
ePressure = electricFluxDensity.^2/(2*epsilon0);
end
%%
function [xCoords, yCoords, pForce, BC, nEdges, nodeIdEdges] = CalculateInitialLoadsBCs(geomtype, xyvalues, x, y, xgrid, ygrid, ePressure_grid)
% Find the grid points on the geometry
if strcmpi(geomtype, 'rectangle')
%
% (x4,y4) ------------------------ (x3,y3)
% | |
% | |
% | |
% (x1,y1) ------------------------ (x2,y2)
%
x1 = xyvalues(1,1); y1 = xyvalues(1,2);
x2 = xyvalues(2,1); y2 = xyvalues(2,2);
x3 = xyvalues(3,1); y3 = xyvalues(3,2);
x4 = xyvalues(4,1); y4 = xyvalues(4,2);
dx = x(2)-x(1);
dy = y(2)-y(1);
% Find points closest to the edge connecting (x1,y1) and (x2,y2)
xid = find(xgrid >= (x1-0.5*dx) & xgrid <= (x2+0.5*dx));
yid = find(ygrid >= (y1-dy) & ygrid <= y1);
e1id = intersect(xid, yid);
% Find points closest to the edge connecting (x2,y2) and (x3,y3)
xid = find(xgrid >= x2 & xgrid <= (x2+dx));
yid = find(ygrid >= (y2-0.5*dy) & ygrid <= (y3+0.5*dy));
e2id = intersect(xid, yid);
% Find points closest to the edge connecting (x3,y3) and (x4,y4)
xid = find(xgrid >= (x1-0.5*dx) & xgrid <= (x2+0.5*dx));
yid = find(ygrid >= y3 & ygrid <= (y3+dy));
e3id = intersect(xid, yid);
% Plot the Electrostic pressure at the boundary of the rectangle
subplot(2,2,4)
pForce = [];
% Force on the edge connecting (x1,y1) to (x2,y2)
% unit vector normal to this edge
mag = sqrt((x1-x2)^2+(y2-y1)^2);
f1_dir = [y2-y1 x1-x2]./mag;
% x,y coordinates of distributed load
f1_xCoord = xgrid(e1id(1))+0.5*dx:dx:xgrid(e1id(end));
f1_yCoord = y1*ones(size(f1_xCoord));
for idx=1:length(e1id)-1
f1_vals(idx) = 0.5*(ePressure_grid(e1id(idx))+ePressure_grid(e1id(idx+1)))*dx; %*1e-6; % Assuming unit micron thickness of the beam
end
for idx=1:length(f1_xCoord)
pForce(idx).location = [f1_xCoord(idx) f1_yCoord(idx)];
pForce(idx).compX = f1_vals(idx)*f1_dir(1);
pForce(idx).compY = f1_vals(idx)*f1_dir(2);
end
% Force on the edge connecting (x2,y2) to (x3,y3)
% unit vector normal to this edge
mag = sqrt((x3-x2)^2+(y3-y2)^2);
f2_dir = [y3-y2 x2-x3]./mag;
% x,y coordinates of distributed load
f2_yCoord = ygrid(e2id(1))+0.5*dy:dy:ygrid(e2id(end));
f2_xCoord = x2*ones(size(f2_yCoord));
for idx=1:length(e2id)-1
f2_vals(idx) = 0.5*(ePressure_grid(e2id(idx))+ePressure_grid(e2id(idx+1)))*dy; %*1e-6; % Assuming unit micron thickness of the beam
end
% Force on the edge connecting (x3,y3) to (x4,y4)
% unit vector normal to this edge
mag = sqrt((x4-x3)^2+(y4-y3)^2);
f3_dir = [y4-y3 x3-x4]./mag;
% x,y coordinates of distributed load
f3_xCoord = xgrid(e3id(1))+0.5*dx:dx:xgrid(e3id(end));
f3_yCoord = y3*ones(size(f3_xCoord));
for idx=1:length(e3id)-1
f3_vals(idx) = 0.5*(ePressure_grid(e3id(idx))+ePressure_grid(e3id(idx+1)))*dx; %*1e-6; % Assuming unit micron thickness of the beam
end
lp = length(pForce);
for idx=1:length(f3_xCoord)
pForce(lp+idx).location = [f3_xCoord(idx) f3_yCoord(idx)];
pForce(lp+idx).compX = f3_vals(idx)*f3_dir(1);
pForce(lp+idx).compY = f3_vals(idx)*f3_dir(2);
end
% Plot the Electrostatic forces on the boundary
plot3(f1_xCoord,f1_yCoord,f1_vals, 'Color', 'g', 'Marker', '*')
hold on
plot3(f2_xCoord,f2_yCoord,f2_vals, 'Color', 'g', 'Marker', '*')
plot3(f3_xCoord,f3_yCoord,f3_vals, 'Color', 'g', 'Marker', '*')
hold off
% Boundary conditions
BC(1).location = [x1 y1];
BC(1).X = 0;
BC(1).Y = 0;
BC(2).location = [x4 y4];
BC(2).X = 0;
BC(2).Y = 0;
% XY coordinates of the geometry
xCoords = xyvalues(1:4,1);
yCoords = xyvalues(1:4,2);
% Number of free edges
nEdges = 3;
% Node IDs on free edges
nodeIdEdges(1).id = e1id;
nodeIdEdges(2).id = e2id;
nodeIdEdges(3).id = e3id;
end
end
%%
function [xCoords, yCoords, pForce, BC] = CalculateLoadsBCs(geomtype, feaPoints, feaTriangles, yDisp, electricFluxDensity_grid, xyvalues, nEdges, nodeIdEdges, xgrid, ygrid, itr)
xCoords = [];
yCoords = [];
pForce = [];
BC = [];
if strcmpi(geomtype, 'rectangle')
e1id = nodeIdEdges(1).id;
e3id = nodeIdEdges(3).id;
% Find the points of the FEA mesh that lie on the boundary
xfea = feaPoints(:,1);
yfea = feaPoints(:,2);
dt = DelaunayTri(xfea,yfea);
bdId = convexHull(dt);
% Coordinates of FEA points on the boundary
xbd = xfea(bdId);
ybd = yfea(bdId);
yDisp_Bd = yDisp(bdId);
% Find the FEA nodes close to rectangular grid
xgrid1 = xgrid(e1id);
xgrid3 = xgrid(e3id);
dx = 0.5*(abs(xgrid1(2)-xgrid1(1)));
yDisp_grid = [];
xbd_grid = [];
ybd_grid = [];
for idx=1:length(xgrid1)
nx = find(xbd >= (xgrid1(idx)-dx) & xbd <= (xgrid1(idx)+dx));
ny = find(ybd >= (xyvalues(1,2)-0.01*dx) & ybd <= (xyvalues(1,2)+0.01*dx));
n1 = intersect(nx, ny);
n11 = nearestNeighbor(dt, [xgrid1(idx), xyvalues(1,2)]);
mindist = sqrt((xgrid1(idx)-dt.X(n11,1))^2+(xyvalues(1,2)-dt.X(n11,2))^2);
xNearest = dt.X(n11,1);
yNearest = dt.X(n11,2);
yDispNearest = yDisp(n11);
for jdx=1:length(n1)
dist = sqrt((xgrid1(idx)-xbd(n1(jdx)))^2+(xyvalues(1,2)-ybd(n1(jdx)))^2);
if dist < mindist
mindist = dist;
xNearest = xbd(n1(jdx));
yNearest = ybd(n1(jdx));
yDispNearest = yDisp_Bd(n1(jdx));
end
end
yDisp_grid(idx,1) = yDispNearest; %sum(yDisp_Bd(n1))/length(n1);
xbd_grid(end+1) = xNearest;
ybd_grid(end+1) = yNearest;
end
for idx=1:length(xgrid3)
nx = find(dt.X(:,1) >= (xgrid3(idx)-dx) & dt.X(:,1) <= (xgrid3(idx)+dx));
ny = find(dt.X(:,2) >= (xyvalues(3,2)-0.01*dx) & dt.X(:,2) <= (xyvalues(3,2)+0.01*dx));
n3 = intersect(nx, ny);
n33 = nearestNeighbor(dt, [xgrid3(idx), xyvalues(3,2)]);
mindist = sqrt((xgrid3(idx)-dt.X(n33,1))^2+(xyvalues(3,2)-dt.X(n33,2))^2);
xNearest = dt.X(n33,1);
yNearest = dt.X(n33,2);
yDispNearest = yDisp(n33);
for jdx=1:length(n3)
dist = sqrt((xgrid3(idx)-dt.X(n3(jdx),1))^2+(xyvalues(3,2)-dt.X(n3(jdx),2))^2);
if dist < mindist
mindist = dist;
xNearest = dt.X(n3(jdx),1);
yNearest = dt.X(n3(jdx),2);
yDispNearest = yDisp(n3(jdx));
end
end
yDisp_grid(idx,2) = yDispNearest; %sum(yDisp_Bd(n3))/length(n3);
xbd_grid(end+1) = xNearest;
ybd_grid(end+1) = yNearest;
end
% Assemble D0 (electric flux density) at the boundary
D0 = [electricFluxDensity_grid(e1id) electricFluxDensity_grid(e3id)];
% New electrostatic pressure at the boundary
% Ddef = D0*(Gap/Gap-y displacement)
Gap = xyvalues(5,1)*ones(size(yDisp_grid));
Ddef = D0.*(Gap./(Gap-yDisp_grid));
ePressure_Def = CalculateElectrostaticPressure(Ddef);
% XY coordinates of the geometry
if length(xbd) < length(xgrid1)
xCoords = [xbd_grid(1:length(xgrid1)) xbd_grid(end:-1:length(xgrid1)+1)];
ytempCoords = ybd_grid + [yDisp_grid(:,1)' yDisp_grid(:,2)'];
yCoords = [ytempCoords(1:length(xgrid1)) ytempCoords(end:-1:length(xgrid1)+1)];
if (xCoords(end) ~= xCoords(1) || yCoords(end) ~= yCoords(1))
xCoords(end+1) = xCoords(1);
yCoords(end+1) = yCoords(1);
end
else
xCoords = xbd;
yCoords = ybd + yDisp_Bd;
end
figure
plot(xCoords, yCoords);
hold on
% Compute forces
% Force on the edge connecting (x1,y1) to (x2,y2)
% x coordinate of the distributed load
f1_xCoord = xgrid1(1)+dx:2*dx:xgrid1(end);
ytemp = xyvalues(1,2)*ones(1,length(e1id))+yDisp_grid(:,1)';
for idx=1:length(e1id)-1
f1_vals(idx) = 0.5*(ePressure_Def(idx,1)+ePressure_Def(idx,1))*2*dx; %*1e-6; % Assuming unit micron thickness of the beam
nx = find(xCoords >= (f1_xCoord(idx)-2*dx) & xCoords <= (f1_xCoord(idx)+2*dx));
ny = find(yCoords >= (ytemp(idx)-0.02*dx) & yCoords <= (ytemp(idx)+0.02*dx));
n1 = intersect(nx, ny);
if length(n1)>1
xdist = abs(xCoords(n1)-f1_xCoord(idx)*ones(size(xCoords(n1))));
[xdist ix] = sort(xdist);
f1_yCoord(idx) = max([yCoords(n1(ix(2))) yCoords(n1(ix(1)))]);
else
f1_yCoord(idx) = yCoords(n1);
end
% unit vector normal to this edge
mag = sqrt((xgrid1(idx)-xgrid1(idx+1))^2+(ytemp(idx+1)-ytemp(idx))^2);
f1_dir = [ytemp(idx+1)-ytemp(idx) xgrid1(idx)-xgrid1(idx+1)]./mag;
pForce(idx).location = [f1_xCoord(idx) f1_yCoord(idx)];
pForce(idx).compX = f1_vals(idx)*f1_dir(1);
pForce(idx).compY = f1_vals(idx)*f1_dir(2);
end
% Force on the edge connecting (x3,y3) to (x4,y4)
% x coordinate of the distributed load
lp = length(pForce);
f3_xCoord = xgrid3(1)+dx:2*dx:xgrid3(end);
ytemp = xyvalues(3,2)*ones(1,length(e3id))+yDisp_grid(:,2)';
for idx=1:length(e3id)-1
f3_vals(idx) = 0.5*(ePressure_Def(idx,2)+ePressure_Def(idx,2))*2*dx; %*1e-6; % Assuming unit micron thickness of the beam
nx = find(xCoords >= (f3_xCoord(idx)-2*dx) & xCoords <= (f3_xCoord(idx)+2*dx));
ny = find(yCoords >= (ytemp(idx)-0.05*dx) & yCoords <= (ytemp(idx)+0.05*dx));
n1 = intersect(nx, ny);
if length(n1)>1
xdist = abs(xCoords(n1)-f3_xCoord(idx)*ones(size(xCoords(n1))));
[xdist ix] = sort(xdist);
f3_yCoord(idx) = min([yCoords(n1(ix(2))) yCoords(n1(ix(1)))]);
else
f3_yCoord(idx) = yCoords(n1);
end
% unit vector normal to this edge
mag = sqrt((xgrid3(idx)-xgrid3(idx+1))^2+(ytemp(idx+1)-ytemp(idx))^2);
f3_dir = [ytemp(idx+1)-ytemp(idx) xgrid3(idx)-xgrid3(idx+1)]./mag;
pForce(lp+idx).location = [f3_xCoord(idx) f3_yCoord(idx)];
pForce(lp+idx).compX = f3_vals(idx)*f3_dir(1);
pForce(lp+idx).compY = f3_vals(idx)*f3_dir(2);
end
plot(f1_xCoord,f1_yCoord, '.', 'markersize', 10);
plot(f3_xCoord,f3_yCoord, '.', 'markersize', 10);
title(['Deformed shape - Iteration ' num2str(itr)]);
hold off
% Boundary conditions
BC(1).location = [xyvalues(1,1) xyvalues(1,2)];
BC(1).X = 0;
BC(1).Y = 0;
BC(2).location = [xyvalues(4,1) xyvalues(4,2)];
BC(2).X = 0;
BC(2).Y = 0;
end
end
|
|
