Visualización de campos escalares y vectoriales en elasticidad (Grupo 16-A)

El presente trabajo tiene por objeto la visualización y representación de los campos escalares y vectoriales sobre un sólido que experimenta una deformación. Se ha considerado como objeto de estudio un sólido rectangular (en 2D) cuya deformación viene dada por la temperatura y los desplazamientos producidos por la acción de una fuerza determinada.

Para su realización se empleará el software de programación y cálculo numérico Matlab.

Trabajo realizado por estudiantes
Título	Visualización de campos escalares y vectoriales en elasticidad (Grupo 16-A)
Asignatura	Teoría de Campos
Curso	2023-24
Autores	Lara Gutiérrez Kreutzer Miguel Moreno Martín Luis Anthony Vera Guapi Emilio Villegas Maroto
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Representación de la placa

En primer lugar, se representa la placa dibujando los puntos interiores del sólido mediante un mallado. Este está definido en la región [math](x, y) ∈ [-1, 1]×[0, 12][/math] y como paso de muestreo se empleará h=2/10 para ambas variables.

Para su realización, se ha comenzado definiendo los vectores x e y con los valores que toma el rectángulo. Dichos vectores se han metido en dos matrices, X e Y, empleando el comando meshgrid() y con el comando mesh() se crea el mallado buscado. Este comando requiere tres matrices, por lo que en la ultima entrada multiplicamos cualquiera de las obtenidas por cero, puesto que se trata de una placa plana. Por último, con el comando axis() se definen los valores que tomarán los ejes de la gráfica [math](x, y) ∈ [-2, 2]×[0, 13][/math] y nombramos tanto la gráfica como los ejes.

A continuación, se muestra el código empleado para la obtención de la placa y la imagen resultante:

Figura 1. Representación del mallado

x=-1:0.2:1;           % Vector x con valores entre -1 y 1 con paso de muestreo h=2/10
y=0:0.2:12;           % Vector y con valores entre 0 y 12 con paso de muestreo h=2/10
[X,Y]=meshgrid(x,y);  % Creación del mallado con los vectores x e y 
mesh(X,Y,0*X);        % Dibujo de la malla pedida
axis equal
axis([-2,2,0,13])     % Con el comando axis se definen los valores mínimos y máximos de los ejes de la gráfica
title("Mallado de la placa rectangular plana")
xlabel("Eje X")       % Título del eje X
ylabel("Eje Y")       % Título del eje Y
view(2)               % Visualiamos el gráfico en dos dimensiones

2 Temperatura de la placa

La temperatura del sólido viene dada por la función: [math]T(x, y) = 3log(1 +(x-1)^2)+log(1+(y-8)^2)[/math]

Primero representaremos las distintas curvas de nivel del campo de temperaturas mediante el siguiente código:

Figura 2. curvas de nivel del campo de temperaturas

h=2/10; %Paso de muestreo
x=-1:h:1;
y=0:h:12; 
[X,Y]=meshgrid(x,y)
T=3*log(1 +(X-1).^2)+log(1+(Y-8).^2); %Función que representa la temperatura
contour(x,y,T,20); % Se dibujan las curvas de nivel 
axis equal;
axis([-1,1,0,12]);
view(2); %Vista en planta
title("Curvas de nivel del campo de temperaturas") %Nombre de la figura
xlabel("Eje X")       % Título del eje X
ylabel("Eje Y")       % Título del eje Y
colorbar

A continuación, mostraremos el campo de temperaturas, que representa la variación de temperatura a lo largo de la placa. Usaremos el siguiente código:

Figura 3. Campo de temperaturas

h=2/10; %Paso de muestreo
 r=-1:h:1; % Discretización del radio
t=0:h:12; %Discretización del ángulo
[X,Y]=meshgrid(r,t); %Generación de la retícula rectangular
T=3*log(1 +(X-1).^2)+log(1+(Y-8).^2); %Función que representa la temperatura
surf(X,Y,T); %Se dibuja la placa
shading interp %Se realiza un degradado mediante interpolación con los colores de relleno
axis equal
axis([-1,1,0,12]);
colorbar
view(2)
title("Campo de temperaturas") %Nombre de la figura
Tmax=max(max(T)); %Se calcula la temperatura máxima alcanzada
fprintf('La temperatura máxima alcanzada en la placa es de: %1.4f ºC\n' , Tmax)

Como resultado obtenemos que la temperatura máxima alcanzada en la placa es de: 9.0027 ºC.

INTERPRETACIÓN DEL CAMPO DE TEMPERATURAS

En la imagen podemos apreciar dos gamas de colores distintos, una primera de colores fríos (azul) que corresponden a unas temperaturas más bajas, mientras que la otra se trata de una gama de colores más cálidos (amarillo) que corresponden a temperaturas mayores.

2.1 Gradiente de la temperatura

En segundo lugar, calcularemos el gradiente de la temperatura, que es un vector cuya primera componente es la derivada del campo escalar de la temperatura respecto de [math]x[/math], multiplicada por el vector [math]\vec{i}[/math] y como segunda componente la derivada del campo escalar respecto de [math]y[/math], multiplicada por el vector [math]\vec{j}[/math].

[math]∇T = \frac{\partial T}{\partial x}\vec{i} + \frac{\partial T}{\partial y}\vec{j} = \frac{6(x-1)}{1+(x-1)^2}\vec{i} + \frac{2(y-8)}{1+(y-8)^2}\vec{j} [/math]

Una vez calculado el gradiente, crearemos las variables [math]dx[/math] y [math]dy[/math] en función de las matrices [math]X[/math] e [math]Y[/math]. Por último, utilizaremos el comando [math]quiver ()[/math] para representar el gradiente de la temperatura.

Al tener que comprobar que las curvas de nivel del campo de temperatura y el gradiente son perprendiculares, nos ayudaremos del comando [math]hold on[/math] para poder representar las curvas de nivel y el gradiente en el mismo gráfico.

Figura 4.Gradiente de Temperatura

h=2/10; %Paso de muestreo
 x=-1:h:1; %Discretización del radio
 y=0:h:12; %Discretización del ángulo
 [X,Y]=meshgrid(x,y); %Generación de la retícula rectangular
 T=3*log(1 +(X-1).^2)+log(1+(Y-8).^2); %Función de temperatura
 dx=(6.*(X-1))./(1+(X-1).^2);
 dy=(2.*(Y-8))./(1+(Y-8).^2);
 axis equal;
 hold on
 quiver(x,y,dx,dy); %Se dibuja el gradiente
 contour(x,y,T,20); %Se dibujan las curvas de nivel de la temperatura
 axis([-1,1,0,12]);
 title("Gradiente de temperatura");
 hold off

Como no se puede apreciar bien el gradiente en la imagen general, hemos realizado la siguiente ampliación donde se pude observar la ortogonalidad del gradiente.

INTERPRETACIÓN DEL GRADIENTE DEL CAMPO DE TEMPERATURAS

El gradiente apunta a los puntos de la gráfica a los cuales la gráfica tiene un mayor incremento de temperatura.

3 Ley de Fourier

De acuerdo a la Ley de Fourier la energía calorífica Q viaja de acuerdo a la fórmula Q = −κ∇T, donde κ es la constante de conductividad térmica de la placa que supondremos κ = 1.

[math]Q=-∇T = -(\frac{\partial T}{\partial x}\vec{i} + \frac{\partial T}{\partial y}\vec{j}) = -\frac{6(x-1)}{1+(x-1)^2}\vec{i} - \frac{2(y-8)}{1+(y-8)^2}\vec{j} [/math]

Figura 5. Ley de Fourier

h=2/10; %Paso de muestreo
 x=-1:h:1; %Discretización del radio
 y=0:h:12; %Discretización del ángulo
 [X,Y]=meshgrid(x,y); %Generación de la retícula rectangular
 T=3*log(1 +(X-1).^2)+log(1+(Y-8).^2); %Función de temperatura
 k=1; %Constante de conductividad térmica
 dx=-(6.*(X-1))./(1+(X-1).^2);
 dy=-(2.*(Y-8))./(1+(Y-8).^2);
 axis equal;
 hold on
 quiver(x,y,dx,dy); %Se dibuja el gradiente
 contour(x,y,T,20); %Se dibujan las curvas de nivel de la temperatura
 axis([-1,1,0,12]);
 title("Ley de Fourier");
 hold off

INTERPRETACIÓN DE LA ENERGÍA CALORÍFICA

Podemos apreciar la similitud con el gradiente. Según la Ley de Fourier el flujo de transferencia de calor es proporcional y de sentido contrario al gradiente de temperatura en esa dirección, tal y como podemos observar en la imagen.

4 Campo de vectores

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

[math]\vec{u}(x, y)=x/3 \vec{i} sin (\pi( 1/12 \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:

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

Primero creamos las matrices X e Y de los puntos del mallado de la placa. Definimos las componentes en la dirección de [math]\vec{i}[/math] y [math]\vec{j}[/math] del campo de desplazamiento y con el comando mesh() dibujamos el mallado sobre el que vamos a poder ver el campo de vectores creados con quiver (). Los comandos hold on y hold off permiten ver los vectores sobre el sólido sin que este último desaparezca.

Figura 6. Campo de vectores en t=0.

clear
close all
x=-1:0.2:1;           
y=0:0.2:12;
[X,Y]=meshgrid(x,y)              
ux=(X.*(1/3).*sin((pi/12).*Y));  % Componentes en la dirección i 
uy=(0.*Y);                       % Componentes en la dirección j 
figure
mesh(X,Y,X*0)                    % Creación del mallado
axis equal                  
hold on
quiver(X,Y,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

INTERPRETACIÓN DEL CAMPO DE VECTORES

El campo de vectores queda representado en color rojo sobre la malla del sólido. Puesto que el campo de vectores se intensifica lateralmente siendo máximo en los cantos, podemos esperar que la deformación provoque un bombeo/ curvatura del sólido.

5 Comparación de la placa antes y después del desplazamiento

La placa sufre una deformación provocada por el campo de vectores. Estudiamos el sólido antes y después de su deformación para 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) después. Para que resulte más sencilla su interpretación, se visualizará la gráfica antes, después y la comparación con ayuda del comando subplot.

Comenzamos representando la placa antes de la deformación como se ha hecho anteriormente, definimos los valores que toma el rectángulo y lo visualizamos con el comando surf(). A continuación, definimos los puntos después del desplazamiento y los representamos de la misma manera. Para la comparación, se crea un gráfico que muestre mediante líneas el sólido sin deformaciones y sobre este, con el comando hold on y hold of, que se muestre con líneas el sólido deformado.

Figura 7. Sólido antes de la deformación y después.

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

6 Divergencia

La divergencia mide la diferencia entre el flujo saliente y el flujo entrante de un campo vectorial sobre la superficie que rodea a un volumen de control, por tanto, si el campo tiene "fuentes" la divergencia sera positiva, y si tiene "sumideros" sera negativa.

Calculamos la divergencia de:

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

Entonces:

[math]\nabla ·\vec u = \frac{\partial u_1}{\partial x} + \frac{\partial u_2}{\partial y} = \frac{1}{3}sin (\frac{\pi*y}{12}) [/math]

Utilizaremos el codigo que se muestra a continuación para obtener los valores máximos y mínimos de la divergencia.

Divergencia 2D

Divergencia 3D

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

En la grafica se puede apreciar que la divergencia máxima es de 0.333 que se produce en los puntos donde existe mayor desplazamiento y la mínima es también la nula, es decir, 0, que coincide con los puntos donde no existe desplazamiento.

7 Rotacional del sólido

En este siguiente apartado calcularemos el rotacional y a continuación lo representaremos.

[math]∇ × \vec{u}= \begin{vmatrix}\vec{i} & \vec{j}& \vec{k}\\ \frac{\partial }{\partial x} & \frac{\partial }{\partial y}& \frac{\partial }{\partial z}\\ u_{1}& u_{2}& u_{3} \end{vmatrix} = \begin{vmatrix} \vec{i} & \vec{j} & \vec{k}\\ \frac{\partial }{\partial x} & \frac{\partial }{\partial y} &\frac{\partial }{\partial z} \\ \frac{x}{3}\sin(\frac{pi*y}{12}) & 0 & 0\end{vmatrix} = -\frac{\pi*x}{36}*cos(\frac{pi*y}{12})\vec{k} [/math].

Módulo del rotacional

clc
clear
close all

%Definición de regiones
h=2/10; 
x=[-1:h:1];          
y=[0:h:12];
%Matriz de x e y
[Mx,My]=meshgrid(x,y);
%Cálculo del módulo del rotacional
rot=-((Mx.*pi)/36).*cos((pi.*My)/12)
%Representación gráfica
surf(Mx,My,rot)
shading flat
axis equal
colorbar
view(2);
axis([-1,1,0,12]);
title('Módulo del rotacional');
xlabel('Eje X');
ylabel('Eje Y')

INTERPRETACIÓN DEL ROTACIONAL

Como podemos observar los puntos que sufren un mayor rotacional son aquellos representados en color amarillo, estos tendrán una mayor capacidad de giro ya que el rotacional nos indica la capacidad que tienen los distintos puntos de nuestra placa para girar sobre otro punto determinado.

8 Tensiones en las direcciones [math]\vec i [/math], [math]\vec j [/math], [math]\vec k [/math]

Definamos [math]e(\vec{u}) = (∇\vec{u} + \triangledown \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 [math]σ_{ij}[/math] a través de la fórmula [math]σ = λ∇ · \vec{u} 1 + 2 μe[/math], donde λ y μ son los conocidos como coeficientes de Lamé que dependen de las propiedades elásticas de cada material. A pesar de que los desplazamientos son planos (es decir [math]\vec{u}[/math] no tiene componente en la dirección de [math]\vec{k}[/math]) las tensiones no tienen por qué ser planas y puede haber tensiones en la dirección ortogonal al plano de la placa. Tomando [math]λ = μ = 1[/math]. Asi que pasamos a dibujar 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].

Primero para calcular e, debemos calcular el gradiente de [math]e(\vec{u})[/math] y su transpuesto.

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

[math] \nabla\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}\frac{1}{3}sin(\frac{π}{12}y) & \frac{πx}{36}cos(\frac{π}{12}y) & 0 \\ 0 & 0 & 0\\ 0 & 0 & 0\end{pmatrix} [/math]

Hallamos su transpuesta.

[math]\nabla \vec u ^t= \begin{pmatrix}\frac{1}{3}sin(\frac{π}{12}y) & 0 & 0 \\ \frac{πx}{36}cos(\frac{π}{12}y) & 0 & 0 \\ 0 & 0 & 0\end{pmatrix} [/math]

Encontramos [math]e(\vec{u})[/math] aplicando la fórmula.

[math]e(\vec{u}) = \begin{pmatrix}\frac{1}{3}sin(\frac{π}{12}y) & \frac{πx}{72}cos(\frac{π}{12}y) & 0 \\ \frac{πx}{72}cos(\frac{π}{12}y) & 0 & 0 \\ 0 & 0 & 0\end{pmatrix} [/math]

Teniendo todos los datos, ya podemos calcular nuestro tensor de tensiones [math]σ_{ij}[/math] a través de la fórmula [math]σ = λ∇ · \vec{u} 1 + 2 μe[/math].

[math]σ = λ\frac{1}{3}sin(\frac{π}{12}y)\begin{pmatrix}1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1\end{pmatrix}+2μ\begin{pmatrix}\frac{1}{3}sin(\frac{π}{12}y) & \frac{πx}{72}cos(\frac{π}{12}y) & 0 \\ \frac{πx}{72}cos(\frac{π}{12}y) & 0 & 0 \\ 0 & 0 & 0\end{pmatrix}[/math]

[math]σ = \begin{pmatrix}\frac{1}{3}sin(\frac{π}{12})(2μ+λ) & \frac{πx}{36}cos(\frac{π}{12}y)(μ) & 0 \\ \frac{πx}{36}cos(\frac{π}{12}y)(μ) &\frac{1}{3}sin(\frac{π}{12}y)(λ) & 0 \\ 0 & 0 &\frac{1}{3}sin(\frac{π}{12}y)(λ)\end{pmatrix}[/math]

Teniendo en cuenta que λ=μ=1, obtenemos:

[math]σ = \begin{pmatrix}sin(\frac{π}{12}y) & \frac{πx}{36}cos(\frac{π}{12}y) & 0 \\ \frac{πx}{36}cos(\frac{π}{12}y) &\frac{1}{3}sin(\frac{π}{12}y) & 0 \\ 0 & 0 &\frac{1}{3}sin(\frac{π}{12}y)\end{pmatrix}[/math]

Por último se procede al cálculo de las tensiones normales a la dirección que marca cada vector de la base ortonormal orientada positivamente adoptada. Dichas tensiones son las que se encuentran en los elementos de la diagonal de la matriz [math]\sigma_{ij}[/math]:

8.1 Tensiones normales en dirección del eje [math]\vec i [/math]

Aplicamos: [math]\vec i \cdot \sigma_{ij} \cdot \vec i = sin(\frac{π}{12}y)[/math]

Figura 11. Tension Normal en i

h=2/10; 
x=-1:h:1;
y=0:h:12;
%Generación de la retícula rectangular
[xx,yy]=meshgrid(x,y);
%tensiones normales en la dirección normal del eje i
tni=sin((pi.*yy)/12);
tnj=(1/3).*sin((pi.*yy)/12);
tnk=(1/3).*sin((pi.*yy)/12);
%Representación tensiones normales
figure
surf(xx,yy,tni)
shading flat
 axis equal
colorbar
axis([-1.5,2,-1,15]);
title('Tensión normal en dirección i');
xlabel('Eje X');
 ylabel('Eje Y');
view(2);

8.2 Tensiones normales en dirección del eje [math]\vec j [/math]

Aplicamos: [math]\vec j \cdot \sigma_{ij} \cdot \vec j = \frac{1}{3}sin(\frac{π}{12}y)[/math]

Figura 12. Tension Normal en j

h=2/10; 
x=-1:h:1;
y=0:h:12;
%Generación de la retícula rectangular
[xx,yy]=meshgrid(x,y);
%tensiones normales en la dirección normal del eje j
tni=sin((pi.*yy)/12);
tnj=(1/3).*sin((pi.*yy)/12);
tnk=(1/3).*sin((pi.*yy)/12);
surf(xx,yy,tnj)
shading flat
axis equal
colorbar
axis([-1.5,2,-1,15]);
title('Tensión normal en dirección j');
xlabel('Eje X');
ylabel('Eje Y');
view(2);

8.3 Tensiones normales en dirección del eje [math]\vec k [/math]

Aplicamos: [math]\vec k \cdot \sigma_{ij} \cdot \vec k = \frac{1}{3}sin(\frac{π}{12}y)[/math]

Figura 13. Tension Normal en k

h=2/10; 
x=-1:h:1;
y=0:h:12;
%Generación de la retícula rectangular
[xx,yy]=meshgrid(x,y);
%tensiones normales en la dirección normal del eje k
tni=sin((pi.*yy)/12);
tnj=(1/3).*sin((pi.*yy)/12);
tnk=(1/3).*sin((pi.*yy)/12);
surf(xx,yy,tnk)
shading flat
axis equal
colorbar
axis([-1.5,2,-1,15]);
title('Tensión normal en dirección k');
xlabel('Eje X');
ylabel('Eje Y');
view(2);

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

Sea [math]\vec{T}[/math] la tensión tangencial respecto al plano ortogonal a [math]\vec{i}[/math], que viene dado por la siguiente expresión:[math]\ \left| \sigma ·\vec i -(\vec i·\sigma·\vec i)\vec i \right| [/math]

Con los datos hallados anteriormente obtenemos:

[math]|\begin{pmatrix} sin(\frac{π}{12}y) & \frac{πx}{36}cos(\frac{π}{12}y) & 0 \\ \frac{πx}{36}cos(\frac{π}{12}y) &\frac{1}{3}sin(\frac{π}{12}y) & 0 \\ 0 & 0 &\frac{1}{3}sin(\frac{π}{12}y)\end{pmatrix}[/math]. [math]\begin{pmatrix}1\\ 0\\ 0\end{pmatrix}-\begin{pmatrix}1\\ 0\\ 0\end{pmatrix}\begin{pmatrix} sin(\frac{π}{12}y) & \frac{πx}{36}cos(\frac{π}{12}y) & 0 \\ \frac{πx}{36}cos(\frac{π}{12}y) &\frac{1}{3}sin(\frac{π}{12}y) & 0 \\ 0 & 0 &\frac{1}{3}sin(\frac{π}{12}y)\end{pmatrix}\begin{pmatrix}1\\ 0\\ 0\end{pmatrix}|=0[/math]

10 Tensión de Von Mises

La tensión de Von Mises es una magnitud física proporcional a la energía de distorsión. Debe su nombre a Richard Edler von Mises quien propuso que un material dúctil sufría fallo elástico cuando la energía de distorsión elástica rebasaba cierto valor. En resumen, se trata de una magnitud escalar que se suele usar como indicador para saber cuando un material inicia un comportamiento plástico y no elástico puro.

Esta viene determinada por la siguiente fórmula:

[math] \sigma_V =\sqrt{\frac{(\sigma_1-\sigma_2)^2+(\sigma_2-\sigma_3)^2+(\sigma_3-\sigma_1)^2}{2}}[/math]
σ1,σ2,σ3 son los autovalores del tensor de tensiones σ(ρ,θ).

[math]σ = \begin{pmatrix}σ1 & 0 & 0 \\ 0 &σ2 & 0 \\ 0 & 0 &σ3\end{pmatrix}[/math] = [math]\begin{pmatrix}sin(\frac{π}{12}y) & 0 & 0 \\ 0 &\frac{1}{3}sin(\frac{π}{12}y) & 0 \\ 0 & 0 &\frac{1}{3}sin(\frac{π}{12}y)\end{pmatrix}[/math]

Utilizaremos el siguiente código para el cálculo de la tensión máxima de Von Mises así como la representación de la misma:

Tensión Von Mises

clear
clc
h=0.2;
x=-1:h:1;           
y=0:h:12; 
[Mx,My]=meshgrid(x,y);   %Creamos el mallado

Ti=sin((pi/12).*My);       % Tension normal en i
Tj=(1/3)*sin((pi/12).*My);       % Tension normal en j
Tk=(1/3)*sin((pi/12).*My);       % Tension normal en k

%Autovalores de la tensión de Von Mises
Mz = ones(size(y,2),size(x,2));
T  = zeros(3);
for m=1:length(y)
    for n=1:length(x)
      %Valores para cada componente del mallado para la función de tensión de Von Mises
      a=Ti(m,n);
      b=Tj(m,n);
      c=Tk(m,n);
      T=[a 0 0; 0,b,0; 0,0,c];
      %Obtención de los autovalores de la tensión de Von Mises
      [P,D]=eig(T);
      Autoval1=D(1,1);
      Autoval2=D(2,2);
      Autoval3=D(3,3);
      %Cálculo de Von Mises para cada componente
      VonMises=sqrt( ((Autoval1-Autoval2)^2+(Autoval2-Autoval3)^2+(Autoval3-Autoval1)^2) * 1/2 );
      Mz(m,n)=VonMises;
    end
end

%Representamos Von Mises en la placa bidimensional
figure
surf(Mx,My,Mz)
colorbar
axis equal
xlabel('Eje X')
ylabel('Eje Y')
title('Tensión de Von Mises')
view(2)

%maxima tensión
maxSIGMA=max(M)

fprintf('El valor máximo de la tensión de Von Mises en la placa es %f\n',maxSIGMA)

El valor máximo de la tensión de Von Mises en la placa es 0.666667

11 Campo de fuerzas [math]\vec F[/math]

El campo de fuerzas [math]\vec{F}[/math] que actúa sobre la placa se aproxima usando la ecuación de la elasticidad lineal

[math]\vec{F}=\frac{∂^2\vec{u}}{∂t^2}-∇· σ[/math]

donde [math]∇ · σ[/math] es el campo vectorial que se obtiene al hacer la divergencia de los vectores cuyas componentes son las filas de la matriz [math]σ[/math].

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

Nuestro vector [math]\vec{u}[/math] es [math]\vec{u}=\vec{a}·sin(pi*k(d·\vec{r0}(x, y)−vt))[/math] siendo [math]\vec{a}=x/3 \vec{i} [/math], [math]\vec{d}=1/12 \vec{j}[/math] y [math]k=1[/math].

Sustituyendo las componentes dadas quedaría [math]\vec{u}=\frac{x}{3}sin (\frac{\pi*y}{12}-vt)\vec{i}.[/math]

Calculando [math]∇ · σ[/math], siendo [math]σ = \begin{pmatrix}sin(\frac{π}{12}y-vt) & \frac{πx}{36}cos(\frac{π}{12}y-vt) & 0 \\ \frac{πx}{36}cos(\frac{π}{12}y-vt) &\frac{1}{3}sin(\frac{π}{12}y-vt) & 0 \\ 0 & 0 &\frac{1}{3}sin(\frac{π}{12}y-vt)\end{pmatrix}[/math] obtenemos que [math]∇ · σ = -\frac{xπ}{36}sin(\frac{y\pi}{12}-vt)[/math]

Por otra parte tendremos que calcular [math]\frac{d^2\vec{u}}{dt^2} [/math]

Derivando:

[math]\frac{∂^2\vec{u}}{∂t^2}= \frac{x}{3}·v^2·(-sin(\frac{\pi*y}{12}-vt))\vec{i}.[/math]

Ahora, como sabemos que [math]\vec{F}[/math] = 0 igualamos [math]\frac{d^2\vec{u}}{dt^2} = \nabla \cdot \sigma[/math]

Y con esto hallaremos finalmente que

[math] \vec{v}= \sqrt{\frac{\pi}{12}}\vec{i}[/math]

12 Módulo del desplazamiento transversal

En este último aparatado calcularemos el módulo del desplazamiento transversal a lo largo del tiempo en el intervalo [math]t∈[0, 10][/math] en el punto (1/2,1). Vamos a usar la fórmula obtenida en el apartado anterior para definir la velocidad de propagación.

Finalmente la fórmula con la que nos quedamos es [math]\vec{u}(x,y)=\frac{1}{6}sen(\frac{π}{12}-\sqrt{\frac{\pi}{12}}t)·[/math]

Desplazamiento transversal en función del tiempo

El código que utilizaremos para la representación es el siguiente:

h= 2/10;
t=0:h:10;
f=@(t) (1/6)*sin((pi/12)-(sqrt((pi/12)).*t));
desplazamiento=f(t);
plot(t,desplazamiento)
title ('Representación de la función desplazamiento')