main.cpp 4.04 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
// This code simulates the harmonic mechanical deflection of a 3D polysilicon - PZT microbilayer when 
// the PZT (a piezoelectric material) is sandwiched between two electrically actuated electrodes.


#include "sparselizardbase.h"


using namespace mathop;

void sparselizard(void)
{	
    // The domain regions as defined in 'bilayer.geo':
    int pztlayer = 1, polysiliconlayer = 2, electrode = 3, ground = 4, clamp = 6, freeside = 7, elecbox = 8, mecabox = 9;
    
    mesh mymesh("bilayer.msh");
    
Alexandre Halbach's avatar
Alexandre Halbach committed
17 18 19 20 21
    int wholedomain = regionunion({pztlayer, polysiliconlayer});
    
    // Harmonic analysis. Set the fundamental frequency [Hz]:
    setfundamentalfrequency(1e4);
    
22
    // Nodal shape functions 'h1' for v (the electric potential) and u 
Alexandre Halbach's avatar
Alexandre Halbach committed
23 24 25 26
    // (the membrane displacement). Three components are used for u.
    // Use harmonic 2, i.e. u(x,t) = U(x)*sin(2pi*f0*t) for u and 
    // v(x,t) = V(x)*sin(2pi*f0*t) for v.
    //
27
    field u("h1xyz",{2}), v("h1",{2});
Alexandre Halbach's avatar
Alexandre Halbach committed
28
    
29 30
    // Use interpolation order 2 on 'wholedomain' for u and 1 for v in the PZT layer:
    u.setorder(wholedomain, 2);
Alexandre Halbach's avatar
Alexandre Halbach committed
31
    v.setorder(pztlayer, 1);
32 33 34
    
    // Clamp on surface 'sur' (i.e. 0 valued-Dirichlet conditions):
    u.setconstraint(clamp);
Alexandre Halbach's avatar
Alexandre Halbach committed
35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
    // To force a displacement and see the resulting electric potential you could use:
    // u.compx().setconstraint(freeside, 1e-8);
    // u.setconstraint(freeside, array3x1(1e-8,1.5e-8,0.8e-8));
    
    v.setconstraint(electrode, 10);	
    v.setconstraint(ground, 0);	
    
    
    // Young's modulus, Poisson's ratio and the density of polysilicon:
    double E = 169e9, nu = 0.22, rhopolysi = 2320;
    
    // Diagonal relative permittivity matrix for PZT:
    expression K(3,3,{1704,1704,1433});
    K = K * 8.854e-12;
    
    // Coupling matrix [C/m^2] for PZT (6 rows, 3 columns):
    expression C(6,3,{0,0,-6.62, 0,0,-6.62, 0,0,23.24, 0,17.03,0, 17.03,0,0, 0,0,0});
    
    // Anisotropic Hooke's matrix [Pa] for PZT. Ordering is [exx,eyy,ezz,gyz,gxz,gxy] (Voigt form).
    // Only the lower triangular part (top to bottom and left to right) is provided since it is symmetric.
    expression H(6,6, {1.27e11, 8.02e10,1.27e11, 8.46e10,8.46e10,1.17e11, 0,0,0,2.29e10, 0,0,0,0,2.29e10, 0,0,0,0,0,2.34e10});
    
    // PZT density [kg/m^3]:
    double rhopzt = 7500;
    
    
61
    formulation piezo;
Alexandre Halbach's avatar
Alexandre Halbach committed
62 63 64 65 66 67 68 69 70 71 72
    
    // Standard isotropic elasticity in the polysilicon (not piezoelectric):
    piezo += integral(polysiliconlayer, predefinedelasticity(dof(u), tf(u), E, nu) );
    // Inertia term:
    piezo += integral(polysiliconlayer, -rhopolysi*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.
73 74
    // strain(u) returns the strain vector [exx,eyy,ezz,gyz,gxz,gxy] of u.
    piezo += integral(pztlayer, -( H*strain(dof(u)) )*strain(tf(u)) - ( C*grad(dof(v)) )*strain(tf(u)) );
Alexandre Halbach's avatar
Alexandre Halbach committed
75 76 77
    // Inertia term for PZT:
    piezo += integral(pztlayer, -rhopzt*dtdt(dof(u))*tf(u) );
    // Define the electrical equations of the problem in the piezo:
78
    piezo += integral(pztlayer, ( transpose(C)*strain(dof(u)) )*grad(tf(v)) - ( K*grad(dof(v)) )*grad(tf(v)) );
Alexandre Halbach's avatar
Alexandre Halbach committed
79
    
80
    piezo.generate();
Alexandre Halbach's avatar
Alexandre Halbach committed
81
    
82
    vec sol = solve(piezo.A(), piezo.b());
Alexandre Halbach's avatar
Alexandre Halbach committed
83
    
84 85
    // Transfer the data from the solution vector to the v and u fields:
    u.setdata(wholedomain, sol);
Alexandre Halbach's avatar
Alexandre Halbach committed
86 87
    v.setdata(pztlayer, sol);
    
88 89
    // Write the deflection and electric potential on the volume boundary.
    u.write(mecabox, "u.pos", 2);
Alexandre Halbach's avatar
Alexandre Halbach committed
90 91
    v.write(elecbox, "v.pos", 1);
    // Harmonic fields can be saved in time. Use 50 time steps:
92
    u.write(mecabox, "utime.pos", 2, 50);
Alexandre Halbach's avatar
Alexandre Halbach committed
93
    v.write(elecbox, "vtime.pos", 1, 50);
94 95
    
    // Code validation line. Can be removed.
Alexandre Halbach's avatar
Alexandre Halbach committed
96
    std::cout << (compz(u.harmonic(2)).integrate(electrode, 5) < -1.3622e-15 && compz(u.harmonic(2)).integrate(electrode, 5) > -1.3623e-15);
97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117
}

int main(void)
{	
    SlepcInitialize(0,{},0,0);

    sparselizard();

    SlepcFinalize();

    return 0;
}