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 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 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 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 committed
31
    v.setorder(pztlayer, 1);
32 33 34
    
    // Clamp on surface 'sur' (i.e. 0 valued-Dirichlet conditions):
    u.setconstraint(clamp);
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 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 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 committed
79
    
80
    piezo.generate();
Alexandre Halbach committed
81
    
82
    vec sol = solve(piezo.A(), piezo.b());
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 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 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 committed
93
    v.write(elecbox, "vtime.pos", 1, 50);
94 95
    
    // Code validation line. Can be removed.
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;
}