Skip to content
Permalink
Branch: master
Find file Copy path
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
305 lines (243 sloc) 14.6 KB
// This code simulates the harmonic mechanical deflection of a 3D PMUT while supposing axisymmetry.
// A PMUT is a microscale ultrasonic transducer with a piezoelectric actuation. In this example
// the PMUT vibrates in air but it is straightforward to use another medium (e.g. water).
//
// The PMUT geometry and mesh in this example are created using sparselizard but GMSH could have been used instead.
// In any case the goal is not to illustrate the geometry and mesh creation but rather the PMUT axisymmetric simulation.
// The PMUT is made up of (from bottom layer to top):
//
// - a substrate (not in this example)
// - a vacuum cavity
// - a membrane layer
// - a bottom metal layer
// - a piezo layer (sandwiched between the two metal layers)
// - a top metal layer
// - a fluid
//
// Refer to the following paper for more detail:
//
// "Display Compatible pMUT Device for Mid Air Ultrasound Gesture Recognition", C.H. Huang et al.
//
// Please note that since this PMUT is being fabricated the materials have been modified for confidentiality reasons.
// Standard literature properties or arbitrary values have been considered for the material properties.
// It should thus be no surprise that the resonance frequencies and deflections do not match the paper values.
#include "sparselizardbase.h"
using namespace mathop;
// Arguments are:
//
// PMUT radius, radius of the fluid region around (should be large enough), thickness of the top metal, piezo layer and
// bottom metal, thickness of the membrane, thickness of the cavity, electrode coverage (in %) and x-length of the pillar under the membrane.
//
mesh createmesh(double r, double rfluid, double thtopelec, double thpiezo, double thbotelec, double thmem, double thcav, double cov, double lpillar);
void sparselizard(void)
{
wallclock clk;
// Define the PMUT geometric dimensions [m]:
double r = 300e-6, rfluid = 6e-3, thtopelec = 100e-9, thpiezo = 500e-9, thbotelec = 100e-9, thmem = 15e-6, thcav = 35e-6, cov = 0.67, lpillar = 300e-6;
// Axisymmetric assumption:
setaxisymmetry();
// The domain regions as defined in 'createmesh':
int piezo = 1, membrane = 2, topelec = 3, botelec = 4, pillar = 5, fluid = 6, clamp = 7, yaxis = 8, electrode = 9, fluidboundary = 10;
// Create the geometry and the mesh:
mesh mymesh = createmesh(r, rfluid, thtopelec, thpiezo, thbotelec, thmem, thcav, cov, lpillar);
// Write the mesh for display:
mymesh.write("pmutaxisym.msh");
// Define additional regions:
int ground = regionintersection({piezo, botelec});
int solid = regionunion({pillar, membrane, botelec, piezo, topelec});
int isotropicsolid = regionunion({pillar, membrane, botelec, topelec});
int pmuttop = regionintersection({topelec,fluid});
// Harmonic analysis. Set the fundamental frequency [Hz]:
setfundamentalfrequency(163.2e3);
// Nodal shape functions 'h1' for v (the electric potential), p (acoustic
// pressure) and u (membrane displacement). Three components are used for u.
// Use harmonic 2 and 3, i.e. u(x,t) = Us(x)*sin(2pi*f0*t) + Uc(x)*cos(2pi*f0*t)
// for u and p and harmonic 2, i.e. v(x,t) = V(x)*sin(2pi*f0*t) for v.
// The y coordinate field is also defined since it might be of interest.
//
field u("h1xyz",{2,3}), p("h1", {2,3}), v("h1",{2}), y("y");
// Use interpolation order 2 for p and u and 1 for v:
u.setorder(solid, 2);
p.setorder(fluid, 2);
v.setorder(piezo, 1);
// Clamp on the clamp region:
u.setconstraint(clamp);
// To avoid the need to refine too much x-wise close to the y axis:
u.compx().setconstraint(yaxis);
v.setconstraint(electrode, 1);
v.setconstraint(ground, 0);
// Young's modulus [Pa], Poisson's ratio [] (for isotropic materials) and the density [kg/m^3]:
parameter E, nu, rho;
E|pillar = 3.1e9;
E|membrane = 3.1e9;
E|botelec = 100e9;
E|topelec = 100e9;
nu|pillar = 0.34;
nu|membrane = 0.34;
nu|botelec = 0.3;
nu|topelec = 0.3;
rho|pillar = 1300;
rho|membrane = 1300;
rho|botelec = 1000;
rho|piezo = 1780;
rho|topelec = 1000;
rho|fluid = 1.2;
// Propagation speed c [m/s] and a scaling factor for numerical conditionning:
double c = 340, scaling = 1e10;
// Diagonal relative permittivity matrix for the piezo:
expression K(3,3,{7.4,9.3,7.6});
K = K * 8.854e-12;
// Coupling matrix [C/m^2] in Voigt notation for the piezo (6 rows, 3 columns).
// The matrix is expressed in the usual coordinates, i.e. with z being the layer growth direction.
expression C(6,3,{0,0,0.024, 0,0,0.001, 0,0,-0.027, 0,0,0, 0,0,0, 0,0,0});
// To take into account that here y is the growth direction y and z must be flipped:
C.reorderrows({0,2,1,3,5,4}); C.reordercolumns({0,2,1});
// Hooke's matrix [Pa] in Voigt notation for the piezo.
// The matrix is expressed in the usual coordinates, i.e. with z being the layer growth direction.
expression H(6,6, {3.8e9, 1.9e9,3.8e9, 0.9e9,0.9e9,1.2e9, 0,0,0,7e8, 0,0,0,0,9e8, 0,0,0,0,0,9e8});
// To take into account that here y is the growth direction y and z must be flipped:
H.reorderrows({0,2,1,3,5,4}); H.reordercolumns({0,2,1,3,5,4});
formulation pmutmodel;
// Standard isotropic elasticity. dof(u) is the unknown field, tf(u) the test function field.
pmutmodel += integral(isotropicsolid, predefinedelasticity(dof(u), tf(u), E, nu) );
// All inertia terms:
pmutmodel += integral(solid, -rho*dtdt(dof(u))*tf(u) );
// The classical weak formulation for piezoelectricity below can be found e.g. in paper:
//
// "A new fnite-element formulation for electromechanical boundary value problems"
// Define the mechanical equations of the problem in the piezo.
// strain(u) returns the strain vector [exx,eyy,ezz,gyz,gxz,gxy] of u.
pmutmodel += integral(piezo, -( H*strain(dof(u)) )*strain(tf(u)) -( C*grad(dof(v)) )*strain(tf(u)) );
// Define the electrical equations of the problem in the piezo:
pmutmodel += integral(piezo, ( transpose(C)*strain(dof(u)) )*grad(tf(v)) - ( K*grad(dof(v)) )*grad(tf(v)) );
// The wave equation is solved in the fluid:
pmutmodel += integral(fluid, -grad(dof(p))*grad(tf(p)) -1/pow(c,2)*dtdt(dof(p))*tf(p));
// A Sommerfeld condition is used on the fluid boundary to have outgoing waves:
pmutmodel += integral(fluidboundary, -1/c*dt(dof(p))*tf(p));
// Elastoacoustic coupling terms.
// The first term is the forces applied by the fluid on the PMUT top.
// The second term is Newton's law: a membrane acceleration creates a
// pressure gradient in the fluid.
//
// To have a good matrix conditionning the pressure source is divided by
// the scaling factor and to compensate it multiplies the pressure force.
// This leads to the correct membrane deflection but a pressure field divided by the scaling factor.
pmutmodel += integral(pmuttop, -dof(p)*tf(compy(u)) * scaling);
pmutmodel += integral(pmuttop, rho*dtdt(dof(compy(u)))*tf(p) / scaling);
// Generate the algebraic matrix A and right handside vector b of 'Ax = b':
pmutmodel.generate();
// Solve 'Ax = b':
vec sol = solve(pmutmodel.A(), pmutmodel.b());
// Transfer the data from the solution vector to the v, u and p fields:
u.setdata(solid, sol);
p.setdata(fluid, sol);
v.setdata(piezo, sol);
// Write the deflection, pressure and electric potential to file (with an order 2 interpolation for p and u).
u.write(solid, "u.pos", 2);
(scaling*p).write(fluid, "p.pos", 2);
v.write(piezo, "v.pos", 1);
// Write p with 50 timesteps for illustration:
(scaling*p).write(fluid, "p.pos", 2, 50);
// Output the peak deflection:
std::cout << "Peak in-phase deflection: " << 1e9*abs(u.compy().harmonic(2)).max(solid, 6)[0] << " nm" << std::endl;
std::cout << "Peak quadrature deflection: " << 1e9*abs(u.compy().harmonic(3)).max(solid, 6)[0] << " nm" << std::endl;
// Output the pressure at 'rfluid' meters above the PMUT center:
double pressureabove = sqrt( pow(scaling*p.harmonic(2),2) + pow(scaling*p.harmonic(3),2) ).interpolate(fluid, {0,rfluid,0})[0];
std::cout << "Pressure at " << 1000*rfluid << " mm above PMUT center: " << pressureabove << " Pa" << std::endl;
clk.print("Total computation time:");
// Code validation line. Can be removed.
std::cout << (pressureabove < 0.577405 && pressureabove > 0.577404);
}
// THE MESH BELOW IS FULLY STRUCTURED AND IS CREATED USING THE (BASIC) SPARSELIZARD GEOMETRY CREATION TOOL.
// THE ADVANTAGE OF IT IS THAT THE CODE ABOVE CAN BE CALLED FOR ANY PMUT DIMENSION WITHOUT NEEDING CALLS TO EXTERNAL MESHING SOFTWARE.
// AS AN ALTERNATIVE, GMSH COULD HAVE BEEN USED TO EASILY DEFINE THE GEOMETRY AND CREATE A DELAUNAY MESH IN THE FLUID.
mesh createmesh(double r, double rfluid, double thtopelec, double thpiezo, double thbotelec, double thmem, double thcav, double cov, double lpillar)
{
// Give names to the physical region numbers:
int piezo = 1, membrane = 2, topelec = 3, botelec = 4, pillar = 5, fluid = 6, clamp = 7, yaxis = 8, electrode = 9, fluidboundary = 10;
// Number of mesh layers:
int nxpmut = 20, nzthick = 10, nzthin = 5, nair = 50;
int nxpillar = nxpmut*(r+lpillar)/r;
// Cavity layer:
double h = -(thcav+thmem+thbotelec+thpiezo+thtopelec);
shape q13("quadrangle", pillar, {r,h,0, r+lpillar,h,0, r+lpillar,h+thcav,0, r,h+thcav,0}, {nxpillar, nzthick, nxpillar, nzthick});
shape clampline1 = q13.getsons()[0];
clampline1.setphysicalregion(clamp);
shape clampline2 = q13.getsons()[1];
clampline2.setphysicalregion(clamp);
// Membrane layer:
h = h+thcav;
shape q21("quadrangle", membrane, {0,h,0, r*cov,h,0, r*cov,h+thmem,0, 0,h+thmem,0}, {int(nxpmut*cov), nzthick, int(nxpmut*cov), nzthick});
shape q22("quadrangle", membrane, {r*cov,h,0, r,h,0, r,h+thmem,0, r*cov,h+thmem,0}, {int(nxpmut*(1-cov)), nzthick, int(nxpmut*(1-cov)), nzthick});
shape q23("quadrangle", membrane, {r,h,0, r+lpillar,h,0, r+lpillar,h+thmem,0, r,h+thmem,0}, {nxpillar, nzthick, nxpillar, nzthick});
// Bottom electrode layer:
h = h+thmem;
shape q31("quadrangle", botelec, {0,h,0, r*cov,h,0, r*cov,h+thbotelec,0, 0,h+thbotelec,0}, {int(nxpmut*cov), nzthin, int(nxpmut*cov), nzthin});
shape q32("quadrangle", botelec, {r*cov,h,0, r,h,0, r,h+thbotelec,0, r*cov,h+thbotelec,0}, {int(nxpmut*(1-cov)), nzthin, int(nxpmut*(1-cov)), nzthin});
shape q33("quadrangle", botelec, {r,h,0, r+lpillar,h,0, r+lpillar,h+thbotelec,0, r,h+thbotelec,0}, {nxpillar, nzthin, nxpillar, nzthin});
// Piezo layer:
h = h+thbotelec;
shape q41("quadrangle", piezo, {0,h,0, r*cov,h,0, r*cov,h+thpiezo,0, 0,h+thpiezo,0}, {int(nxpmut*cov), nzthin, int(nxpmut*cov), nzthin});
shape q42("quadrangle", piezo, {r*cov,h,0, r,h,0, r,h+thpiezo,0, r*cov,h+thpiezo,0}, {int(nxpmut*(1-cov)), nzthin, int(nxpmut*(1-cov)), nzthin});
shape q43("quadrangle", piezo, {r,h,0, r+lpillar,h,0, r+lpillar,h+thpiezo,0, r,h+thpiezo,0}, {nxpillar, nzthin, nxpillar, nzthin});
// Top electrode layer:
h = h+thpiezo;
shape q51("quadrangle", topelec, {0,h,0, r*cov,h,0, r*cov,h+thtopelec,0, 0,h+thtopelec,0}, {int(nxpmut*cov), nzthin, int(nxpmut*cov), nzthin});
shape q52("quadrangle", topelec, {r*cov,h,0, r,h,0, r,h+thtopelec,0, r*cov,h+thtopelec,0}, {int(nxpmut*(1-cov)), nzthin, int(nxpmut*(1-cov)), nzthin});
shape q53("quadrangle", topelec, {r,h,0, r+lpillar,h,0, r+lpillar,h+thtopelec,0, r,h+thtopelec,0}, {nxpillar, nzthin, nxpillar, nzthin});
shape electrodeline = q51.getsons()[0];
electrodeline.setphysicalregion(electrode);
// Fluid:
shape q61("quadrangle", fluid, {0,0,0, r*cov,0,0, r*cov,r+lpillar,0, 0,r+lpillar,0}, {int(nxpmut*cov), int(nxpmut*(r+lpillar)/r), int(nxpmut*cov), int(nxpmut*(r+lpillar)/r)});
shape q62("quadrangle", fluid, {r*cov,0,0, r,0,0, r,r+lpillar,0, r*cov,r+lpillar,0}, {int(nxpmut*(1-cov)), int(nxpmut*(r+lpillar)/r), int(nxpmut*(1-cov)), int(nxpmut*(r+lpillar)/r)});
shape q63("quadrangle", fluid, {r,0,0, r+lpillar,0,0, r+lpillar,r+lpillar,0, r,r+lpillar,0}, {nxpillar, int(nxpmut*(r+lpillar)/r), nxpillar, int(nxpmut*(r+lpillar)/r)});
shape l1 = q61.getsons()[2];
shape l2("line", -1, {r*cov,r+lpillar,0, r*cov,sqrt(pow(rfluid,2)-pow(r*cov,2)),0}, nair);
shape l3("arc", -1, {r*cov,sqrt(pow(rfluid,2)-pow(r*cov,2)),0, 0,rfluid,0, 0,0,0}, int(nxpmut*cov));
shape l4("line", -1, {0,r+lpillar,0, 0,rfluid,0}, nair);
shape q71("quadrangle", fluid, {l1,l2,l3,l4});
shape l5 = q62.getsons()[2];
shape l6("line", -1, {r,r+lpillar,0, r,sqrt(pow(rfluid,2)-pow(r,2)),0}, nair);
shape l7("arc", -1, {r,sqrt(pow(rfluid,2)-pow(r,2)),0, r*cov,sqrt(pow(rfluid,2)-pow(r*cov,2)),0, 0,0,0}, int(nxpmut*(1-cov)));
shape q72("quadrangle", fluid, {l5,l6,l7,l2});
shape l8 = q63.getsons()[2];
shape l9("line", -1, {r+lpillar,r+lpillar,0, rfluid*sqrt(2)/2,rfluid*sqrt(2)/2,0}, nair);
shape l10("arc", -1, {rfluid*sqrt(2)/2,rfluid*sqrt(2)/2,0, r,sqrt(pow(rfluid,2)-pow(r,2)),0, 0,0,0}, nxpillar);
shape q73("quadrangle", fluid, {l8,l9,l10,l6});
shape l11("line", -1, {r+lpillar,0,0, rfluid,0,0}, nair);
shape l12("arc", -1, {rfluid,0,0, rfluid*sqrt(2)/2,rfluid*sqrt(2)/2,0, 0,0,0}, int(nxpmut*(r+lpillar)/r));
shape l13 = q63.getsons()[1];
shape q74("quadrangle", fluid, {l11,l12,l9,l13});
// Z axis:
shape yaxis2 = q21.getsons()[3];
yaxis2.setphysicalregion(yaxis);
shape yaxis3 = q31.getsons()[3];
yaxis3.setphysicalregion(yaxis);
shape yaxis4 = q41.getsons()[3];
yaxis4.setphysicalregion(yaxis);
shape yaxis5 = q51.getsons()[3];
yaxis5.setphysicalregion(yaxis);
shape yaxis6 = q61.getsons()[3];
yaxis6.setphysicalregion(yaxis);
shape yaxis7 = q71.getsons()[3];
yaxis7.setphysicalregion(yaxis);
// Fluid boundary:
shape fb1 = q71.getsons()[2];
fb1.setphysicalregion(fluidboundary);
shape fb2 = q72.getsons()[2];
fb2.setphysicalregion(fluidboundary);
shape fb3 = q73.getsons()[2];
fb3.setphysicalregion(fluidboundary);
shape fb4 = q74.getsons()[1];
fb4.setphysicalregion(fluidboundary);
// Provide to the mesh all shapes of interest:
mesh mymesh({q13,q21,q22,q23,q31,q32,q33,q41,q42,q43,q51,q52,q53,q61,q62,q63,q71,q72,q73,q74, clampline1,clampline2, yaxis2,yaxis3,yaxis4,yaxis5,yaxis6,yaxis7, electrodeline, fb1,fb2,fb3,fb4});
return mymesh;
}
int main(void)
{
SlepcInitialize(0,{},0,0);
sparselizard();
SlepcFinalize();
return 0;
}
You can’t perform that action at this time.