US simulador

Se incluye código en FreeFem++:

// Free vibrations of 2.5cm x 1cm PZT5A cylindrical disc analysed in half of its rectangular cross-section
// Versión COMPLEJA para obtener la impedancia compleja con pérdidas

// Definir la unidad imaginaria y los potenciales como complejos
complex Iunit = 1.0i;
complex V0 = 0, V1 = 1;   // V0 en el electrodo inferior, V1 en el electrodo superior

// ------ Tabla de frecuencias ------
// Vector de 100 frecuencias equiespaciadas entre 50 kHz y 200 kHz
int Nf = 2000; 
real fmin = 20e3, fmax = 500e3;
real[int] ff(Nf);

for (int i = 0; i < Nf; i++) {
    ff[i] = fmin + i * (fmax - fmin) / (Nf - 1);
}


// ------ Geometría -------
real a = 0.025/2, l = 0.01;      // dimensiones del disco
int MM = 15, NN = 12;            // resolución de la malla
mesh Sp = square(MM, NN, [a*x, l*y]);   // generación de la malla
plot(Sp, ps="ff_mesh.eps");             // guardar malla

// ------ Constantes del material ------
// Parámetros mecánicos y piezoeléctricos (reales)
real rho = 7750;                      // densidad
real e31 = -5.4, e33 = 15.8, e15 = 12.3;    // constantes piezoeléctricas
real c11 = 120e9, c12 = 75.2e9, c13 = 75.1e9, c33 = 110e9, c44 = 21.1e9;  // constantes elásticas
// Parámetros dieléctricos (reales)
real eps11S = 8.1e-9, eps33S = 7.3e-9;       // constantes dieléctricas

// Factores de pérdida (adimensionales)
real etam = 0.05;  // 5% de pérdida mecánica
real etad = 0.01;  // 1% de pérdida dieléctrica

// Matriz de rigidez y acoplamiento con pérdidas:
// Se incluyen pérdidas mecánicas en los coeficientes elásticos y pérdidas dieléctricas en los parámetros dieléctricos.
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 ------
// L: operador aplicado a las variables de la solución (sin conjugación)
macro L(ur,uz,phi) [ dx(ur), ur/x, dy(uz), dy(ur)+dx(uz), -dx(phi), -dy(phi) ] //

// LC: misma función pero con conjugación de los términos de prueba
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)) ] //

// Abrir archivo para guardar (frecuencia, Re(Z), Im(Z))
ofstream impedFile("impedancias_complejas.txt");

// Declaramos el espacio de elementos finitos para funciones reales...
fespace Vh3(Sp, P1);
// ... y luego las variables del espacio como funciones complejas:
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 -------
    solve Piezo2D([ur, uz, phi], [vr, vz, w])
         = int2d(Sp)( x * rho * w0^2 * ( ur * conj(vr) + uz * conj(vz) ) )
         - int2d(Sp)( x * ( LC(vr,vz,w)' * C * L(ur,uz,phi) ) )
         + on(1, phi = V0)
         + on(3, phi = V1)
         + on(4, ur = 0);

    // ------ Cálculo de la respuesta eléctrica ------
    complex Q = int1d(Sp, 3)( x * ( e33*dy(uz) - eps33S*dy(phi) ) );
    // La corriente en régimen armónico: I = Iunit * w0 * Q
    complex Ival = Iunit * w0 * Q;
    // Impedancia: Z = V1 / I (con V1 = 1 V)
    complex Z = V1 / Ival;

    impedFile << f0 << " " << real(Z) << " " << imag(Z) << "\n";

    // ------- Visualización (opcional) -------
    real c2 = 100000;
    mesh Sp2 = movemesh(Sp, [x + c2*real(ur), y + c2*real(uz)]);
    // Para plotear, extraemos la parte real de 'phi' en un espacio real
    fespace Vhr(Sp, P1);
    Vhr phir = real(phi);
    plot(Sp, Sp2, phir, cmm="f0 = ", fill=true, ps="ff_complex.eps");
}

El resultado es un archivo llamado impedancias_complejas.txt

Algunos resultados muestran la grafica de deformaciones a las diferentes frecuencias inspeccionadas:

También se realizan las gráficas de impedancias del archivo impedancia_complejas.txt en octave con el siguiente programa:

close all
clear all
clc

% 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

% Graficar los resultados
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);
xlabel('Frecuencia (kHz)');
ylabel('Im(Z) (\Omega)');
title('Parte Imaginaria de la Impedancia');
grid minor;
axis([0, 300, 10, 1e5]);

% Mostrar la gráfica
print -dpng 'impedancia_compleja.png';  % Guardar la gráfica como imagen

A continuación se muestran los resultados de las gráficas: