US simulador

// Análisis de disco cilindrico de PZT5A de 2.5cm x 1cm en vibracion libre
// a la mitad de su seccion rectangular
// Marek Moszynski 30.03.2020

// ------ Variables ------
real[int] ff = [72e3, 73e3, 128e3, 129e3, 156e3,
157e3, 164e3, 165e3, 189e3, 190e3]; // frecuencias inspeccionadas

// ------ Geometria -------
real a = 0.025/2; // radio del disco
real l = 0.01; // espesor del disco
int MM = 15, NN = 12; // resolución de la malla
mesh Sp = square(MM, NN, [ax, ly]); // generación de la malla
plot(Sp, wait=1, ps="ff_mesh.eps"); // graficar y guardar malla

// ------ Constantes ------
real V0 = 0, V1 = 1; // potenciales en los electrodos
real rho = 7750; // densidad
real e31 = -5.4, e33 = 15.8, e15 = 12.3, // constantes piezoeléctricas
eps11S = 8.1e-9, eps33S = 7.3e-9, // constantes dieléctricas
c11 = 120e9, c12 = 75.2e9, c13 = 75.1e9,
c33 = 110e9, c44 = 21.1e9; // constantes elásticas

func C = [[c11, c12, c13, 0, 0, -e31], // matriz de "rigidez" y acoplamiento
[c12, c11, c13, 0, 0, -e31],
[c13, c13, c33, 0, 0, -e33],
[0, 0, 0, c44, -e15, 0],
[0, 0, 0, e15,eps11S, 0],
[e31, e31, e33, 0, 0, eps33S]];

// ------ Macro para el operador diferencial ------
macro L(ur,uz,phi) [ dx(ur), ur/x, dy(uz), dy(ur)+dx(uz), -dx(phi), -dy(phi) ] // diff op

// Abrir archivo para guardar impedancias (frecuencia e impedancia)
ofstream impedFile("impedancias.txt");

for(int ii = 0; ii < ff.n; ii++) {
real f0 = ff[ii], w0 = 2pif0;
cout << f0/1e3 << " kHz" << endl;

// ------ Espacio de elementos finitos -------
fespace Vh3(Sp, P1);
Vh3 ur, uz, phi, vr, vz, w; // variables varicionales

// ------ Problema variacional -------
solve Piezo2D([ur, uz, phi], [vr, vz, w])
= int2d(Sp)( x * rho * w0^2 * [vr, vz]' * [ur, uz] )
- int2d(Sp)( x * L(vr, vz, w)' * C * L(ur, uz, phi) )
+ on(1, phi = V0) // BC: electrodo inferior
+ on(3, phi = V1) // BC: electrodo superior
+ on(4, ur = 0); // BC: eje de simetría

// ------ Cálculo de la impedancia eléctrica ------
// En el electrodo superior (frontera 3) la componente normal es (0,1)
// Se calcula la carga Q integrada (considerando la ponderación x para coordenadas axiales)
real Q = int1d(Sp, 3)( x * ( e33dy(uz) - eps33Sdy(phi) ) );
real I = w0 * Q; // corriente (amplitud) en régimen armónico
real Z = (abs(I) > 1e-12 ? V1 / I : 1e12); // evitar división por cero

// Guardar en el archivo: frecuencia [Hz] y impedancia (valor calculado)
impedFile << f0 << " " << Z << "\n";

// ------- Plot de la malla deformed ----------
real c2 = 100000; // coeficiente de escalado para visualización
mesh Sp2 = movemesh(Sp, [x + c2ur, y + c2uz]);
plot(Sp, Sp2, phi, cmm="f0 = " + f0/1e3 + " kHz", fill=true,
ps="ff_" + (f0/1e3) + "kHz.eps");
}

Crea imagenes eps y un archivo "impedancias.txt".

La primer imagen es de la malla utilizada:

Las imagenes subsecuentes son del campo de desplazamientos: