Visualización de campos escalares y vectoriales en elasticidad Grupo 41

Trabajo realizado por estudiantes
Título	Visualización de campos escalares y vectoriales en elasticidad Grupo 41
Asignatura	Teoría de Campos
Curso	2023-24
Autores	Sergio Menéndez Fernandez Alfonso José Fernández Sarmiento Jaime Rico Arroyo Mateo Navarro
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

En este trabajo consideraremos un placa rectangular situada en el espacio, representaremos la temperatura y la dirección de de los desplazamientos causados por la fuerza. Supondremos que la fuerza aplicada sobre la placa ha provocado un desplazamiento ondulatorio de los puntos de la misma dado por el vector U (expresado mas adelante). Gracias a todos los datos aportados nuestro objetivo será solucionar diferentes problemas que nos proponen y sacar alguna conclusión sobre el trabajo para mas tarde, en clase, poder exponerlo. Para ello utilizaremos diferentes mecanismos y formulas dados en clase.

1 Mallado de la placa

Lo primero que debemos hacer es representar el mallado de la superficie sobre la que trabajaremos, siguiendo asi las instrucciones indicadas de una toma de ejes (x,y)∈[-1,1] x [0,12] y un paso de muestreo h=2/10 para las variables x e y. Sobre esta superficie y una temperatura dada en el enunciado se basaran todas las operaciones (operadores diferenciales) y sus correspondientes representaciones de una forma gráfica realizadas con lenguaje m. Primero se realiza la definición de las variables, siendo estas vectores que comienzan en un valor y siguen una serie con el intervalo h hasta llegar al extremo del intervalo. Después se crea una malla formada por todos los puntos de x e y y posteriormente se calcula la z, siendo esta nula. Luego se representa la superficie y se observa como vista una representación en 2D y ya las ediciones y etiquetas a los distintos componentes de un gráfico.

Representación del mallado.

clear;clc;

%parametrización con muestreo h=2/10
x=linspace(-1,1,10);
y=linspace(0,12,60);
%mallado
[Mx,My]= meshgrid(x,y);             
Z=zeros(size(Mx));  
surf(Mx,My,Z)
%vista en 2D
view(2)
%nombramos los ejes
axis([-1.5,1.5,-0.5,12.5])
xlabel('Eje X')
ylabel('Eje Y')

2 Temperatura de la placa.

Este apartado consiste en un diagrama que muestra por colores la variación de temperatura a través de la placa y sus distintas curvas de nivel.

Representación del mallado.

clear;clc;
r=linspace(-1,1,70);
t=linspace(0,12,45);
[x,y]= meshgrid(r,t);      
          
T=3.*log(1.+(x+1).^2)+log(1.+(y-2).^2);

hold on
subplot(1,2,1)                      %División de la ventana en dos                        
surf(x,y,T)
view(2)
axis([-1.5,1.5,-0.5,12.5])
xlabel('Eje X')
ylabel('Eje Y')
colorbar                                    %Mostramos la escala
subplot(1,2,2)                              %Trabajamos en la segunda ventana
contour(x,y,T,30)
title('Curvas de Nivel','Fontsize',16)      %Líneas de nivel
colorbar                                    %Mostramos las escala
axis([-1.5,1.5,-0.5,12.5])
xlabel('Eje X')
ylabel('Eje Y')
hold off

En la grafica de la izquierda podemos observar que a medida que nos desplazamos hacia el suroeste la temperatura disminuye, por lo que la zona mas cálida se encuentra en la zona de x máxima e y máxima.

3 Gradiente de la temperatura y su representación.

El siguiente código muestra como hemos calculado el gradiente del campo escalar temperatura calculado en el apartado anterior. A la derecha vemos representado dicho campo en sus curvas de nivel con el gradiente siendo perfectamente ortogonal a las curvas isotermas. Dicho gradiente previamente calculado, se representa con el comando quiver, el cual crea las flechas representadas en la gráfica.

Representación del mallado.

clear;clc;
%---Divergencia del campo de temperaturas---
r=linspace(-1,1,10);
t=linspace(0,12,45);
[x,y]= meshgrid(r,t);                           %hasta aquí es para incializar lel programa.

u=(-1).*(3.*((2.*x)+2)./(1+(x+1).^2));              %componente i
v=(-1).*(((2.*y)-4)./(1+(y-2).^2));                   %componente j

T=3.*log(1+(x+1).^2)+log(1+(y-2).^2);

hold on                                     %Campo escalar de temperatura            
contour(x,y,T,30)                           %Líneas de nivel
colorbar                                    %Mostramos las escala
w=quiver(x,y,u,v,0.40,'k');
title ('Gradiente de la Temperatura','Fontsize',18);
set(w,'maxheadsize')
axis([-1.5,1.5,-0.5,12.5])
xlabel('Eje X')
ylabel('Eje Y')
hold off

4 Campo de deformaciones para el instante inicial.

El campo de vectores viene dado por [math]\vec{u}(x, y, t)=\vec{a}sin(k\pi(\vec{d}·\vec{r_{0}}(x,y)-vt)).[/math] asi mismo la representación de este para t=0 se hará tomando como valores para la amplitud, [math]\vec{a}= 1/3 \vec{j}[/math] ; el número de onda, k=1 ; y el vector unitario de propagación, [math]\vec{d}= 1/3 \vec{j}[/math]. Por lo que el campo de vectores vendrá dado por la ecuación:

[math]\vec{u}(x, y)=1/3 \vec{j} sin (\pi( 1/3 \vec{j}·\vec{r_{0}}(x,y))).[/math]

Siendo [math]\vec{r_{0}}(x, y)= x \vec{i} + y \vec{j} [/math]

Por lo que operando nos quedaría:

[math]\vec{u}(x, y)=\frac{1}{3}sin ((pi*y)/3) \vec j\ [/math]

Representación del mallado.

clear;clc
  x=-1:0.2:1;           
  y=0:0.2:12;
  [MX,MY]=meshgrid(x,y)              
  Ux=(0.*X);  % Componentes en la dirección i 
  Uy=((1/3).*sin((1/3).*pi.*Y));                       % Componentes en la dirección j 
  figure
  mesh(MX,MY,Y*0)                    % Creación del mallado
 axis equal                  
 hold on
 quiver(MX,MY,Ux,Uy,"r")            % Campo de vectores, "r" define el color en la grafica 
 title("Campo de desplazamientos")  
 xlabel("Eje X")
 ylabel("Eje Y")
 view(2)
 hold off

5 Placa antes y después del desplazamiento

En la placa se produce una deformación provocada por el campo vectorial, por lo que estudiamos el sólido antes y después de su deformación en t = 0, conociendo el vector de posición de los puntos de la placa antes [math]\vec{r_{0}}(x, y)= x \vec{i} + y \vec{j} [/math] y la posición de cada punto (x,y) despues. Para poder compararlos usaremos el comando subplot y visualizaremos asi la la gráfica antes y después.

Representación del mallado.

clear;clc
  x=-1:0.2:1;           
  y=0:0.2:12;           
  [X,Y]=meshgrid(x,y);  
  % Definimos la posición final 
  rx=(0.*X)+X;
  ry=((1/3).*sin((1/3).*pi.*Y))+Y;
  % Antes de desplazamiento
 subplot(1,2,1)
 surf(X,Y,0*Y)
 view(2)
 title("Antes del desplazamiento") 
 xlabel("Eje X")
 ylabel("Eje Y")
 % Después del desplazamiento
 subplot(1,2,2)
  surf(rx,ry,ry*0)
 view(2)
 title("Después del desplazamiento") 
 xlabel("Eje X")
ylabel("Eje y")

6 Divergencia del campo vectorial sobre la capa y su representación grafica

La divergencia nos proporciona la diferencia entre el flujo saliente y el flujo entrante de un campo vectorial sobre la superficie que rodea al volumen, por lo que si el campo dispone de “fuentes” esta será positiva y al contrario, si dispone de “sumideros” será negativa.

Procederemos a calcular así la divergencia del campo:

[math]\vec{u}(x, y)=\frac{1}{3}sin ((π*y)/3) \vec j\ [/math]

Entonces:

[math]\ \nabla \cdot \vec u =\frac{\partial \vec u_1}{\partial x} + \frac{\partial \vec u_2}{\partial y} + \frac{\partial \vec u_3}{\partial z} = \frac{π}{9} \cos((π * y)/3) [/math]

Figura 2. curvas de nivel del campo de temperaturas

Figura 2. curvas de nivel del campo de temperaturas

clear;clc
 h=2/10; 
  x=-1:h:1;
  y=0:h:12;
  [xx,yy]=meshgrid(x,y);
  diver=(pi/9).*cos((pi/3).*yy);
  %Representación
  surf(xx,yy,diver)
  shading interp
  colorbar
 title('Divergencia')
 xlabel('x')
 ylabel('y')
 zlabel('z')
 %Dibujo de la divergencia
 figure
 pcolor(xx,yy,diver)
 shading interp
 hold on
 contour(xx,yy,diver,'k')
 title('Divergencia')
 xlabel('x')
 ylabel('y')
 colorbar
 axis([-1.5,1.5,-0.5,12.5])
view(2)
 hold off
 maximo=max(max(diver))
 minimo=min(min(diver))

7 Calcular |∇ × [math]\vec{u}[/math] | en todos los puntos del solido en t = 0

En el instante inicial obtendremos el siguiente campo: [math]\vec{u}({x},{y},{0}) = (\vec{a}*sin(π*k(\vec{d}*(x\vec{i}+y\vec{j}))) [/math]

Tomaremos en particular : [math]\vec{a} =\vec{d} = (1/3)\vec{j}[/math], k=1.

Obteniendo el siguiente campo : [math]\vec{u}(x,y) = 1/3sin((π/3)*y)\vec{j}[/math]

Para calcular el rotacional de un campo de desplazamiento ondulatorio, se aplica la fórmula del producto vectorial en los ejes del sistema de coordenadas.

[math]\ \bigtriangledown \times \vec{u} = \begin{vmatrix}\vec{i} & \vec{j} &\vec{k} \\ d/dx & d/dy & d/dt\\ x & y & t\end{vmatrix}= \begin{vmatrix}\vec{i} & \vec{j} & \vec{k} \\ d/dx & d/dy & d/dt \\ 0 & 1/3sin((π/3)×y) & 0\end{vmatrix}=(0-0)\vec{i}+(0-0)\vec{j}+(0-0)\vec{k}=0 [/math]

Por lo tanto, se trata de un campo conservativo es decir, el rotacional en todos los puntos de la placa es nulo, lo que implica: [math]\ \left | \bigtriangledown \times \vec{u} \right | = 0 [/math]

El rotacional indica la dirección y velocidad de giro del campo [math]\vec{u} [/math] en cada punto, si nos fijamos en la representación gráfica, se puede apreciar que no existe ninguna rotación.

Representación del mallado.

8 Tensor de deformaciones

Siendo [math]ε(\vec u)=\frac {\nabla \vec u + \nabla \vec u ^t}{2}[/math], la parte simétrica del tensor gradiente de [math]\vec u[/math] conocido como tensor de deformaciones. En un medio elástico lineal, isótropo y homogéneo los desplazamientos permiten escribir el tensor de tensiones a través de la fórmula [math]σ=λ\nabla·\vec u 1+2με[/math].

[math]ε(\vec u)=\nabla \vec u= \begin{pmatrix} 0 & 0 & 0\\ 0 & \frac{π}{9}cos(\frac{π}{3}y) & 0 \\ 0 & 0 & 0 \end{pmatrix}[/math]

[math]σ=λ\nabla·\vec u 1+2με = λ\frac{π}{9}cos(\frac{π}{3}y)·\begin{pmatrix} 1 & 0 & 0\\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} + 2μ\begin{pmatrix}0 & 0 & 0 \\ 0 & \frac{π}{9}cos(\frac{π}{3}y) & 0 \\ 0 & 0 & 0 \end{pmatrix} = 2λ\frac{cosx}{5}·I+\frac{π}{9}cos(\frac{π}{3}y)\begin{pmatrix} 0 & 0 & 0\\ 0 & 2μ & 0 \\ 0 & 0 & 0\end{pmatrix} = [/math]

[math]= \frac{π}{9}cos(\frac{π}{3}y) (\begin{pmatrix} λ & 0 & 0\\ 0 & λ & 0 \\ 0 & 0 & λ \end{pmatrix} + \begin{pmatrix} 0 & 0 & 0\\ 0 & 2μ & 0 \\ 0 & 0 & 0 \end{pmatrix}) = \frac{π}{9}cos(\frac{π}{3}y)\begin{pmatrix} λ & 0 & 0\\ 0 & 2μ+λ & 0 \\ 0 & 0 & λ \end{pmatrix}→{λ=μ=1}→\frac{π}{9}cos(\frac{π}{3}y)\begin{pmatrix} 1 & 0 & 0\\ 0 & 3 & 0 \\ 0 & 0 & 1 \end{pmatrix}[/math]

A continuación, calculamos las tensiones normales en la dirección que marca el eje [math]\vec i[/math], es decir [math]\vec i · σ ·\vec i[/math], las tensiones normales en la dirección que marca el eje
[math]\vec j[/math], es decir [math]\vec j · σ · \vec j [/math] y las correspondientes al eje [math]\vec k[/math], es decir [math]\vec k · σ ·\vec k[/math].

• [math]\vec i · σ ·\vec i = \begin{pmatrix} 1 & 0 & 0 \end{pmatrix}·\begin{pmatrix} \frac{π}{9}cos(\frac{π}{3}y) & 0 & 0\\ 0 & \frac{π}{3}cos(\frac{π}{3}y) & 0 \\ 0 & 0 & \frac{π}{9}cos(\frac{π}{3}y) \end{pmatrix}·\begin{pmatrix} 1 \\0 \\ 0 \end{pmatrix} = \frac{π}{9}cos(\frac{π}{3}y) [/math]

• [math]\vec j · σ ·\vec j = \begin{pmatrix} 0 & 1 & 0 \end{pmatrix}·\begin{pmatrix} \frac{π}{9}cos(\frac{π}{3}y) & 0 & 0\\ 0 & \frac{π}{3}cos(\frac{π}{3}y) & 0 \\ 0 & 0 & \frac{π}{9}cos(\frac{π}{3}y) \end{pmatrix}·\begin{pmatrix} 0 \\1 \\ 0 \end{pmatrix} = \frac{π}{3}cos(\frac{π}{3}y) [/math]

• [math]\vec k · σ ·\vec k = \begin{pmatrix} 0 & 0 & 1 \end{pmatrix}·\begin{pmatrix}\frac{π}{9}cos(\frac{π}{3}y) & 0 & 0\\ 0 & \frac{π}{3}cos(\frac{π}{3}y) & 0 \\ 0 & 0 & \frac{π}{9}cos(\frac{π}{3}y) \end{pmatrix}·\begin{pmatrix} 0 \\0 \\ 1 \end{pmatrix} = \frac{π}{9}cos(\frac{π}{3}y)[/math]

Representación del mallado.

%____Tensor de deformaciones____
x=linspace(-1,1,10);
y=linspace(0,12,45);    
[Mx,My]= meshgrid(x,y); 
Ti=(pi/9).*cos((pi/3)).*My;       % Tension normal en i
Tj=(pi/3).*cos((pi/3)).*My;      % Tension normal en j
Tk=(pi/9).*cos((pi/3)).*My;        % Tension normal en k
figure
subplot(1,3,1)           
pcolor(Mx,My,Ti)
colorbar
shading flat
axis([-1.5,1.5,-0.5,12.5])
xlabel('Eje X')
ylabel('Eje Y') 
zlabel('Eje Z')

title('Tensiones normales en i')  
subplot(1,3,2) 
pcolor(Mx,My,Tj)
colorbar
shading flat
axis([-1.5,1.5,-0.5,12.5])
xlabel('Eje X')
ylabel('Eje Y') 
zlabel('Eje Z')

title('Tensiones normales en j')
subplot(1,3,3)
pcolor(Mx,My,Tk)
colorbar
shading flat
axis([-1.5,1.5,-0.5,12.5])
xlabel('Eje X')
ylabel('Eje Y') 
zlabel('Eje Z')
title('Tensiones normales en k')

9 Tensiones tangenciales

La tensión tangencial respecto al plano ortogonal a [math]\vec{i}[/math] (en t=0), viene dada por la siguiente expresión:

[math]\ \left| \sigma ·\vec i -(\vec i·\sigma·\vec i)\vec i \right| [/math]

[math] \sigma ·\vec i = \begin{pmatrix} \frac{π}{9}cos(\frac{π}{3}y) & 0 & 0\\ 0 & \frac{π}{3}cos(\frac{π}{3}y) & 0 \\ 0 & 0 & \frac{π}{9}cos(\frac{π}{3}y) \end{pmatrix}·\begin{pmatrix} 1 \\0 \\ 0 \end{pmatrix} = \frac{π}{9}cos(\frac{π}{3}y)\vec i [/math]

[math](\vec i · σ ·\vec i)\vec i = \begin{pmatrix} 1 & 0 & 0 \end{pmatrix}·\begin{pmatrix} \frac{π}{9}cos(\frac{π}{3}y) & 0 & 0\\ 0 & \frac{π}{3}cos(\frac{π}{3}y) & 0 \\ 0 & 0 & \frac{π}{9}cos(\frac{π}{3}y) \end{pmatrix}·\begin{pmatrix} 1 \\0 \\ 0 \end{pmatrix} = \frac{π}{9}cos(\frac{π}{3}y)\vec i [/math]

Por lo tanto:

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

Lo que quiere decir que la tensión tangencial respecto al plano que contiene a los vectores [math]\vec{j}[/math] y [math]\vec{k}[/math] es nula.

10 Tensión de Von Mises

La tensión de Von Mises es un concepto crucial en el campo de la ingeniería y la mecánica de materiales. Esta medida se utiliza para evaluar el estado de esfuerzos en un material, especialmente cuando este se somete a cargas externas. Cuando un material se encuentra bajo la influencia de fuerzas, ya sean tensiones o compresiones, los diferentes puntos en el interior del material experimentan distintos niveles de esfuerzos.

La fórmula de la tensión de Von Mises combina las tensiones principales, que son los esfuerzos máximos en direcciones específicas, para proporcionar un único valor que representa el esfuerzo equivalente. Esto es esencial porque permite simplificar y analizar el comportamiento del material de manera más eficiente.

La tensión de Von Mises en nuestro ejercicio puede obtenerse gracias a las tensiones principales del tensor tensión en un punto de un sólido deformable, con la ayuda de la siguiente formula:

[math]σ=\sqrt{\frac {(σ_1-σ_2)^2+(σ_2-σ_3)^2+(σ_3-σ_1)^2}{2}}[/math]

Donde [math]σ_1,σ_2,σ_3[/math] son los autovalores de [math]σ=\frac{π*cos(\frac{π*y}{3})}{9}\begin{pmatrix} 1 & 0 & 0\\ 0 & 3 & 0 \\ 0 & 0 & 1 \end{pmatrix}[/math].

Gracias a esto podemos obtener que en nuestro ejercicio la tension seria: 0.8533

Representación del mallado.

x=linspace(-1,1,50);                   % Vector x con valores entre 0 y 10    
y=linspace(0,12,50);                     % Vector x con valores entre 0 y 2      
[X,Y]=meshgrid(x,y);                     % Obtención del mallado      
% Definimos nuestra función                                               
fx=inline('(pi/9)*cos(pi*y/3)','x','y');        
fy=inline('pi/3)*cos(pi*y/3)','x','y');
fz=inline('(pi/9)*cos(pi*y/3)','x','y');
% Inicio del bucle
for i= 1:length(x)                              
   for j= 1:length(y)                             
% Asignación de valores para cada componente del mallado                                                                    
xx=fx(X(i,j),Y(i,j));                    % Respecto a fx
yy=fy(X(i,j),Y(i,j));                    % Respecto a fy
zz=fz(X(i,j),Y(i,j));                    % Respecto a fz
% Creación de un vector con las componentes [xx, yy, zz]
v=[xx yy zz];                            
M=diag(v);                               % Diagonalización del vector       
autov=eig(M);                            % Obtención de los autovalores      

% Cálculo de la tensión de Von Mises para cada componente                                                
VonM=sqrt(((autov(1)-autov(2))^2+((autov(2)-autov(3))^2+((autov(3)-autov(1))^2))*1/2));  
Z(i,j)=VonM;
    end
end
% Fin del bucle                                             

subplot(2,1,1)
mesh(X,Y,Z)                                     % Representación superficial con líneas entrecruzadas
title('TENSIÓN DE VON MISES')                   % Título del gráfico
xlabel('Eje X')                                 % Título del eje x
ylabel('Eje Y')                                 % Título del eje y
zlabel('Eje Z')                                 % Título del eje z

subplot(2,1,2)
plot(Y,Z,'*r')                                  % Representación superficial con líneas entrecruzadas
title('TENSIÓN DE VON MISES EN 2D')             % Título del gráfico
xlabel('Eje X')                                 % Título del eje x
ylabel('Eje Z')                                 % Título del eje y

maximo=max(max(Z));                           % Obtención del valor máximo

11 Velocidad de propagación v en términos de Lamé

La velocidad de propagación v en un material elástico se puede expresar en términos de los parámetros de Lamé, que son dos constantes elásticas λ (lambda) y μ (mu). Estos parámetros son fundamentales para caracterizar las propiedades elásticas de un material.

Para poder calcular la velocidad de propagación de las ondas v en términos de las constantes de Lamé necesitaremos suponer que F=0.

El desplazamiento en la placa viene dado por las fuerzas [math]\vec{F}[/math] que actúan sobre dicha placa, las cuales se pueden aproximar con la siguiente función de elasticidad:

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

en nuestro caso el vector u seria:

[math]\vec{u}[/math] = [math]\frac{1}{3}\sin(\frac{πy}{3}-vt)\vec{j}[/math]

siendo 'v' la velocidad de propagación y apoyandonos en los calculos que hemos realizado en el apartado 8, [math]\sigma[/math] será :

[math]\sigma = \begin{pmatrix} \frac{π}{9}cos(\frac{πy}{3}-vt) & 0 & 0 \\ 0 &\frac{π}{3}cos(\frac{πy}{3}-vt) & 0 \\ 0 & 0 & \frac{π}{9}cos(\frac{πy}{3}-vt)\end{pmatrix}[/math]

Gracias a esto sabemos que la divergencia [math]\nabla \cdot \sigma[/math] sería de la siguiente manera:

[math]\nabla \cdot \sigma[/math]=[math]\frac{-\pi^2}{9}sen(\frac{\pi y}{3}-vt)[/math]

Por otra parte tendremos que calcular

[math]\frac{d^2\vec{u}}{dt^2} [/math]

Haciendo todas las operaciones hallaremos

[math]\frac{d^2\vec{u}}{dt^2}= \frac{-v^2}{3}sen(\frac{\pi y}{3}-vt)\vec{j} [/math]

Ahora, sabiendo que [math]\vec{F}[/math] = 0 igualaremos [math]\frac{d^2\vec{u}}{dt^2} = \nabla \cdot \sigma[/math] Y con esto hallaremos finalmente que

[math] \vec{v}= \frac{\pi}{\sqrt{3}}\vec{j}[/math]

12 Módulo de desplazamiento en dirección [math]\vec{i}[/math] y representacion

La onda es longitudinal, es decir, que se desplaza en el mismo sentido que la amplitud, por tanto no habrá desplazamiento en otro eje distinto al eje en el que se transmite dicha onda, de aquí la demostración:

Adoptando el campo de desplazamientos: [math]\vec u (x,y,t)=\vec a\ sin(π k (\vec d\ \vec r\ (x,y) - v t))[/math]

Fijando el punto [math](x,y)=(1/2,1)[/math] se tiene: [math]\vec u (1/2,1,t)=\frac{1}{3}\ sen(π/3-vπt) \vec j\ [/math]donde el parámetro [math]t[/math], tiempo, pertence al intervalo [math][0,10][/math]

En la dirección de [math]\vec i [/math] se obtiene que: [math]\vec u \cdot \vec i = \frac{1}{3}\ sen(π/3-vπt) \vec j\ \vec i\ =\vec 0\ [/math]

Es decir, al ser una onda longitudinal en la que la dirección de propagación es la misma que la amplitud, no se producen desplazamientos en el eje [math]\vec{i}[/math] , demostrando así lo dicho arriba.

clear;clc;
% Definición de parámetros
My = 1; % Puedes ajustar este valor según sea necesario, sirve para ajustar la frecuencia de la funcion seno, es decir, las oscilaciones.
v=1;    %Lo defino como un vector unitario de modulo 1 que insertare en la direccion j.

% Creación de un mallado rectangular
x = linspace(-1, 1, 30); % valores x
t = linspace(0, 10, 30); % valores de tiempo t
[X, T] = meshgrid(x, t); % creación de la malla

% Cálculo de la función de desplazamiento
j = (1/3) * sin((pi/3) * My - pi * T);
v = cos((pi/3) .* My - v*pi * T); % vector unitario en dirección j

% Visualización del resultado
figure;
quiver(X, T, zeros(size(v)), v, 'AutoScale', 'off', 'Color', 'r'); % quiver plot
title('Desplazamiento en dirección j');
xlabel('Eje x');
ylabel('t');
axis([-1.5 1.5 -0.5 10.5]); % ajustar los ejes según sea necesario
grid on;

Representación del mallado.

Tras mucho tiempo dedicado a este apartado, sencillo en los números pero más complejo en la programación; creemos que puede haber un fallo al no salir como creemos que tiene que salir. Hemos estado investigando y concretando y no hemos llegado a la solución deseada.