US simulador

El modelo FEM es el siguiente:

Solo se considera un piezoeléctrico tipo PZT4 modelado como axial simétrico

Se obtiene una buena respuesta, que se muestra a continuación:

El código utilizado en FreeFem++ fue:

// Mejorado basado en diss.pdf y KocbachLundeVestrheim-UiBSciTechRep1999
// Piezoelectric transducer (PZT-4) de 50 mm de diámetro y 2.49 mm de espesor, simulado en vacío
// Se incluye:
// - PZT en aire
// - Electrodos de plata

// 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 2896 frecuencias equiespaciadas entre 1 kHz y 300 kHz
int Nf = 290; // Para efectos de prueba se modifica a 290 frecuencias
real fminima = 10e3, fmaxima = 299.6e3;
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) ------
// PIEZO
// Dominio del piezoeléctrico
// (0,0), (pRad,0), (pRad,pHigh) y (0,pHigh)
real pRad = 0.05/2, pHigh = 0.00251;

// AIRE
// Dominio completo lleno de aire
// (0,0), (aireRadio,0), (aireRadio,aireAlto) y (0,aireAlto)
real aRad = 0.04, aHigh = 0.04;

// ELECTRODO DE PLATA
// Dominio de electrodo
real eRad = 0.05/2, eHigh = 0.000010;

// Piezo
border pB( t = 0, pRad  ){ x = t;        y = -pHigh/2;     label = 15; }
border pR( t = 0, pHigh ){ x = pRad;     y = -pHigh/2 + t; label = 11; }
border pT( t = 0, pRad  ){ x = pRad - t; y = pHigh/2;      label = 14; }
border pL( t = 0, pHigh ){ x = 0;        y = pHigh/2 - t;  label = 6; }
//mesh piezoM = buildmesh( pB(16) + pR(8) + pT(16) + pL(8) );
//plot( piezoM, wait = true, ps = "piezoM.eps");

// Electrodo 1
border e1B( t = 0, eRad  ){ x = t;        y = -pHigh/2 - eHigh;     label = 9; }
border e1R( t = 0, eHigh ){ x = eRad;     y = -pHigh/2 - eHigh + t; label = 10; }
border e1T( t = 0, eRad  ){ x = eRad - t; y = -pHigh/2;             label = 20; }
border e1L( t = 0, eHigh ){ x = 0;        y = -pHigh/2 - t;         label = 7; }
//mesh electrodo1M = buildmesh( e1B(16) + e1R(4) + e1T(16) + e1L(4) );
//plot( electrodo1M, wait = true, ps = "electrodo1M.eps");

// Electrodo 2
border e2B( t = 0, eRad  ){ x = t;        y = pHigh/2 ;            label = 21; }
border e2R( t = 0, eHigh ){ x = eRad;     y = pHigh/2 + t;         label = 12; }
border e2T( t = 0, eRad  ){ x = eRad - t; y = pHigh/2 + eHigh;     label = 13; }
border e2L( t = 0, eHigh ){ x = 0;        y = pHigh/2 + eHigh - t; label = 5; }
mesh electrodo2M = buildmesh( e2B(16) + e2R(4) + e2T(16) + e2L(4) );
plot( electrodo2M, wait = true, ps = "electrodo2M.eps");

// Area completa rellenada de aire
border mtB( t = 0, aRad     ){ x = t;         y = -aHigh/2;     label = 1; }
border mtR( t = 0, aHigh    ){ x = aRad;      y = -aHigh/2 + t; label = 2; }
border mtT( t = 0, aRad     ){ x = aRad - t;  y = aHigh/2;      label = 3; }
border mtL1( t = 0, aHigh/2-pHigh/2-eHigh ){ x = 0; y = aHigh/2 - t;      label = 4; }
border mtL2( t = 0, pHigh+2*eHigh ){ x = 0;  y = pHigh/2+eHigh - t;       label = 22; }
border mtL3( t = 0, aHigh/2-pHigh/2-eHigh ){ x = 0; y = -pHigh/2-eHigh-t; label = 8; }
mesh mTotal = buildmesh( mtB(16) + mtR(16) + mtT(16) + mtL1(8) + e1L(4) + e2L(4) + mtL3(8) + // parte externa
                         e1B(16) + e1R(4) + // Electrodo 1
                         e2R(4) + e2T(16) + // Electrodo 2
                         pB(16) + pR(8) + pT(16) + pL(8) ); // piezo
plot( mTotal, wait = true, ps = "mTotal.eps");

// Se obtienen las regiones adquiriendo un punto interior de la malla
int piezo = mTotal( pRad/2, 0 ).region;
int electrodo1  = mTotal( eRad/2, pHigh/2 + eHigh/2 ).region;
int electrodo2  = mTotal( eRad/2, -pHigh/2 - eHigh/2 ).region;
int aire  = mTotal( aRad/2, 0.95*aHigh ).region;

// ------ Constantes del material para PZT-4 ------
// Valores aproximados basados en diversas fuentes para PZT‑4.
real rho = 7000.0; // densidad [kg/m^3]

/* Constantes piezoeléctricas [C/m^2] */
real e31 = -4.5, e33 = 12.5, e15 = 9.0;  

/* Constantes elásticas [N/m^2] */
real c11 = 13.0e10, c12 = 7.8e10, c13 = 7.8e10, c33 = 12.0e10, c44 = 2.2e10;

/* Constantes dieléctricas a tensión constante;
   se calculan con ϵ0 = 8.8541878176e-12 */
real eps11S = 500*8.8541878176e-12, eps33S = 600*8.8541878176e-12;

/* Factores de pérdida: se usan tan(δ)=0.02 para ambos (mecánico y dieléctrico) */
real etam = 0.002;  // 0.005 pérdida mecánica
real etad = 0.002;  // 0.005 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‑4),
// no se define C66 pues C66 = (C11 - C12)/2.
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) ]];

// ------ Propiedades del aire (región 2) ------  
// Propiedades mecánicas del aire (simulación de medio fluido)  
real rhoAir = 1.2;          // densidad [kg/m^3]  
real EAir = 85;             // módulo de Young [N/m^2] (calculado a partir del bulk modulus del aire)  
real nuAir = 0.4999;        // coeficiente de Poisson, casi 0.5 para simular un fluido  
real muAir = EAir/(2*(1+nuAir));  
real lambdaAir = EAir*nuAir/((1+nuAir)*(1-2*nuAir));  
real C11Air = lambdaAir + 2*muAir;  
real C12Air = lambdaAir;

// Matriz de rigidez para el aire (solo parte mecánica, sin acoplamiento eléctrico)
// Se asume que en la región 2 la contribución eléctrica es nula (conductor ideal)
func CAir = [
  [ C11Air, C12Air, C12Air, 0,    0, 0 ],
  [ C12Air, C11Air, C12Air, 0,    0, 0 ],
  [ C12Air, C12Air, C11Air, 0,    0, 0 ],
  [ 0,     0,     0,     muAir, 0, 0 ],
  [ 0,     0,     0,     0,    0, 0 ],
  [ 0,     0,     0,     0,    0, 0 ]
];

// ------ Propiedades del electrodo de plata ------
// Propiedades mecánicas de la plata
real rhoAg = 8700;  // densidad [kg/m^3]
real EAg = 60e9;     // módulo de Young [N/m^2]
real nuAg = 0.30;    // coeficiente de Poisson
real muAg = EAg/(2*(1+nuAg));
real lambdaAg = EAg*nuAg/((1+nuAg)*(1-2*nuAg));
real C11Ag = lambdaAg + 2*muAg;
real C12Ag = lambdaAg;
real etamAg = 0.001;  // 0.01 pérdida mecánica
real etadAg = 0.001;  // 0.01 pérdida dieléctrica

// ------ 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(mTotal, 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;

    // En la región 2 se modela la propiedad eléctrica mediante una permittividad efectiva:
    // epsAgEff = eps0 - i*sigmaAg/w0
    real sigmaAg = 6.3e6; // Propiedad eléctrica del electrodo: conductividad de la plata [S/m]
    complex epsAgEff = 8.8541878176e-12 - Iunit*sigmaAg/w0;
    
    // Se define la matriz de rigidez y acoplamiento para el electrodo de plata (región 2)
    // Definición de la matriz de rigidez y acoplamiento para el electrodo de plata, incluyendo pérdidas mecánicas
func CAg = [
  [ C11Ag*(1+1.0i*etamAg), C12Ag*(1+1.0i*etamAg), C12Ag*(1+1.0i*etamAg), 0,    0,                          0 ],
  [ C12Ag*(1+1.0i*etamAg), C11Ag*(1+1.0i*etamAg), C12Ag*(1+1.0i*etamAg), 0,    0,                          0 ],
  [ C12Ag*(1+1.0i*etamAg), C12Ag*(1+1.0i*etamAg), C11Ag*(1+1.0i*etamAg), 0,    0,                          0 ],
  [ 0,                     0,                     0,                     muAg*(1+1.0i*etamAg), 0,              0 ],
  [ 0,                     0,                     0,                     0,    epsAgEff*(1+1.0i*etadAg),     0 ],
  [ 0,                     0,                     0,                     0,    0,                          epsAgEff*(1+1.0i*etadAg) ] ];

    // ------ 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( mTotal, region == piezo )(
               2*pi*x*( rho*w0^2*( ur*conj(vr) + uz*conj(vz) )
                        - LC(vr, vz, w)' * C * L(ur, uz, phi) ) )
         + int2d( mTotal, region == aire )(
               2*pi*x*( rhoAir*w0^2*( ur*conj(vr) + uz*conj(vz) )
                        - LC(vr, vz, w)' * CAir * L(ur, uz, phi) ) )
         + int2d( mTotal, region == electrodo1 )(
               2*pi*x*( rhoAg*w0^2*( ur*conj(vr) + uz*conj(vz) )
                        - LC(vr, vz, w)' * CAg * L(ur, uz, phi) ) )
         + int2d( mTotal, region == electrodo2 )(
               2*pi*x*( rhoAg*w0^2*( ur*conj(vr) + uz*conj(vz) )
                        - LC(vr, vz, w)' * CAg * L(ur, uz, phi) ) )
         + on( 15, phi = V0 ) // Potencial V0 en el electrodo inferior
         + on( 14, phi = V1 ) // Potencial V1 en el electrodo superior
         + on( 4, ur = 0 )
         + on( 5, ur = 0 )
         + on( 6, ur = 0 )
         + on( 7, ur = 0 )
         + on( 8, 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 (borde 14), 
    // multiplicando la integral 1D por 2*pi*x.
    complex Q = int1d( mTotal, 14 )( 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 del campo de presión en el agua -------
    fespace Vhr( mTotal, P1 );
    Vhr presion = real( ur + uz );
    real pmin = presion[].min;
    real pmax = presion[].max;
    fespace Vhr2( mTotal, P1 );
    Vhr2 presionNormalizada = ( presion - pmin ) / ( pmax - pmin );
    //plot(mTotal, presionNormalizada, cmm="Presion normalizada en agua para f0 = " + f0, fill=true, value=true, ps="ff_presion_agua.eps");
   
    // ------- Visualización del campo de desplazamiento -------
    real c2 = 1e6;
    //mesh Sp2 = movemesh( mTotal, [ x + c2*real(ur), y + c2*real(uz) ] );
    Vhr phir = real( phi );
    plot( mTotal, phir, cmm="f0 = " + f0, fill=true, value=true, ps="ff_complex_revised.eps" );
}

y el código en Octave para graficar los resultados:

close all
clear all
clc

% Leer archivos de datos reales del piezo1
piezo1 = dlmread('piezo1_data.s1p');


% Leer el archivo de datos
data = dlmread('impedancias_complejas.txt');

% Extraer columnas
frecuencia = data(:, 1);  % Frecuencia en Hz
ReZ = data(:, 2);         % Parte real de la impedancia
ImZ = data(:, 3);         % Parte imaginaria de la impedancia

%%%%%%%%%%%
Z0 = 50;
frec_piezo1 = piezo1(:, 1);
S11_real = piezo1(:, 2);
S11_imag = piezo1(:, 3);

# Crear el número complejo S11
S11 = S11_real + sqrt(-1) * S11_imag;

# Calcular la impedancia a partir de S11
Z = Z0 * (1 + S11) ./ (1 - S11);

figure;
subplot(2,1,1);
plot(frecuencia / 1e3, -ReZ, 'b', 'LineWidth', 2);
xlabel('Frecuencia (kHz)');
ylabel('Re(Z) (\Omega)');
title('Parte Real de la Impedancia');
grid minor;
axis tight;

subplot(2,1,2);
plot(frecuencia / 1e3, ImZ, 'r', 'LineWidth', 2);
xlabel('Frecuencia (kHz)');
ylabel('Im(Z) (\Omega)');
title('Parte Imaginaria de la Impedancia');
grid minor;
axis tight;

figure;
semilogy(frecuencia/1e3, sqrt(ImZ.^2+ReZ.^2), 'r', 'LineWidth', 2);
hold on;
semilogy(frec_piezo1/1e3, abs(Z));
xlabel('Frecuencia (kHz)');
ylabel('Abs(Z) (\Omega)');
title("Magnitud de la Impedancia");
legend("Simulacion FEM", "Medicion");
grid minor;
axis([10, 300, 1, 10^5]);

figure;
plot(frecuencia/1e3, atan2(ImZ, -ReZ)); %/max(unwrap(atan2(ImZ, -ReZ))), 'r', 'LineWidth', 2);
hold on;
plot(frec_piezo1/1e3, angle(Z)); %/max(unwrap(angle(Z))));
xlabel('Frecuencia (kHz)');
ylabel('Fase(Z) (°)');
title("Fase de la Impedancia");
legend("Simulacion FEM", "Medicion");
grid minor;
axis tight;
%axis([10, 300, 1, 10^5]);