US simulador

Con el siguiente programa se logra mejorar el matching de la grafica de impedancias compleja del piezoeléctrico:

// Mejorado basado en diss.pdf y KocbachLundeVestrheim-UiBSciTechRep1999
// Piezoelectric transducer (PZT-5A) de 49.96 mm de diámetro y 10.02 mm de espesor, simulado en vacío

// Definición de parámetros complejos y potenciales
complex Iunit = 1.0i;
complex V0 = 0, V1 = 1;   // Potencial: V0 en el electrodo inferior, V1 en el electrodo superior

// ------ Tabla de frecuencias ------
// Se generan 1000 frecuencias equiespaciadas entre 1 kHz y 300 kHz
int Nf = 200;
real fminima = 1e3, fmaxima = 300e3;
real[int] ff(Nf);
for (int i = 0; i < Nf; i++) {
    ff[i] = fminima + i*(fmaxima - fminima)/(Nf - 1);
}

// ------ Geometría y mapeo para 3D (modelo axisimétrico) ------
// Parámetros geométricos:
// a: radio del disco (la mitad del diámetro, es decir, 49.96 mm/2)
// l: espesor del disco (10.02 mm)
real a = 0.04996/2, l = 0.01002;

// Se recomienda una malla más refinada para mejorar la precisión en los modos resonantes
// (véase discusiones en diss.pdf). Aquí se usa una malla 30x24.
mesh Sp = square(30, 24, [ a*x, l*y ]);
plot(Sp, ps="ff_mesh_revised.eps");

// ------ Constantes del material ------
// Valores basados en la Tabla 2.1 de diss.pdf para PZT-5A.
real rho = 7750;   // densidad [kg/m^3]
real e31 = -5.4, e33 = 15.8, e15 = 12.3;  // constantes piezoeléctricas [C/m^2]
real c11 = 12.1e10, c12 = 7.54e10, c13 = 7.52e10, c33 = 11.1e10, c44 = 2.11e10; // constantes elásticas [N/m^2]
// Constantes dieléctricas a tensión constante; se calculan con ϵ0=8.8541878176e-12
real eps11S = 916*8.8541878176e-12, eps33S = 830*8.8541878176e-12;

// Factores de pérdida: se usan tan(δ)=0.02 para ambos (mecánico y dieléctrico), acorde a diss.pdf.
real etam = 0.015;  // pérdida mecánica
real etad = 0.01;  // pérdida dieléctrica

// Matriz de rigidez y acoplamiento (con pérdidas)
// En la notación reducida para un material en la clase de simetría 6mm (PZT-5A), no se define C66
// pues C66 = (C11 - C12)/2. Los parámetros se definen directamente.
func C =  [[ c11*(1+1.0i*etam), c12*(1+1.0i*etam), c13*(1+1.0i*etam), 0,                     0,                     -e31 ],
           [ c12*(1+1.0i*etam), c11*(1+1.0i*etam), c13*(1+1.0i*etam), 0,                     0,                     -e31 ],
           [ c13*(1+1.0i*etam), c13*(1+1.0i*etam), c33*(1+1.0i*etam), 0,                     0,                     -e33 ],
           [ 0,                 0,                 0,                 c44*(1+1.0i*etam),     -e15,                 0 ],
           [ 0,                 0,                 0,                 e15,                  eps11S*(1+1.0i*etad), 0 ],
           [ e31,               e31,               e33,               0,                     0,                     eps33S*(1+1.0i*etad) ]];

// ------ Macros para los operadores diferenciales ------
// Se utilizan las macros originales para el caso axisimétrico en FreeFEM++
// En esta formulación, x representa el radio físico (r)
macro L(ur,uz,phi) [ dx(ur), ur/x, dy(uz), dy(ur)+dx(uz), -dx(phi), -dy(phi) ] //
macro LC(vr,vz,w) [ conj(dx(vr)), conj(vr)/x, conj(dy(vz)), conj(dy(vr)+dx(vz)), -conj(dx(w)), -conj(dy(w)) ] //

// Nota: Para obtener la cantidad total en 3D se debe multiplicar la integración en el dominio 2D (y en el de borde)
// por el factor 2π (la integración en la variable angular). Esto es crucial para la normalización.
// Se han modificado las integrales para incluir el factor 2*pi.

// Abrir archivo para guardar (frecuencia, Re(Z), Im(Z)) para análisis de resonancia
ofstream impedFile("impedancias_complejas.txt");

// Espacio de elementos finitos (elementos P1)
fespace Vh3(Sp, P1);
Vh3<complex> ur, uz, phi, vr, vz, w;

for (int ii = 0; ii < ff.n; ii++) {
    real f0 = ff[ii], w0 = 2*pi*f0;
    cout << f0/1e3 << " kHz" << endl;
    
    // ------ Formulación variacional COMPLEJA ------
    // Se utiliza la variable x (que representa el radio físico) y se multiplica por 2*pi para integrar en el dominio completo.
    solve Piezo2D([ur, uz, phi], [vr, vz, w])
         = int2d(Sp)( 
               2*pi*x * ( rho * w0^2 * ( ur * conj(vr) + uz * conj(vz) ) )
             - 2*pi*x * ( LC(vr,vz,w)' * C * L(ur,uz,phi) )
           )
         + on(1, phi = V0)
         + on(3, phi = V1)
         + on(4, ur = 0);  // Se fija el desplazamiento radial en x=0 (borde fijo o eje de revolución)
         
    // ------ Cálculo de la respuesta eléctrica ------
    // Se calcula la carga eléctrica en la interfaz superior, multiplicando la integral 1D por 2*pi*x.
    complex Q = int1d(Sp, 3)( 2*pi*x * ( e33*dy(uz) - eps33S*dy(phi) ) );
    // La corriente en régimen armónico es I = Iunit * w0 * Q y la impedancia se define como Z = V1 / I.
    complex Ival = Iunit * w0 * Q;
    complex Z = V1 / Ival;
    impedFile << f0 << " " << real(Z) << " " << imag(Z) << "\n";
    
    // ------- Visualización (opcional) -------
    // Se genera una malla movida para visualizar el desplazamiento real.
    real c2 = 100;
    mesh Sp2 = movemesh(Sp, [ x + c2*real(ur), y + c2*real(uz) ]);
    fespace Vhr(Sp, P1);
    Vhr phir = real(phi);
    plot(Sp, Sp2, phir, cmm="f0 = ", fill=true, ps="ff_complex_revised.eps");
}