2.3.2 Validación

Como se mencionó en otros tutoriales, los valores nodales de deplazamientos y presiones son almacenados en el archivo binario en una matriz UNOD. Esta matriz matrices posee 3 índices:

• número de nodo
• número de grado de libertad (7 campos): $u_x$, $u_y$, $u_z$, $r_{xy}$, $r_{xz}$, $r_{yz}$ y $u_w$
• paso de tiempo guardado

En este ejemplo, leeremos los desplazamientos y velocidades en los nodo para compararlos con los valores analíticos entregados en §2. Por ejemplo para los desplazamientos, podemos emplear el código Matlab siguiente:

clear all
close all

%--------------------------------------------------------------------Test P
%GEFDyn
filename='tut4_parax';
MODEL=gefread([filename,'.in']);
NodeID=feutil('findnode y==0',MODEL);
[Z,I]=sort(MODEL.Node(NodeID,7),'descend'); NodeID=NodeID(I);
UNOD=gefread([filename,'_SAVE UNOD']);
VNOD=gefread([filename,'_SAVE VNOD']);
TIME=gefread([filename,'_SAVE TIME']);

%Theoretical
rho=1600;   E=10e6; cp=sqrt(E/rho);
for step=1:length(TIME)
   I=find(-Z/cp<TIME(step));
   THEOu(1:length(Z),step)=0;
   THEOv(1:length(Z),step)=0;
   if ~isempty(I)
    THEOu(I,step)=0.02/(rho*cp)*(TIME(step)+Z(I)/cp);
    THEOv(I,step)=0.02/(rho*cp);
   end
end

%Plots
figure %Displacement
hold on
Pas_obs=[20 50 100 150];
for pas=1:length(Pas_obs)
    h1=plot(-squeeze(UNOD(NodeID,3,Pas_obs(pas))),Z,'-r')
    h2=plot(THEOu(:,Pas_obs(pas)),Z,'k--')
end
hold off
xlim([0 2.5e-8])
legend([h1 h2],'Computed','Theoretical','Location','SouthEast')
title('Top P wave')
xlabel('u_z [m]')
ylabel('z [m]')
format_figures

y se obtendrá

Figura 3: Perfiles de $u_z$ para $t=0.02, 0.05, 0.10 \mbox{ y } 0.15$ [s]

que corresponden a los perfiles de desplazamiento vertical para 0.02 [s], 0.05 [s], 0.10 [s] y 0.15 [s ] (ordenados de izquierda a la derecha en la Fig.3). Para los dos primero instantes la onda vertical no ha alcanzado el límite inferior de la malla, mientras que para los dos últimos la onda ha llegado al borde inferior y ha sido absorbida por el elemento paraxial (no hay reflexiones). Se recomienda repetir el cálculo anterior eliminando el elemento paraxial inferior en cuyo caso se verán claramente las reflexiones en dicho borde.

En forma completamente análoga se puede chequear también el perfil de velocidades verticales (Fig.4):

figure %velocity
hold on
Pas_obs=[50];
for pas=1:length(Pas_obs)
    h1=plot(-squeeze(VNOD(NodeID,3,Pas_obs(pas))),Z,'-r')
    h2=plot(THEOv(:,Pas_obs(pas)),Z,'k--')
end
text(0,-2,'t=0.05 [s]','FontSize',14)
hold off
%xlim([0 2.5e-8])
legend([h1 h2],'Computed','Theoretical','Location','NorthWest')
title('Top P wave')
xlabel('v_z [m/s]')
ylabel('z [m]')
format_figures

Figura: Perfiles de $\dot{u}_z$ para $t=0.05$ [s]

En este caso si bien la aproximación del perfil de velocidades es de menor calidad, la tendencia observada es la misma. En la parte superior al frente de onda la velocidad es aproximadamente constante, mientras que por debajo del frente de onda la velocidad debería ser nula. La calidad de la aproximación puede ser mejorada empleando un mallado más fino, un paso de tiempo más fino o bien integrando los elementos paraxiales implícitamente manteniendo los demás parámetros (no disponible en esta versión de GEFDyn).