Visualización de Campos Escalares y Vectoriales en elasticidad, Grupo(2)

Trabajo realizado por estudiantes
Título	Visualización de Campos Escalares y Vectoriales en elasticidad, Grupo(2)
Asignatura	Teoría de Campos
Curso	2023-24
Autores	Soto, Pablo Mateo, Nicolás García, Andrea Villaverde, Rodrigo Pérez, Jaime
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

El objetivo de este artículo es visualizar los distintos campos escalares y campos vectoriales en elasticidad. Para ello, consideramos una placa rectangular plana y dos magnitudes físicas, la temperatura, expresada por un campo escalar y el desplazamiento ondulatorio producido por una fuerza determinada. Supondremos que el desplazamiento se trata de una onda longitudinal en la que la dirección de propagación es ortogonal a la amplitud. Haremos uso de programas de cálculo matemático como Octave y MATLAB para realizar cálculos complejos y representar visualmente los distintos campos.

1 Representación de la placa rectangular

Para empezar el ejercicio deberemos crear un mallado acorde a los parámetros establecidos en el enunciado. Para ello deberemos tomar los ejes empleando el comando 'axis' en el rectángulo definido por [math] (x, y) ∈ [-1,1] × [0,12] [/math] y empleando un paso de muestreo [math] h=\frac{2}{10}[/math] para las variables [math] x [/math] e [math] y [/math]. Haremos uso del siguiente código:

clc
clear
h=2/20; %Definimos el paso de muestro
x=[-1:h:1]; %Definimos los límites de la placa
y=[0:h:12];
[X,Y]=meshgrid(x,y); %Generamos el mallado de los parámetros establecidos
Z=0.*X.*Y; 
mesh(X,Y,Z); %Creamos la placa
axis equal %Los ejes mantienen las mismas proporciones
title('Placa rectangular plana')
xlabel('eje x') %Nombramos los ejes
ylabel('eje y')
view(2); %Establecemos el ángulo de vista

2 Curvas de nivel de la temperatura

La temperatura viene definida por la expresión [math] T(x,y)=3log(1+(x+1)^2)+log(1+(y-2)^2) [/math]. Además de esto, deberemos indicar cuales son los puntos de la gráfica en los que la temperatura es máxima por lo que calcularemos el gradiente del campo y lo representaremos graficamente mediante MATLAB.

clc
clear
subplot(1,2,1) %Representamos la gráfica de la primera columna
h=2/20; %Definimos el paso de muestro
x=[-1:h:1]; %Definimos los límites de la placa
y=[0:h:12];
[X,Y]=meshgrid(x,y); %Generamos el mallado de los parámetros establecidos
Z=3*log(1+(X+1).^2)+log(1+(Y-2).^2); %Ponemos la función de la T en función del eje Z
mesh(X,Y,Z) %Creamos la placa
axis equal %Los ejes mantienen las proporciones
title("Curvas de nivel de la temperatura en 3D") %Título de la gráfica y los ejes
xlabel('eje x')
ylabel('eje y')
zlabel('T(x,y)')
subplot(1,2,2) %Representamos la gráfica de la segunda columna
pcolor(X,Y,Z)
shading flat %Quitamos el mallado negro
hold on
contour(X,Y,Z,20,'k') %Creamos las curvas de nivel en color negro
hold off
axis equal %Los ejes mantienen las proporciones
title("Curvas de nivel de la temperatura en 2D") %Título de la gráfica y los ejes
xlabel('eje x')
ylabel('eje y')
zlabel('T(x,y)')

2.1 Temperatura máxima

Según la gráfica, podemos observar que el punto donde más temperatura se alcanza es en el punto [math](1,12)[/math]. Para saber cual es el valor máximo de la temperatura utilizaremos el siguiente comando al final del código anterior, este código nos dará el valor máximo que alcanza la temperatura en la región:

Tmax=max(max(Z))

El programa nos indicará que el valor de la temperatura máxima es 9,4434ºC.

2.2 Cálculo de ▽T y su representación

Para calcular el gradiente de un campo escalar, que en este caso es la temperatura dada en función de [math]T= 3log(1+(x+1)^2)+log(1+(y-2)^2)[/math], deberemos realizar la derivada parcial de cada coordenada reflejada sobre la base ortonormal orientada positivamente [math]\vec{i},\vec{j},\vec{k}[/math], sin embargo, al ser una placa en 2 dimensiones no se realizará la derivada parcial de la coordenada [math]\vec{k}[/math]. Por lo que la divergencia de la temperatura será [math]▽T=\frac{\partial{T}}{\partial{x}}\vec{i}+\frac{\partial{T}}{\partial{y}}\vec{j}[/math], que con nuestra función temperatura es [math]▽T=(\frac{6x+6}{1+(x+1)^2})\vec{i}+(\frac{2y-4}{1+(y-2)^2})\vec{j}[/math]. Una vez calculada la divergencia podemos representarlo en la gráfica. Utilizaremos el siguiente código:

clc
clear
h=2/20;
f=inline('3*log(1+(x+1).^2)+log(1+(y-2).^2)','x','y'); %Creamos las funciones del campo
fx=inline('(6*(x+1))./(1+(x+1).^2)','x','y');
fy=inline('(2*(y-2))./(1+(y-2).^2)','x','y');
x=[-1:h:1]; %Discretizamos la placa
y=[0:h:12];
[X,Y]=meshgrid(x,y); %Generamos el mallado de los parámetros establecidos
Z=f(X,Y); 
contour(X,Y,Z,20,'g') %Creamos las curvas de nivel del campo
hold on
k=0.3; %Establecemos el tamaño de paso
xx=[-1:k:1];
yy=[0:k:12];
[XX,YY]=meshgrid(xx,yy);
U=fx(XX,YY);
V=fy(XX,YY);
quiver(XX,YY,U,V,'k') %Representamos el campo vectorial
axis equal %Los ejes mantienen las proporciones
title("Representación de ▽T")
xlabel('eje x') %Título de la gráfica y los ejes
ylabel('eje y')
zlabel('T(x,y)')

3 Energía calorífica [math] \vec{Q} [/math]

De acuerdo con la Ley de Fourier, la energía calorífica [math] \vec{Q} [/math] viaja según la fórmula [math] \vec{Q}=-k∇T [/math], donde κ es la constante de conductividad térmica de la placa que supondremos [math]k = 1[/math]. Una vez calculada la energía calorífica [math] \vec{Q}= -(\frac{6x+6}{1+(x+1)^2})\vec{i}+(\frac{2y-4}{1+(y-2)^2})\vec{j}[/math] deberemos representarlo en la gráfica utilizando el siguiente código:

clc
clear
h=2/20;
f=inline('3*log(1+(x+1).^2)+log(1+(y-2).^2)','x','y'); %Creamos las funciones del campo
fx=inline('-(6*(x+1))./(1+(x+1).^2)','x','y');
fy=inline('-(2*(y-2))./(1+(y-2).^2)','x','y');
x=[-1:h:1]; %Discretizamos la placa
y=[0:h:12];
[X,Y]=meshgrid(x,y); %Generamos el mallado de los parámetros establecidos
Z=f(X,Y); 
contour(X,Y,Z,20,'g') %Creamos las curvas de nivel del campo
hold on
k=0.3; %Establecemos el tamaño de paso
xx=[-1:k:1];
yy=[0:k:12];
[XX,YY]=meshgrid(xx,yy);
U=fx(XX,YY);
V=fy(XX,YY);
quiver(XX,YY,U,V,'k') %Representamos el campo vectorial
axis equal %Los ejes mantienen las proporciones
title("Representación de ▽T")
xlabel('eje x') %Título de la gráfica y los ejes
ylabel('eje y')
zlabel('T(x,y)')

Se puede observar que el campo generado por la energía calorífica tiene sentido opuesto a ▽T.

4 Campo vectorial para [math]t=0[/math]

El campo vectorial viene definido por la expresión [math] \vec{u}(x,y,t)= \vec{a}sin(πk(\vec{d}*\vec{r_o}(x,y)-vt))[/math]. Sustituyendo los valores de la ecuación indicados en el enunciado, [math] \vec{a}=\vec{d}=\frac{1}{3}\vec{j} [/math], [math]k=1[/math], y por último [math]r_o=x\vec{i}+y\vec{j}[/math], para [math]t=0 [/math] nos queda que [math]\vec{u}(x,y,0)= \frac{1}{3}\vec{j}sin(π(\frac{1}{3}\vec{j}*(x\vec{i}+y\vec{j})))[/math]. Si continuamos operando la expresión nos queda [math]\vec{u}(x,y,0)= (0,\frac{1}{3},0)sin(x,y(0,\frac{1}{3},0)(1,1,0))[/math]. Fianlmente obtenemos que [math]\vec{u}(x,y,0)=\frac{1}{3}sin(\frac{1}{3},x,y)\vec{j}[/math]. Utilizaremos el siguiente código para representar el campo vectorial:

clc
clear
k=0.3; %Establecemos un tamaño de paso
u=inline('0.*x.*y','x','y'); %Aplicamos la función en u
v=inline('(1/3)*sin(pi*y/3)','x','y'); %Aplicamos la función en v
x=[-1:k:1];% Discretización de x
y=[0:k:12]; % Discretización de  y
[X,Y]=meshgrid(x,y); %Generamos el mallado de los parámetros establecidos
xx=[-1:k:1];
yy=[0:k:12];
[XX,YY]=meshgrid(xx,yy); %Creamos el mallado del campo
U=u(XX,YY); %Aplicamos la funcion u a XX e YY
V=v(XX,YY);% Aplicamos la funcion v a XX e YY
quiver(XX,YY,U,V) %Representamos el campo vectorial
axis image %Ajustar los limites con la figura y mantiene la proporción de los ejes
title('Campo Vectorial')
xlabel('eje x')
ylabel('eje y')

5 Representación de las deformaciones de la placa

El campo vectorial [math]\ \vec{u}(x,y,0)=\frac{1}{3}sin(\frac{π}{3}y)\vec{j}[/math] es el campo de fuerzas que actúan sobre nuestra placa en el instante t=0. La representación de las posiciones en los puntos de la placa mediante MATLAB es la siguiente:

clc
clear
k=0.2;
x=[-1:k:1];% Discretización x
y=[0:k:12];% Discretizción y
[X,Y]=meshgrid(x,y); %Generamos el mallado de los parámetros establecidos
u=0.*X.*Y; %Función de u
v=(1/3)*sin(pi*Y./3);% Función de v
Z=zeros(size(X)); %Creamos una matriz de ceros de tamaño X
mesh(X,Y,Z); %Creamos el mallado del campo
subplot(1,2,1)
pcolor(X,Y,Z)
axis equal %Los ejes mantienen las proporciones
xlabel('eje x') %Nombramos los ejes
ylabel('eje y')
title ('antes desplazamiento')
subplot(1,2,2) %en esta figura sumamos los vectores posicion antes de hacer acuar el campo mas los vectores que se generan una vez acutado el campo
XX=X+u; %Suma de u las posiciones de x
YY=Y+v; %Suma de ven las posiciones de y
pcolor(XX,YY,Z);
axis equal
xlabel('eje x')
ylabel('eje y')
title('despues del desplazamiento')

6 Divergencia de un campo de desplazamientos

Dado el campo vectorial [math]\ \vec{u}(x,y,0)=\frac{1}{3}sin(\frac{π}{3}.y)\vec{j}[/math], la divergencia de [math]\ \vec{u}[/math] viene dada por [math]▽\vec{u}=\frac{\partial{u1}}{\partial{x}}+\frac{\partial{u2}}{\partial{y}}[/math], siendo [math] \ u_1=0\ [/math] y [math] \ u_2=\frac{1}{3}sin(\frac{π}{3}y)\ [/math]. Nos queda la siguiente exptresión [math]\ ▽\vec{u}=\frac{π}{9}cos(\frac{π}{3}y)\vec{j}[/math].

La Divergencia[math]\ ▽\vec{u} [/math] nos muestra el flujo del campo por unidad de volumen.

clc
clear
h=2/20;
x=[-1:h:1]; %Discretización x
y=[0:h:12]; %Discretización y
[X,Y]=meshgrid(x,y); %Generamos el mallado de los parámetros establecidos
Z=pi/9*cos((pi/3).*Y); %Función que define la divergencia con dos variables x e y
subplot(1,2,1) %Representación de la primera figura
surf(X,Y,Z); %Representación de la superficie
axis equal
xlabel('x');
ylabel('y')
zlabel('z')
title ('Divergencia')
shading flat
subplot(1,2,2) %Representación de la segunda figura
hold on
axis equal
pcolor(X,Y,Z);
contour(X,Y,Z,10,'K') %Creamos las curvas de nivel
colorbar
shading flat % Quitamos el mallado negro
xlabel('x');
ylabel('y')
zlabel('Divergencia')
axis equal

Como se puede observar en la gráfica, los puntos en los que la divergencia es máxima es en los puntos donde el eje y equivale a 0,6 y 12. Los puntos en los que la divergencia es mínima son aquellos donde el eje y equivale a 3 y 9. Por último, los puntos donde la divergencia es nula equivalen a aquellos puntos donde la gráfica no crece ni decrece, es decir, en los máximos y mínimos. Por otra parte, si se puede observar como la divergencia es un cambio de volumen.

7 Rotacional de un campo de desplazamientos

El rotacional de un campo vectorial [math]\vec{u}(x_1,x_2,x_3)[/math], expresado en coordenadas cartesianas, se define como la tendencia del campo a inducir rotación alrededor de un punto, y lo designamos por[math]\ \bigtriangledown \times \vec{u}[/math] .

Sea [math] \vec{u}(x_1,x_2,x_3)=u_1(x_1,x_2,x_3)\vec{i}+u_2(x_1,x_2,x_3)\vec{j}+u_3(x_1,x_2,x_3)\vec{k} [/math], expresado en la base canónica de [math]ℝ^3[/math], el rotacional de [math]\vec{u}[/math] se define como: [math]\ \bigtriangledown \times \vec{u} = \begin{vmatrix}\vec{i} & \vec{j} &\vec{k} \\ \frac{\partial }{\partial x_1} & \frac{\partial }{\partial x_2} & \frac{\partial }{\partial x_3}\\ u_1 & u_2 & u_3\end{vmatrix}[/math].

Empleando en nuestro caso, el campo vectorial [math]\vec{u}(x,y,z)=\frac{1}{3}sen(π/3y)\vec{j}[/math], el rotacional vendría definido como [math]\ \bigtriangledown \times \vec{u} = \begin{vmatrix}\vec{i} & \vec{j} &\vec{k} \\ \frac{\partial }{\partial x} & \frac{\partial }{\partial y} & \frac{\partial }{\partial 3}\\ 0 & \frac{1}{3}sen(π/3y)\vec{j} & 0\end{vmatrix}=(0-0)\vec{i}+(0-0)\vec{j}+(0-0)\vec{k}=\vec{0}[/math].

Al obtener como solucion del rotacional el vector nulo, podemos ver que nos encontramos con que [math]\vec{u}(x_1,x_2,x_3)[/math] es un campo irrotacional ó conservativo en el que se cumple que [math]\ \left | \bigtriangledown \times \vec{u} \right | = 0 [/math]. Al ser irrotacional no habrá puntos que sufran un mayor rotacional que otros.

8 Tensiones normales en la dirección que marcan los ejes [math]\vec{i}, \vec{j}[/math] y [math]\vec{k}[/math]

Tenemos el campo de tensiones [math]\sigma=λ\bigtriangledown \cdot \vec{u}1 + 2µe[/math] en un medio elástico, isótropo y homogéneo, donde la parte simétrica del gradiente del campo vectorial [math]\vec u[/math] es [math]e=\frac{\bigtriangledown\vec{u}+\bigtriangledown\vec{u}^t}{2}[/math]. Los coeficientes de Lamé son [math]λ=µ=1[/math] y la divergencia del campo vectorial [math]\vec u[/math] es [math]\bigtriangledown \cdot \vec{u}[/math] . Para calcular las deformaciones en las direcciones de los vectores de la base ortonormal orientada positivamente [math]{\vec i,\vec j,\vec k}[/math], se comienza con estos cálculos:

[math]\bigtriangledown\vec{u}=\begin{pmatrix}\frac{\partial u_{1}}{\partial x}&\frac{\partial u_{1}}{\partial y}&\frac{\partial u_{1}}{\partial z}\\\frac{\partial u_{2}}{\partial x}&\frac{\partial u_{2}}{\partial y}&\frac{\partial u_{2}}{\partial z}\\\frac{\partial u_{3}}{\partial x}&\frac{\partial u_{3}}{\partial y}&\frac{\partial u_{3}}{\partial z} \end{pmatrix}=\begin{pmatrix}0&0&0\\0&\frac{\pi}{9}cos(\frac{\pi}{3}y)&0\\0&0&0\end{pmatrix}= \frac{\pi}{9}cos(\frac{\pi}{3}y) \vec j \otimes \vec j[/math]

Podemos ver que [math]\bigtriangledown\vec{u}[/math] es una matriz simétrica, por lo que su matriz traspuesta y ella misma son idénticas y se representarán de la misma forma.

[math]e=\frac{\bigtriangledown\vec{u}+\bigtriangledown\vec{u}^t}{2}=\frac{\frac{\pi}{9}cos(\frac{\pi}{3}y) \vec j \otimes \vec j + \frac{\pi}{9}cos(\frac{\pi}{3}y) \vec j \otimes \vec j}{2}=\frac{\pi}{9}cos(\frac{\pi}{3}y) \vec j \otimes \vec j[/math]

Continuamos calculando el tensor de deformaciones:

[math]\sigma=\frac{\pi}{9}cos(\frac{\pi}{3}y)(\vec i \otimes \vec i + \vec j \otimes \vec j + \vec k \otimes \vec k) + 2\cdot\frac{\pi}{9}cos(\frac{\pi}{3}y)(\vec j \otimes \vec j) = \frac{\pi}{9}cos(\frac{\pi}{3}y) (\vec i \otimes \vec i) +\frac{\pi}{3}cos(\frac{\pi}{3}y) (\vec j \otimes \vec j) +\frac{\pi}{9}cos(\frac{\pi}{3}y) (\vec k \otimes \vec k) [/math]

[math]\sigma_{ij} = \begin{pmatrix} \frac{\pi}{9}cos(\frac{\pi}{3}y)\ & 0 & 0 \\ 0 & \frac{\pi}{3}cos(\frac{\pi}{3}y)\ & 0 \\ 0 & 0 & \frac{\pi}{9}cos(\frac{\pi}{3}y)\ \end{pmatrix}[/math]

Finalmente calculamos las tensiones normales. Estas tensiones son las que se encuentran en los elementos de la diagonal de la matriz [math]\sigma_{ij}[/math]:

En la dirección [math]\vec{i}[/math]: [math]\:\:\vec i \cdot \sigma \cdot \vec i[/math] = [math]\frac{\pi}{9}cos(\frac{\pi}{3}y)[/math]

En la dirección [math]\vec{j}[/math]: [math]\:\:\vec j \cdot \sigma \cdot \vec j[/math] = [math]\frac{\pi}{3}cos(\frac{\pi}{3}y)[/math]

En la dirección [math]\vec{k}[/math]: [math]\:\:\vec k \cdot \sigma \cdot \vec k[/math] = [math]\frac{\pi}{9}cos(\frac{\pi}{3}y)[/math]

Tras los cálculos lo representamos en Matlab utilizando el código:

clc
clear
h=0.2;

%Discretizamos y creamos mallado
x=[-1:h:1];
y=[0:h:12];       
[X,Y]=meshgrid(x,y);  

%Tensiones normales en i,j y k
i=(pi/9).*cos((pi/3).*Y);
j=(pi/3).*cos((pi/3).*Y);       
k=(pi/9).*cos((pi/3).*Y);  

%Creamos las tres gráficas
figure
hold on
subplot(1,3,1)           
pcolor(X,Y,i)
colorbar
shading flat
xlabel('Eje X')
ylabel('Eje Y') 
zlabel('Eje Z')
axis equal
title('Tensiones normales en i')  

subplot(1,3,2) 
pcolor(X,Y,j)
colorbar
shading flat
xlabel('Eje X')
ylabel('Eje Y') 
zlabel('Eje Z')
axis equal
title('Tensiones normales en j')

subplot(1,3,3)
pcolor(X,Y,k)
colorbar
shading flat
xlabel('Eje X')
ylabel('Eje Y') 
zlabel('Eje Z')
axis equal
title('Tensiones normales en k')
hold off

9 Tensiones tangenciales respecto al plano ortogonal a [math]\ \vec i [/math]

Sea [math]\ \left| \sigma ·\vec i -(\vec i·\sigma·\vec i)\vec i \right| [/math] la expresión de la tensión tangencial [math] \sigma_{T} [/math] respecto al plano ortogonal a [math]\vec{i}[/math] en [math]t=0[/math]

[math]\sigma_{T}=\ \left| \begin{pmatrix}\frac{\pi}{9}cos(\frac{\pi}{3}·y)&0&0\\0&\frac{\pi}{3}cos(\frac{\pi}{3}·y)&0\\0&0&\frac{\pi}{9}cos(\frac{\pi}{3}·y)\end{pmatrix}·\begin{pmatrix}1 & 0 & 0\end{pmatrix}-\left(\begin{pmatrix}1 \\0 \\0 \\ \end{pmatrix}·\begin{pmatrix}\frac{\pi}{9}cos(\frac{\pi}{3}·y)&0&0\\0&\frac{\pi}{3}cos(\frac{\pi}{3}·y)&0\\0&0&\frac{\pi}{9}cos(\frac{\pi}{3}·y)\end{pmatrix}·\begin{pmatrix}1 & 0 & 0\end{pmatrix}\right)\vec i \right| [/math]

Obtenemos que:

[math]\sigma_{T} =\ \left|\frac{\pi}{9}cos(\frac{\pi}{3}·y)\vec i - \frac{\pi}{9}cos(\frac{\pi}{3}·y)\vec i \right| = \vec{0}[/math]

Quedando demostrado que las tensiones tangenciales respecto al plano ortogonal a [math]\vec{i}[/math] (que contiene a los vectores [math]\vec{j}[/math] y [math]\vec{k}[/math]) son nulas.

10 Tensión de Von Mises

Sea [math] \sigma_{VM}[/math] la Tensión de Von Mises una magnitud escalar usada en Cálculo de Estructuras como indicador para conocer el momento en el que un material termina un comportamiento elástico y comienza un comportamiento plástico. Dicha tensión se cálcula a partir de los valores de las tensiones principales de cada punto del espacio, las cuales corresponden con los autovalores [math]( \sigma_1,\sigma_2,\sigma_3)[/math] de [math] \sigma[/math]. La Tensión de Von Mises queda definida mediante la siguiente expresión:

[math] \sigma_{VM} = \sqrt{ \frac { ( \sigma_1 - \sigma_2 )^2 + ( \sigma_2 - \sigma_3 )^2 + ( \sigma_3 - \sigma_1 )^2 }{2} }[/math]

%discretizamos y creamos mallado
h=0.2;
x=-1:h:1;          
y=0:h:12;           
[X,Y]=meshgrid(x,y);  

%tensiones normales en i,j y k
i=(pi/9).*cos((pi/3).*Y);      
j=(pi/3).*cos((pi/3).*Y);      
k=(pi/9).*cos((pi/3).*Y);       
 
%Autovalores Von Mises
Z=ones(size(y,2),size(x,2));
T=zeros(3);
for m=1:length(y)
   for n=1:length(x)
      %Nombramos los valores de cada componente para la tensión Von Mises
      v1=i(m,n);
      v2=j(m,n);
      v3=k(m,n);
      A=[v1,0,0;0,v2,0;0,0,v3];
      %Obtenemos de los autovalores
      [V,D]=eig(A);
      A1=D(1,1);
      A2=D(2,2);
      A3=D(3,3);
      %Calculamos para cada componente
      TVM=sqrt(((A1-A2)^2+(A2-A3)^2+(A3-A1)^2)*1/2);
      Z(m,n)=TVM;
    end
end
 
%Representación en nuestra placa 2D
subplot(1,2,1)
surf(X,Y,Z)
shading flat
colorbar
axis equal
xlabel('Eje X')
ylabel('Eje Y')
title('Tensión de Von Mises')
view(2)
 
%Representación en 3D
subplot(1,2,2)
surf(X,Y,Z)
shading flat
colorbar
axis equal
xlabel('Eje X')
ylabel('Eje Y')
zlabel('Tensión de Mises')
title('Tensión de Von Mises representada en el eje Z')

Representación gráfica del campo escalar de Tensiones de Von Mises en la placa bidimensional

Queda, por tanto, claro que el mayor valor que alcanza la Tensión de Von Mises en la placa se alcanza periódicamente de forma sinusoidal.

11 Velocidad de propagación de las ondas en términos de las constantes de Lamé

En este caso consideramos que el tiempo [math]t[/math] en la ecuación del desplazamiento ondulatorio es distinto de cero [math]\vec{u} = \frac{1}{3}\cdot \sin( \frac{\pi}{3}y-v·t)\vec j[/math] por tanto, la Tensión [math]\sigma[/math] tendrá que volver a definirse, teniendo en cuenta esta vez, que el valor que adquiere [math]\epsilon[/math] aplicando el nuevo criterio de tiempo será el siguiente:

[math]\bigtriangledown\vec{u}=\ \begin{pmatrix}0&0&0\\0&\frac{\pi}{9}cos(\frac{\pi}{3}·y-v·t)&0\\0&0&0\end{pmatrix}[/math] , [math]\epsilon=\frac{\bigtriangledown\vec{u}\ +\bigtriangledown\vec{u}}{2}=\ \frac{\pi}{9}·cos( \ \frac{\pi}{3}y-v·t) \ \vec j \otimes \vec j[/math]

La Tensión será, por tanto:

[math]\sigma=λ\cdot\bigtriangledown \cdot \vec{u} \cdot I + 2µ\epsilon[/math][math]\ = λ \cdot \frac{\pi}{9}cos( \frac{\pi}{3}y-v·t) \cdot I + 2µ\ \frac{\pi}{9}·cos( \ \frac{\pi}{3}y-v·t) \ \vec j \otimes \vec j[/math]

[math]\sigma=\begin{pmatrix} \frac{\pi}{9}cos( \frac{\pi}{3}y-v·t)\cdotλ &0&0\\0& \frac{\pi}{9}cos( \frac{\pi}{3}y-v·t)\cdot( λ + 2µ)&0\\0&0& \frac{\pi}{9}cos( \frac{\pi}{3}y-v·t)\cdotλ \end{pmatrix}[/math]

Para obtener el valor de la velocidad de propagación, primero obtendremos los parámetros de la ecuación de la elasticidad lineal que rige el campo de fuerzas [math]\vec{F}[/math] que actúa sobre la placa y causa los desplazamientos observados. Dicha formula se define de la siguiente manera:

[math]\vec{F} = \frac{\partial^2\vec{u}}{\partial t^2} - \bigtriangledown \cdot \sigma[/math]

Se procede pues, a calcular ambos términos de la ecuación, [math] \frac{\partial^2\vec{u}}{\partial t^2}[/math] y [math]\bigtriangledown \cdot \sigma[/math], y posteriormente, igualar la ecuación a cero como ha sido indicado [math](\vec{F} \ =\ \vec{0})[/math].

[math]\bigtriangledown \cdot \sigma =\begin{pmatrix}0 \\ -\frac{\pi}{27}\cdot\sin( \frac{\pi}{3}y-v·t)\cdot( λ + 2µ) \\0 \\ \end{pmatrix}[/math], [math]\ \ \frac{\partial^2\vec{u}}{\partial t^2}= -\frac{v^2}{3}\cdot\sin( \frac{\pi}{3}y-v·t) \ \vec j[/math]

sustituyendo en la expresión se obtiene:

[math]\ \vec 0=-\frac{v^2}{3}\cdot\sin( \frac{\pi}{3}y-v·t) \ \vec j + \frac{\pi}{27}\cdot\sin( \frac{\pi}{3}y-v·t)\cdot( λ + 2µ) \ \vec j[/math]

Por último, operando la ecuación se obtiene fácilmente el valor (y dirección) de la velocidad de propagación.

[math] v = \sqrt{\frac{\pi}{9} \cdot ( λ + 2µ)} \ \vec j[/math]

11.1 Velocidad de propagación siendo la dirección de desplazamiento transversal

Si la onda de desplazamientos fuera transversal (tomando una amplitud [math]a= \frac{1}{3}\vec i [/math]) , la onda de desplazamientos adoptaría la siguiente expresión [math]\vec{u} = \frac{1}{3}\cdot \sin( \frac{\pi}{3}y-v·t)\vec i[/math]

la velocidad de propagación se calculará de la siguiente manera:

[math]\vec{F} = \frac{\partial^2\vec{u}}{\partial t^2} - \bigtriangledown \cdot \sigma[/math]

[math]\frac{\partial^2\vec{u}}{\partial t^2}= -\frac{v^2}{3}\cdot\sin( \frac{\pi}{3}y-v·t) \ \vec i[/math]

Mientras que el termino [math]\bigtriangledown \cdot \sigma[/math] sigue siendo el mismo,

[math]\bigtriangledown \cdot \sigma = -\frac{\pi}{27}\cdot\sin( \frac{\pi}{3}y-v·t)\cdot( λ + 2µ) \ \vec j[/math]

sustituyendo e igualando [math]\vec{F} =\ \vec 0[/math] ,

[math]\ \vec 0=-\frac{v^2}{3}\cdot\sin( \frac{\pi}{3}y-v·t) \ \vec i + \frac{\pi}{27}\cdot\sin( \frac{\pi}{3}y-v·t)\cdot( λ + 2µ) \ \vec j[/math]

Igualando ambas expresiones se llega a una igualdad vectorial que solo tendía sentido si el valor de la velocidad de propagación v fuera nulo [math] \ \vec v= \ \vec 0[/math], quedando demostrado por tanto que la velocidad de propagación [math] \ v[/math] difiere según el sentido de propagación de la onda.

12 Calculo del módulo de desplazamiento transversal

La fuerza que se aplica sobre la placa genera unos desplazamientos ondulatorios siguiendo la expresión de la forma [math]\vec{u} (x,y,t)=\frac{1}{3}sen(\frac{π}{3}-vt)\vec{j} [/math] correspondiendose asi con una onda longitudinal. Como se observa, la dirección de propagación de la onda corresponde con la dirección de la amplitud [math]\vec{a}=\frac{1}{3}\vec{j}[/math]. Por ello, los desplazamientos producidos en la placa serán únicamente en la dirección de propagación de la onda.

Para graficar el resultado mediante MATLAB, se ha empleado como velocidad de propagación la velocidad obtenida en el apartado 11 con los parámetros de Lamé correspondientes al apartado 8.

clear;
clc;
% Definir variables
%Definicion del punto de cálculo
x = 1/2; 
y=1;
%Velocidad de propagación 
v =sqrt(pi/3);   
%Intervalo de tiempo para analizar
t=0:0.1:10;   

% Calculo del módulo del vector
modulo_vector = abs((1/3) * sin((pi/3)*y - v*t));

% Graficación del resultado
axis equal
plot(t, modulo_vector, 'LineWidth', 2);
xlabel('Tiempo (s)');
ylabel('Módulo');
title('Módulo de desplazamiento');
grid on;

Representación gráfica de la función del modulo de desplazamientos en la dirección longitudinal