MateWiki - Contribuciones del usuario [es]

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

2023-12-13T19:07:43Z

MiguelMoreno: /* Campo de fuerzas \vec F */

{{ TrabajoED | Visualización de campos escalares y vectoriales en elasticidad (Grupo 16-A) | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC23/24|2023-24]] | Lara Gutiérrez Kreutzer Miguel Moreno Martín Luis Anthony Vera Guapi Emilio Villegas Maroto}}

[[Categoría:Teoría de Campos]]
[[Categoría:TC23/24]]

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.

==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:
[[Archivo:Mallado16a.png|miniaturadeimagen|derecha|Figura 1. Representación del mallado]]

{{matlab|codigo=

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
}}

==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:

[[Archivo:Curvas_nivel_temp1.png|miniaturadeimagen|Figura 2. curvas de nivel del campo de temperaturas ]]

{{matlab|codigo=
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:
[[Archivo:Campo_temp.png|miniaturadeimagen|Figura 3. Campo de temperaturas ]]

{{matlab|codigo=
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.'''

<center> '''INTERPRETACIÓN DEL CAMPO DE TEMPERATURAS''' </center>
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.

=== 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>.

<center><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></center>

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.

[[Archivo:Grad90.png|miniaturadeimagen|350px|Figura 4.Gradiente de Temperatura]]

{{matlab|codigo=

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.
[[Archivo:Grad90_ampliación.png|centre|400px|Figura 5. Ampliación Gradiente de Temperatura]]

<center> '''INTERPRETACIÓN DEL GRADIENTE DEL CAMPO DE TEMPERATURAS''' </center>

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

==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.

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("Gradiente de temperatura");
hold off
}}

==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:
<center><math>\vec{u}(x, y)=x/3 \vec{i} sin (\pi( 1/12 \vec{j}·\vec{r_{0}}(x,y))).</math></center>

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

Por lo que:
<center><math>\vec{u}(x, y)=\frac{x}{3}sin (\frac{\pi*y}{12})\vec{i}.</math></center>

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.

[[Archivo:Campodedesplazamientos16a.png|miniaturadeimagen|derecha|425px|Figura 6. Campo de vectores en t=0.]]

{{matlab|codigo=
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
}}

<center> '''INTERPRETACIÓN DEL CAMPO DE VECTORES''' </center>

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.

==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 <math>\vec{r_{d}}(x, y)=\vec{r_{0}}(x, y)+\vec{u}(x, y).</math>.
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.

[[Archivo:Comparacion16a.png|derecha|miniaturadeimagen|400px|Figura 7. Sólido antes de la deformación, después y comparación.]]

{{matlab|codigo=
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,3,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")
}}

==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 será positiva, y si tiene "sumideros", la divergencia será negativa.

Calculamos la divergencia de <math>\vec u </math> :

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

Entonces:

<center> <math>\nabla ·\vec u = \frac{\partial u_1}{\partial x} + \frac{\partial u_2}{\partial y} = \frac{1}{3}sin(\frac{π}{12}y) </math></center>

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

[[Archivo: Diverge.png|miniaturadeimagen|Divergencia 2D]]
[[Archivo:Diverge2.png|miniaturadeimagen|Divergencia 3D]]

{{matlab|codigo=
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 figura se puede apreciar que la divergencia máxima es de 0.333 y la mínima es también la nula, es decir, 0.'''

== Rotacional del sólido ==
Calcular <math>|∇ × \vec{u}|</math> en todos los puntos del sólido y dibujarlo. ¿Qué puntos sufren un mayor rotacional?

=== Cálculo del rotacional ===
<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>.

=== Representación del módulo del rotacional ===
[[Archivo:Ejercicio27.gif|500px|miniatura de imagen|Módulo del rotacional (campo escalar)]]
{{matlab|codigo=
%Limpieza de programas anteriores
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.5,1.5,-0.5,12.5]);
title('Módulo del rotacional');
xlabel('Eje X');
ylabel('Eje Y');

}}

==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.
<center><math>\vec u(x,y) = \frac{x}{3}\vec{i}sin(\pi(\frac{1}{12}\vec{j}·\vec{r_{0}}(x,y)))</math></center>

Siendo <math>\vec{r_{0}}(x, y)= x \vec{i} + y \vec{j} </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>:

=== 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>

[[Archivo:Tension normal en i.png|400px|miniaturadeimagen|right|Figura 11. Tension Normal en i]]

{{matlab|codigo=
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);
}}

=== 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>

[[Archivo:Tension normal en j.png|400px|miniaturadeimagen|right|Figura 12. Tension Normal en j]]

{{matlab|codigo=
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);
}}

=== 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>

[[Archivo:Tension normal en k1.png|400px|miniaturadeimagen|right|Figura 13. Tension Normal en k]]

{{matlab|codigo=
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);
}}

==Tensiones tangenciales respecto al plano ortogonal a==

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>. \begin{pmatrix}
1\\
0\\
0\end{pmatrix} -sin(frac{π}{12}y)\begin{pmatrix}
1\\
0\\
0\end{pmatrix}| =\frac{πx}{36}cos(frac{π}{12}y)
</math>

==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.

[[Archivo:Von_misesp.jpg|100px|derecha|Richard Edler von Mises]]

Esta viene determinada por la siguiente fórmula:
<center> <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 σ(ρ,θ).</center>

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:

[[Archivo:Tensionvonmises16a.png|400px|miniaturadeimagen|derecha|Tensión Von Mises]]

{{matlab|codigo=
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'''

==Campo de fuerzas <math>\vec F</math>==
El campo de fuerzas <math>\vec{F}</math> que actúa sobre la placa (y que son las causantes del desplazamiento observado) se aproxima usando la ecuación de la elasticidad lineal
<center><math>\vec{F}=\frac{∂^2\vec{u}}{∂t^2}-∇· σ</math></center>
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>.
Calcular la velocidad de propagación de las ondas <math>v</math> en términos de las constantes de Lamé, suponiendo que <math>\vec{F} = 0</math>. Si la onda fuera longitudinal, es decir, tomando <math>\vec{a} = x/3\vec{i}</math>, ¿cuál sería la velocidad de propagación? Comprobar que sobre un mismo medio las ondas transversales y longitudinales no viajan a la misma velocidad, tal y como se observa en la transmisión de ondas sísmicas.

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> queda <math>∇ · σ = \frac{pi}{36}cos (\frac{\pi*y}{12})\vec{j}</math>

Derivando:

<math>\frac{∂\vec{u}}{∂t}= 2/5 \vec{i}·(-v)·cos (\vec{i}+\vec{r0}(x, y) − vt)</math>

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

Sustituyendo en la ecuación de la elasticidad lineal y despejando la velocidad:

<math>\vec{F}=\frac{∂^2\vec{u}}{∂t^2}-∇· σ=[x/3 ·v^2·(-sin(\frac{\pi*y}{12}-vt))\vec{i}-\frac{pi}{36}cos (\frac{\pi*y}{12})\vec{j}=0</math>

==Módulo del desplazamiento transversal==

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

2023-12-13T19:06:04Z

MiguelMoreno: /* Campo de fuerzas \vec F */

{{ TrabajoED | Visualización de campos escalares y vectoriales en elasticidad (Grupo 16-A) | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC23/24|2023-24]] | Lara Gutiérrez Kreutzer Miguel Moreno Martín Luis Anthony Vera Guapi Emilio Villegas Maroto}}

[[Categoría:Teoría de Campos]]
[[Categoría:TC23/24]]

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.

==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:
[[Archivo:Mallado16a.png|miniaturadeimagen|derecha|Figura 1. Representación del mallado]]

{{matlab|codigo=

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
}}

==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:

[[Archivo:Curvas_nivel_temp1.png|miniaturadeimagen|Figura 2. curvas de nivel del campo de temperaturas ]]

{{matlab|codigo=
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:
[[Archivo:Campo_temp.png|miniaturadeimagen|Figura 3. Campo de temperaturas ]]

{{matlab|codigo=
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.'''

<center> '''INTERPRETACIÓN DEL CAMPO DE TEMPERATURAS''' </center>
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.

=== 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>.

<center><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></center>

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.

[[Archivo:Grad90.png|miniaturadeimagen|350px|Figura 4.Gradiente de Temperatura]]

{{matlab|codigo=

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.
[[Archivo:Grad90_ampliación.png|centre|400px|Figura 5. Ampliación Gradiente de Temperatura]]

<center> '''INTERPRETACIÓN DEL GRADIENTE DEL CAMPO DE TEMPERATURAS''' </center>

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

==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.

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("Gradiente de temperatura");
hold off
}}

==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:
<center><math>\vec{u}(x, y)=x/3 \vec{i} sin (\pi( 1/12 \vec{j}·\vec{r_{0}}(x,y))).</math></center>

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

Por lo que:
<center><math>\vec{u}(x, y)=\frac{x}{3}sin (\frac{\pi*y}{12})\vec{i}.</math></center>

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.

[[Archivo:Campodedesplazamientos16a.png|miniaturadeimagen|derecha|425px|Figura 6. Campo de vectores en t=0.]]

{{matlab|codigo=
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
}}

<center> '''INTERPRETACIÓN DEL CAMPO DE VECTORES''' </center>

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.

==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 <math>\vec{r_{d}}(x, y)=\vec{r_{0}}(x, y)+\vec{u}(x, y).</math>.
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.

[[Archivo:Comparacion16a.png|derecha|miniaturadeimagen|400px|Figura 7. Sólido antes de la deformación, después y comparación.]]

{{matlab|codigo=
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,3,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")
}}

==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 será positiva, y si tiene "sumideros", la divergencia será negativa.

Calculamos la divergencia de <math>\vec u </math> :

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

Entonces:

<center> <math>\nabla ·\vec u = \frac{\partial u_1}{\partial x} + \frac{\partial u_2}{\partial y} = \frac{1}{3}sin(\frac{π}{12}y) </math></center>

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

[[Archivo: Diverge.png|miniaturadeimagen|Divergencia 2D]]
[[Archivo:Diverge2.png|miniaturadeimagen|Divergencia 3D]]

{{matlab|codigo=
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 figura se puede apreciar que la divergencia máxima es de 0.333 y la mínima es también la nula, es decir, 0.'''

== Rotacional del sólido ==
Calcular <math>|∇ × \vec{u}|</math> en todos los puntos del sólido y dibujarlo. ¿Qué puntos sufren un mayor rotacional?

=== Cálculo del rotacional ===
<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>.

=== Representación del módulo del rotacional ===
[[Archivo:Ejercicio27.gif|500px|miniatura de imagen|Módulo del rotacional (campo escalar)]]
{{matlab|codigo=
%Limpieza de programas anteriores
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.5,1.5,-0.5,12.5]);
title('Módulo del rotacional');
xlabel('Eje X');
ylabel('Eje Y');

}}

==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.
<center><math>\vec u(x,y) = \frac{x}{3}\vec{i}sin(\pi(\frac{1}{12}\vec{j}·\vec{r_{0}}(x,y)))</math></center>

Siendo <math>\vec{r_{0}}(x, y)= x \vec{i} + y \vec{j} </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>:

=== 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>

[[Archivo:Tension normal en i.png|400px|miniaturadeimagen|right|Figura 11. Tension Normal en i]]

{{matlab|codigo=
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);
}}

=== 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>

[[Archivo:Tension normal en j.png|400px|miniaturadeimagen|right|Figura 12. Tension Normal en j]]

{{matlab|codigo=
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);
}}

=== 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>

[[Archivo:Tension normal en k1.png|400px|miniaturadeimagen|right|Figura 13. Tension Normal en k]]

{{matlab|codigo=
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);
}}

==Tensiones tangenciales respecto al plano ortogonal a==

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>. \begin{pmatrix}
1\\
0\\
0\end{pmatrix} -sin(frac{π}{12}y)\begin{pmatrix}
1\\
0\\
0\end{pmatrix}| =\frac{πx}{36}cos(frac{π}{12}y)
</math>

==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.

[[Archivo:Von_misesp.jpg|100px|derecha|Richard Edler von Mises]]

Esta viene determinada por la siguiente fórmula:
<center> <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 σ(ρ,θ).</center>

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:

[[Archivo:Tensionvonmises16a.png|400px|miniaturadeimagen|derecha|Tensión Von Mises]]

{{matlab|codigo=
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'''

==Campo de fuerzas <math>\vec F</math>==
El campo de fuerzas <math>\vec{F}</math> que actúa sobre la placa (y que son las causantes del desplazamiento observado) se aproxima usando la ecuación de la elasticidad lineal
<center><math>\vec{F}=\frac{∂^2\vec{u}}{∂t^2}-∇· σ</math></center>
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>.
Calcular la velocidad de propagación de las ondas <math>v</math> en términos de las constantes de Lamé, suponiendo que <math>\vec{F} = 0</math>. Si la onda fuera longitudinal, es decir, tomando <math>\vec{a} = x/3\vec{i}</math>, ¿cuál sería la velocidad de propagación? Comprobar que sobre un mismo medio las ondas transversales y longitudinales no viajan a la misma velocidad, tal y como se observa en la transmisión de ondas sísmicas.

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> queda <math>∇ · σ = \frac{x}{36}cos (\frac{\pi*y}{12})\vec{j}</math>

Derivando:

<math>\frac{∂\vec{u}}{∂t}= 2/5 \vec{i}·(-v)·cos (\vec{i}+\vec{r0}(x, y) − vt)</math>

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

Sustituyendo en la ecuación de la elasticidad lineal y despejando la velocidad:

<math>\vec{F}=\frac{∂^2\vec{u}}{∂t^2}-∇· σ=[x/3 \vec{i}·v^2·(-sin(\frac{\pi*y}{12}-vt))\vec{i}-\frac{x}{36}cos (\frac{\pi*y}{12})\vec{j}=0</math>

==Módulo del desplazamiento transversal==

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

2023-12-13T19:03:49Z

MiguelMoreno: /* Campo de fuerzas \vec F */

{{ TrabajoED | Visualización de campos escalares y vectoriales en elasticidad (Grupo 16-A) | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC23/24|2023-24]] | Lara Gutiérrez Kreutzer Miguel Moreno Martín Luis Anthony Vera Guapi Emilio Villegas Maroto}}

[[Categoría:Teoría de Campos]]
[[Categoría:TC23/24]]

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.

==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:
[[Archivo:Mallado16a.png|miniaturadeimagen|derecha|Figura 1. Representación del mallado]]

{{matlab|codigo=

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
}}

==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:

[[Archivo:Curvas_nivel_temp1.png|miniaturadeimagen|Figura 2. curvas de nivel del campo de temperaturas ]]

{{matlab|codigo=
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:
[[Archivo:Campo_temp.png|miniaturadeimagen|Figura 3. Campo de temperaturas ]]

{{matlab|codigo=
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.'''

<center> '''INTERPRETACIÓN DEL CAMPO DE TEMPERATURAS''' </center>
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.

=== 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>.

<center><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></center>

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.

[[Archivo:Grad90.png|miniaturadeimagen|350px|Figura 4.Gradiente de Temperatura]]

{{matlab|codigo=

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.
[[Archivo:Grad90_ampliación.png|centre|400px|Figura 5. Ampliación Gradiente de Temperatura]]

<center> '''INTERPRETACIÓN DEL GRADIENTE DEL CAMPO DE TEMPERATURAS''' </center>

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

==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.

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("Gradiente de temperatura");
hold off
}}

==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:
<center><math>\vec{u}(x, y)=x/3 \vec{i} sin (\pi( 1/12 \vec{j}·\vec{r_{0}}(x,y))).</math></center>

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

Por lo que:
<center><math>\vec{u}(x, y)=\frac{x}{3}sin (\frac{\pi*y}{12})\vec{i}.</math></center>

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.

[[Archivo:Campodedesplazamientos16a.png|miniaturadeimagen|derecha|425px|Figura 6. Campo de vectores en t=0.]]

{{matlab|codigo=
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
}}

<center> '''INTERPRETACIÓN DEL CAMPO DE VECTORES''' </center>

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.

==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 <math>\vec{r_{d}}(x, y)=\vec{r_{0}}(x, y)+\vec{u}(x, y).</math>.
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.

[[Archivo:Comparacion16a.png|derecha|miniaturadeimagen|400px|Figura 7. Sólido antes de la deformación, después y comparación.]]

{{matlab|codigo=
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,3,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")
}}

==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 será positiva, y si tiene "sumideros", la divergencia será negativa.

Calculamos la divergencia de <math>\vec u </math> :

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

Entonces:

<center> <math>\nabla ·\vec u = \frac{\partial u_1}{\partial x} + \frac{\partial u_2}{\partial y} = \frac{1}{3}sin(\frac{π}{12}y) </math></center>

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

[[Archivo: Diverge.png|miniaturadeimagen|Divergencia 2D]]
[[Archivo:Diverge2.png|miniaturadeimagen|Divergencia 3D]]

{{matlab|codigo=
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 figura se puede apreciar que la divergencia máxima es de 0.333 y la mínima es también la nula, es decir, 0.'''

== Rotacional del sólido ==
Calcular <math>|∇ × \vec{u}|</math> en todos los puntos del sólido y dibujarlo. ¿Qué puntos sufren un mayor rotacional?

=== Cálculo del rotacional ===
<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>.

=== Representación del módulo del rotacional ===
[[Archivo:Ejercicio27.gif|500px|miniatura de imagen|Módulo del rotacional (campo escalar)]]
{{matlab|codigo=
%Limpieza de programas anteriores
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.5,1.5,-0.5,12.5]);
title('Módulo del rotacional');
xlabel('Eje X');
ylabel('Eje Y');

}}

==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.
<center><math>\vec u(x,y) = \frac{x}{3}\vec{i}sin(\pi(\frac{1}{12}\vec{j}·\vec{r_{0}}(x,y)))</math></center>

Siendo <math>\vec{r_{0}}(x, y)= x \vec{i} + y \vec{j} </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>:

=== 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>

[[Archivo:Tension normal en i.png|400px|miniaturadeimagen|right|Figura 11. Tension Normal en i]]

{{matlab|codigo=
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);
}}

=== 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>

[[Archivo:Tension normal en j.png|400px|miniaturadeimagen|right|Figura 12. Tension Normal en j]]

{{matlab|codigo=
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);
}}

=== 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>

[[Archivo:Tension normal en k1.png|400px|miniaturadeimagen|right|Figura 13. Tension Normal en k]]

{{matlab|codigo=
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);
}}

==Tensiones tangenciales respecto al plano ortogonal a==

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>. \begin{pmatrix}
1\\
0\\
0\end{pmatrix} -sin(frac{π}{12}y)\begin{pmatrix}
1\\
0\\
0\end{pmatrix}| =\frac{πx}{36}cos(frac{π}{12}y)
</math>

==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.

[[Archivo:Von_misesp.jpg|100px|derecha|Richard Edler von Mises]]

Esta viene determinada por la siguiente fórmula:
<center> <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 σ(ρ,θ).</center>

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:

[[Archivo:Tensionvonmises16a.png|400px|miniaturadeimagen|derecha|Tensión Von Mises]]

{{matlab|codigo=
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'''

==Campo de fuerzas <math>\vec F</math>==
El campo de fuerzas <math>\vec{F}</math> que actúa sobre la placa (y que son las causantes del desplazamiento observado) se aproxima usando la ecuación de la elasticidad lineal
<center><math>\vec{F}=\frac{∂^2\vec{u}}{∂t^2}-∇· σ</math></center>
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>.
Calcular la velocidad de propagación de las ondas <math>v</math> en términos de las constantes de Lamé, suponiendo que <math>\vec{F} = 0</math>. Si la onda fuera longitudinal, es decir, tomando <math>\vec{a} = x/3\vec{i}</math>, ¿cuál sería la velocidad de propagación? Comprobar que sobre un mismo medio las ondas transversales y longitudinales no viajan a la misma velocidad, tal y como se observa en la transmisión de ondas sísmicas.

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> queda <math>∇ · σ = \frac{x}{36}cos (\frac{\pi*y}{12})\vec{j}</math>

Derivando:

<math>\frac{∂\vec{u}}{∂t}= 2/5 \vec{i}·(-v)·cos (\vec{i}+\vec{r0}(x, y) − vt)</math>

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

==Módulo del desplazamiento transversal==

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

2023-12-13T19:02:43Z

MiguelMoreno: /* Campo de fuerzas \vec F */

{{ TrabajoED | Visualización de campos escalares y vectoriales en elasticidad (Grupo 16-A) | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC23/24|2023-24]] | Lara Gutiérrez Kreutzer Miguel Moreno Martín Luis Anthony Vera Guapi Emilio Villegas Maroto}}

[[Categoría:Teoría de Campos]]
[[Categoría:TC23/24]]

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.

==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:
[[Archivo:Mallado16a.png|miniaturadeimagen|derecha|Figura 1. Representación del mallado]]

{{matlab|codigo=

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
}}

==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:

[[Archivo:Curvas_nivel_temp1.png|miniaturadeimagen|Figura 2. curvas de nivel del campo de temperaturas ]]

{{matlab|codigo=
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:
[[Archivo:Campo_temp.png|miniaturadeimagen|Figura 3. Campo de temperaturas ]]

{{matlab|codigo=
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.'''

<center> '''INTERPRETACIÓN DEL CAMPO DE TEMPERATURAS''' </center>
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.

=== 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>.

<center><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></center>

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.

[[Archivo:Grad90.png|miniaturadeimagen|350px|Figura 4.Gradiente de Temperatura]]

{{matlab|codigo=

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.
[[Archivo:Grad90_ampliación.png|centre|400px|Figura 5. Ampliación Gradiente de Temperatura]]

<center> '''INTERPRETACIÓN DEL GRADIENTE DEL CAMPO DE TEMPERATURAS''' </center>

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

==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.

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("Gradiente de temperatura");
hold off
}}

==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:
<center><math>\vec{u}(x, y)=x/3 \vec{i} sin (\pi( 1/12 \vec{j}·\vec{r_{0}}(x,y))).</math></center>

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

Por lo que:
<center><math>\vec{u}(x, y)=\frac{x}{3}sin (\frac{\pi*y}{12})\vec{i}.</math></center>

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.

[[Archivo:Campodedesplazamientos16a.png|miniaturadeimagen|derecha|425px|Figura 6. Campo de vectores en t=0.]]

{{matlab|codigo=
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
}}

<center> '''INTERPRETACIÓN DEL CAMPO DE VECTORES''' </center>

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.

==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 <math>\vec{r_{d}}(x, y)=\vec{r_{0}}(x, y)+\vec{u}(x, y).</math>.
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.

[[Archivo:Comparacion16a.png|derecha|miniaturadeimagen|400px|Figura 7. Sólido antes de la deformación, después y comparación.]]

{{matlab|codigo=
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,3,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")
}}

==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 será positiva, y si tiene "sumideros", la divergencia será negativa.

Calculamos la divergencia de <math>\vec u </math> :

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

Entonces:

<center> <math>\nabla ·\vec u = \frac{\partial u_1}{\partial x} + \frac{\partial u_2}{\partial y} = \frac{1}{3}sin(\frac{π}{12}y) </math></center>

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

[[Archivo: Diverge.png|miniaturadeimagen|Divergencia 2D]]
[[Archivo:Diverge2.png|miniaturadeimagen|Divergencia 3D]]

{{matlab|codigo=
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 figura se puede apreciar que la divergencia máxima es de 0.333 y la mínima es también la nula, es decir, 0.'''

== Rotacional del sólido ==
Calcular <math>|∇ × \vec{u}|</math> en todos los puntos del sólido y dibujarlo. ¿Qué puntos sufren un mayor rotacional?

=== Cálculo del rotacional ===
<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>.

=== Representación del módulo del rotacional ===
[[Archivo:Ejercicio27.gif|500px|miniatura de imagen|Módulo del rotacional (campo escalar)]]
{{matlab|codigo=
%Limpieza de programas anteriores
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.5,1.5,-0.5,12.5]);
title('Módulo del rotacional');
xlabel('Eje X');
ylabel('Eje Y');

}}

==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.
<center><math>\vec u(x,y) = \frac{x}{3}\vec{i}sin(\pi(\frac{1}{12}\vec{j}·\vec{r_{0}}(x,y)))</math></center>

Siendo <math>\vec{r_{0}}(x, y)= x \vec{i} + y \vec{j} </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>:

=== 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>

[[Archivo:Tension normal en i.png|400px|miniaturadeimagen|right|Figura 11. Tension Normal en i]]

{{matlab|codigo=
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);
}}

=== 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>

[[Archivo:Tension normal en j.png|400px|miniaturadeimagen|right|Figura 12. Tension Normal en j]]

{{matlab|codigo=
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);
}}

=== 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>

[[Archivo:Tension normal en k1.png|400px|miniaturadeimagen|right|Figura 13. Tension Normal en k]]

{{matlab|codigo=
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);
}}

==Tensiones tangenciales respecto al plano ortogonal a==

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>. \begin{pmatrix}
1\\
0\\
0\end{pmatrix} -sin(frac{π}{12}y)\begin{pmatrix}
1\\
0\\
0\end{pmatrix}| =\frac{πx}{36}cos(frac{π}{12}y)
</math>

==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.

[[Archivo:Von_misesp.jpg|100px|derecha|Richard Edler von Mises]]

Esta viene determinada por la siguiente fórmula:
<center> <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 σ(ρ,θ).</center>

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:

[[Archivo:Tensionvonmises16a.png|400px|miniaturadeimagen|derecha|Tensión Von Mises]]

{{matlab|codigo=
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'''

==Campo de fuerzas <math>\vec F</math>==
El campo de fuerzas <math>\vec{F}</math> que actúa sobre la placa (y que son las causantes del desplazamiento observado) se aproxima usando la ecuación de la elasticidad lineal
<center><math>\vec{F}=\frac{∂^2\vec{u}}{∂t^2}-∇· σ</math></center>
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>.
Calcular la velocidad de propagación de las ondas <math>v</math> en términos de las constantes de Lamé, suponiendo que <math>\vec{F} = 0</math>. Si la onda fuera longitudinal, es decir, tomando <math>\vec{a} = x/3\vec{i}</math>, ¿cuál sería la velocidad de propagación? Comprobar que sobre un mismo medio las ondas transversales y longitudinales no viajan a la misma velocidad, tal y como se observa en la transmisión de ondas sísmicas.

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> queda <math>∇ · σ = frac{x}{36}cos (\frac{\pi*y}{12})\vec{j}</math>

Derivando:

<math>\frac{∂\vec{u}}{∂t}= 2/5 \vec{i}·(-v)·cos (\vec{i}+\vec{r0}(x, y) − vt)</math>

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

==Módulo del desplazamiento transversal==

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

2023-12-13T18:59:41Z

MiguelMoreno: /* Campo de fuerzas \vec F */

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

2023-12-13T18:59:07Z

MiguelMoreno: /* Campo de fuerzas \vec F */

{{ TrabajoED | Visualización de campos escalares y vectoriales en elasticidad (Grupo 16-A) | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC23/24|2023-24]] | Lara Gutiérrez Kreutzer Miguel Moreno Martín Luis Anthony Vera Guapi Emilio Villegas Maroto}}

[[Categoría:Teoría de Campos]]
[[Categoría:TC23/24]]

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.

==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:
[[Archivo:Mallado16a.png|miniaturadeimagen|derecha|Figura 1. Representación del mallado]]

{{matlab|codigo=

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
}}

==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:

[[Archivo:Curvas_nivel_temp1.png|miniaturadeimagen|Figura 2. curvas de nivel del campo de temperaturas ]]

{{matlab|codigo=
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:
[[Archivo:Campo_temp.png|miniaturadeimagen|Figura 3. Campo de temperaturas ]]

{{matlab|codigo=
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.'''

<center> '''INTERPRETACIÓN DEL CAMPO DE TEMPERATURAS''' </center>
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.

=== 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>.

<center><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></center>

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.

[[Archivo:Grad90.png|miniaturadeimagen|350px|Figura 4.Gradiente de Temperatura]]

{{matlab|codigo=

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.
[[Archivo:Grad90_ampliación.png|centre|400px|Figura 5. Ampliación Gradiente de Temperatura]]

<center> '''INTERPRETACIÓN DEL GRADIENTE DEL CAMPO DE TEMPERATURAS''' </center>

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

==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.

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("Gradiente de temperatura");
hold off
}}

==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:
<center><math>\vec{u}(x, y)=x/3 \vec{i} sin (\pi( 1/12 \vec{j}·\vec{r_{0}}(x,y))).</math></center>

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

Por lo que:
<center><math>\vec{u}(x, y)=\frac{x}{3}sin (\frac{\pi*y}{12})\vec{i}.</math></center>

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.

[[Archivo:Campodedesplazamientos16a.png|miniaturadeimagen|derecha|425px|Figura 6. Campo de vectores en t=0.]]

{{matlab|codigo=
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
}}

<center> '''INTERPRETACIÓN DEL CAMPO DE VECTORES''' </center>

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.

==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 <math>\vec{r_{d}}(x, y)=\vec{r_{0}}(x, y)+\vec{u}(x, y).</math>.
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.

[[Archivo:Comparacion16a.png|derecha|miniaturadeimagen|400px|Figura 7. Sólido antes de la deformación, después y comparación.]]

{{matlab|codigo=
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,3,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")
}}

==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 será positiva, y si tiene "sumideros", la divergencia será negativa.

Calculamos la divergencia de <math>\vec u </math> :

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

Entonces:

<center> <math>\nabla ·\vec u = \frac{\partial u_1}{\partial x} + \frac{\partial u_2}{\partial y} = \frac{1}{3}sin(\frac{π}{12}y) </math></center>

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

[[Archivo: Diverge.png|miniaturadeimagen|Divergencia 2D]]
[[Archivo:Diverge2.png|miniaturadeimagen|Divergencia 3D]]

{{matlab|codigo=
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 figura se puede apreciar que la divergencia máxima es de 0.333 y la mínima es también la nula, es decir, 0.'''

== Rotacional del sólido ==
Calcular <math>|∇ × \vec{u}|</math> en todos los puntos del sólido y dibujarlo. ¿Qué puntos sufren un mayor rotacional?

=== Cálculo del rotacional ===
<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>.

=== Representación del módulo del rotacional ===
[[Archivo:Ejercicio27.gif|500px|miniatura de imagen|Módulo del rotacional (campo escalar)]]
{{matlab|codigo=
%Limpieza de programas anteriores
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.5,1.5,-0.5,12.5]);
title('Módulo del rotacional');
xlabel('Eje X');
ylabel('Eje Y');

}}

==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.
<center><math>\vec u(x,y) = \frac{x}{3}\vec{i}sin(\pi(\frac{1}{12}\vec{j}·\vec{r_{0}}(x,y)))</math></center>

Siendo <math>\vec{r_{0}}(x, y)= x \vec{i} + y \vec{j} </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>:

=== 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>

[[Archivo:Tension normal en i.png|400px|miniaturadeimagen|right|Figura 11. Tension Normal en i]]

{{matlab|codigo=
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);
}}

=== 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>

[[Archivo:Tension normal en j.png|400px|miniaturadeimagen|right|Figura 12. Tension Normal en j]]

{{matlab|codigo=
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);
}}

=== 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>

[[Archivo:Tension normal en k1.png|400px|miniaturadeimagen|right|Figura 13. Tension Normal en k]]

{{matlab|codigo=
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);
}}

==Tensiones tangenciales respecto al plano ortogonal a==

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>. \begin{pmatrix}
1\\
0\\
0\end{pmatrix} -sin(frac{π}{12}y)\begin{pmatrix}
1\\
0\\
0\end{pmatrix}| =\frac{πx}{36}cos(frac{π}{12}y)
</math>

==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.

[[Archivo:Von_misesp.jpg|100px|derecha|Richard Edler von Mises]]

Esta viene determinada por la siguiente fórmula:
<center> <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 σ(ρ,θ).</center>

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:

[[Archivo:Tensionvonmises16a.png|400px|miniaturadeimagen|derecha|Tensión Von Mises]]

{{matlab|codigo=
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'''

==Campo de fuerzas <math>\vec F</math>==
El campo de fuerzas <math>\vec{F}</math> que actúa sobre la placa (y que son las causantes del desplazamiento observado) se aproxima usando la ecuación de la elasticidad lineal
<center><math>\vec{F}=\frac{∂^2\vec{u}}{∂t^2}-∇· σ</math></center>
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>.
Calcular la velocidad de propagación de las ondas <math>v</math> en términos de las constantes de Lamé, suponiendo que <math>\vec{F} = 0</math>. Si la onda fuera longitudinal, es decir, tomando <math>\vec{a} = x/3\vec{i}</math>, ¿cuál sería la velocidad de propagación? Comprobar que sobre un mismo medio las ondas transversales y longitudinales no viajan a la misma velocidad, tal y como se observa en la transmisión de ondas sísmicas.

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></center>

==Módulo del desplazamiento transversal==

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

2023-12-13T18:57:46Z

MiguelMoreno: /* Campo de fuerzas \vec F */

{{ TrabajoED | Visualización de campos escalares y vectoriales en elasticidad (Grupo 16-A) | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC23/24|2023-24]] | Lara Gutiérrez Kreutzer Miguel Moreno Martín Luis Anthony Vera Guapi Emilio Villegas Maroto}}

[[Categoría:Teoría de Campos]]
[[Categoría:TC23/24]]

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.

==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:
[[Archivo:Mallado16a.png|miniaturadeimagen|derecha|Figura 1. Representación del mallado]]

{{matlab|codigo=

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
}}

==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:

[[Archivo:Curvas_nivel_temp1.png|miniaturadeimagen|Figura 2. curvas de nivel del campo de temperaturas ]]

{{matlab|codigo=
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:
[[Archivo:Campo_temp.png|miniaturadeimagen|Figura 3. Campo de temperaturas ]]

{{matlab|codigo=
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.'''

<center> '''INTERPRETACIÓN DEL CAMPO DE TEMPERATURAS''' </center>
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.

=== 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>.

<center><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></center>

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.

[[Archivo:Grad90.png|miniaturadeimagen|350px|Figura 4.Gradiente de Temperatura]]

{{matlab|codigo=

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.
[[Archivo:Grad90_ampliación.png|centre|400px|Figura 5. Ampliación Gradiente de Temperatura]]

<center> '''INTERPRETACIÓN DEL GRADIENTE DEL CAMPO DE TEMPERATURAS''' </center>

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

==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.

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("Gradiente de temperatura");
hold off
}}

==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:
<center><math>\vec{u}(x, y)=x/3 \vec{i} sin (\pi( 1/12 \vec{j}·\vec{r_{0}}(x,y))).</math></center>

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

Por lo que:
<center><math>\vec{u}(x, y)=\frac{x}{3}sin (\frac{\pi*y}{12})\vec{i}.</math></center>

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.

[[Archivo:Campodedesplazamientos16a.png|miniaturadeimagen|derecha|425px|Figura 6. Campo de vectores en t=0.]]

{{matlab|codigo=
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
}}

<center> '''INTERPRETACIÓN DEL CAMPO DE VECTORES''' </center>

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.

==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 <math>\vec{r_{d}}(x, y)=\vec{r_{0}}(x, y)+\vec{u}(x, y).</math>.
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.

[[Archivo:Comparacion16a.png|derecha|miniaturadeimagen|400px|Figura 7. Sólido antes de la deformación, después y comparación.]]

{{matlab|codigo=
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,3,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")
}}

==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 será positiva, y si tiene "sumideros", la divergencia será negativa.

Calculamos la divergencia de <math>\vec u </math> :

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

Entonces:

<center> <math>\nabla ·\vec u = \frac{\partial u_1}{\partial x} + \frac{\partial u_2}{\partial y} = \frac{1}{3}sin(\frac{π}{12}y) </math></center>

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

[[Archivo: Diverge.png|miniaturadeimagen|Divergencia 2D]]
[[Archivo:Diverge2.png|miniaturadeimagen|Divergencia 3D]]

{{matlab|codigo=
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 figura se puede apreciar que la divergencia máxima es de 0.333 y la mínima es también la nula, es decir, 0.'''

== Rotacional del sólido ==
Calcular <math>|∇ × \vec{u}|</math> en todos los puntos del sólido y dibujarlo. ¿Qué puntos sufren un mayor rotacional?

=== Cálculo del rotacional ===
<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>.

=== Representación del módulo del rotacional ===
[[Archivo:Ejercicio27.gif|500px|miniatura de imagen|Módulo del rotacional (campo escalar)]]
{{matlab|codigo=
%Limpieza de programas anteriores
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.5,1.5,-0.5,12.5]);
title('Módulo del rotacional');
xlabel('Eje X');
ylabel('Eje Y');

}}

==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.
<center><math>\vec u(x,y) = \frac{x}{3}\vec{i}sin(\pi(\frac{1}{12}\vec{j}·\vec{r_{0}}(x,y)))</math></center>

Siendo <math>\vec{r_{0}}(x, y)= x \vec{i} + y \vec{j} </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>:

=== 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>

[[Archivo:Tension normal en i.png|400px|miniaturadeimagen|right|Figura 11. Tension Normal en i]]

{{matlab|codigo=
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);
}}

=== 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>

[[Archivo:Tension normal en j.png|400px|miniaturadeimagen|right|Figura 12. Tension Normal en j]]

{{matlab|codigo=
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);
}}

=== 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>

[[Archivo:Tension normal en k1.png|400px|miniaturadeimagen|right|Figura 13. Tension Normal en k]]

{{matlab|codigo=
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);
}}

==Tensiones tangenciales respecto al plano ortogonal a==

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>. \begin{pmatrix}
1\\
0\\
0\end{pmatrix} -sin(frac{π}{12}y)\begin{pmatrix}
1\\
0\\
0\end{pmatrix}| =\frac{πx}{36}cos(frac{π}{12}y)
</math>

==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.

[[Archivo:Von_misesp.jpg|100px|derecha|Richard Edler von Mises]]

Esta viene determinada por la siguiente fórmula:
<center> <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 σ(ρ,θ).</center>

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:

[[Archivo:Tensionvonmises16a.png|400px|miniaturadeimagen|derecha|Tensión Von Mises]]

{{matlab|codigo=
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'''

==Campo de fuerzas <math>\vec F</math>==
El campo de fuerzas <math>\vec{F}</math> que actúa sobre la placa (y que son las causantes del desplazamiento observado) se aproxima usando la ecuación de la elasticidad lineal
<center><math>\vec{F}=\frac{∂^2\vec{u}}{∂t^2}-∇· σ</math></center>
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>.
Calcular la velocidad de propagación de las ondas <math>v</math> en términos de las constantes de Lamé, suponiendo que <math>\vec{F} = 0</math>. Si la onda fuera longitudinal, es decir, tomando <math>\vec{a} = x/3\vec{i}</math>, ¿cuál sería la velocidad de propagación? Comprobar que sobre un mismo medio las ondas transversales y longitudinales no viajan a la misma velocidad, tal y como se observa en la transmisión de ondas sísmicas.

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}=x/3\vec{i}·sin(\vec{i}+\vec{r0}(x, y)−vt)</math>

==Módulo del desplazamiento transversal==

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

2023-12-13T18:44:41Z

MiguelMoreno: /* Campo de fuerzas \vec F */

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

2023-12-13T18:14:45Z

MiguelMoreno: /* Representación del módulo del rotacional */

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

2023-12-13T18:12:01Z

MiguelMoreno: /* Ley de Fourier */

{{ TrabajoED | Visualización de campos escalares y vectoriales en elasticidad (Grupo 16-A) | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC23/24|2023-24]] | Lara Gutiérrez Kreutzer Miguel Moreno Martín Luis Anthony Vera Guapi Emilio Villegas Maroto}}

[[Categoría:Teoría de Campos]]
[[Categoría:TC23/24]]

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.

==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:
[[Archivo:Mallado16a.png|miniaturadeimagen|derecha|Figura 1. Representación del mallado]]

{{matlab|codigo=

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
}}

==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:

[[Archivo:Curvas_nivel_temp1.png|miniaturadeimagen|Figura 2. curvas de nivel del campo de temperaturas ]]

{{matlab|codigo=
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:
[[Archivo:Campo_temp.png|miniaturadeimagen|Figura 3. Campo de temperaturas ]]

{{matlab|codigo=
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.'''

<center> '''INTERPRETACIÓN DEL CAMPO DE TEMPERATURAS''' </center>
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.

=== 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>.

<center><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></center>

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.

[[Archivo:Grad90.png|miniaturadeimagen|350px|Figura 4.Gradiente de Temperatura]]

{{matlab|codigo=

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.
[[Archivo:Grad90_ampliación.png|centre|400px|Figura 5. Ampliación Gradiente de Temperatura]]

<center> '''INTERPRETACIÓN DEL GRADIENTE DEL CAMPO DE TEMPERATURAS''' </center>

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

==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.

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("Gradiente de temperatura");
hold off
}}

==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:
<center><math>\vec{u}(x, y)=x/3 \vec{i} sin (\pi( 1/12 \vec{j}·\vec{r_{0}}(x,y))).</math></center>

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

Por lo que:
<center><math>\vec{u}(x, y)=\frac{x}{3}sin (\frac{\pi*y}{12})\vec{i}.</math></center>

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.

[[Archivo:Campodedesplazamientos16a.png|miniaturadeimagen|derecha|425px|Figura 6. Campo de vectores en t=0.]]

{{matlab|codigo=
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
}}

<center> '''INTERPRETACIÓN DEL CAMPO DE VECTORES''' </center>

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.

==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 <math>\vec{r_{d}}(x, y)=\vec{r_{0}}(x, y)+\vec{u}(x, y).</math>.
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.

[[Archivo:Comparacion16a.png|derecha|miniaturadeimagen|400px|Figura 7. Sólido antes de la deformación, después y comparación.]]

{{matlab|codigo=
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,3,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")
}}

==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 será positiva, y si tiene "sumideros", la divergencia será negativa.

Calculamos la divergencia de <math>\vec u </math> :

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

Entonces:

<center> <math>\nabla ·\vec u = \frac{\partial u_1}{\partial x} + \frac{\partial u_2}{\partial y} = \frac{1}{3}sin(\frac{π}{12}y) </math></center>

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

[[Archivo: Diverge.png|miniaturadeimagen|Divergencia 2D]]
[[Archivo:Diverge2.png|miniaturadeimagen|Divergencia 3D]]

{{matlab|codigo=
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 figura se puede apreciar que la divergencia máxima es de 0.333 y la mínima es también la nula, es decir, 0.'''

== Rotacional del sólido ==
Calcular <math>|∇ × \vec{u}|</math> en todos los puntos del sólido y dibujarlo. ¿Qué puntos sufren un mayor rotacional?

=== Cálculo del rotacional ===
<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>.

=== Representación del módulo del rotacional ===
[[Archivo:Ejercicio27.gif|500px|miniatura de imagen|Módulo del rotacional (campo escalar)]]
{{matlab|codigo=
%Limpieza de programas anteriores
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.*My)/36).*cos((pi*My)/12);
%Representación gráfica
surf(Mx,My,rot)
shading flat
axis equal
colorbar
view(2);
axis([-1.5,1.5,-0.5,12.5]);
title('Módulo del rotacional');
xlabel('Eje X');
ylabel('Eje Y');

}}

==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.
<center><math>\vec u(x,y) = \frac{x}{3}\vec{i}sin(\pi(\frac{1}{12}\vec{j}·\vec{r_{0}}(x,y)))</math></center>

Siendo <math>\vec{r_{0}}(x, y)= x \vec{i} + y \vec{j} </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>:

=== 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>

[[Archivo:Tension normal en i.png|400px|miniaturadeimagen|right|Figura 11. Tension Normal en i]]

{{matlab|codigo=
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);
}}

=== 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>

[[Archivo:Tension normal en j.png|400px|miniaturadeimagen|right|Figura 12. Tension Normal en j]]

{{matlab|codigo=
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);
}}

=== 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>

[[Archivo:Tension normal en k1.png|400px|miniaturadeimagen|right|Figura 13. Tension Normal en k]]

{{matlab|codigo=
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);
}}

==Tensiones tangenciales respecto al plano ortogonal a==

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>. \begin{pmatrix}
1\\
0\\
0\end{pmatrix} -sin(frac{π}{12}y)\begin{pmatrix}
1\\
0\\
0\end{pmatrix}| =\frac{πx}{36}cos(frac{π}{12}y)
</math>

==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.

[[Archivo:Von_misesp.jpg|100px|derecha|Richard Edler von Mises]]

Esta viene determinada por la siguiente fórmula:
<center> <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 σ(ρ,θ).</center>

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:

[[Archivo:Tensionvonmises16a.png|400px|miniaturadeimagen|derecha|Tensión Von Mises]]

{{matlab|codigo=
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'''

==Campo de fuerzas <math>\vec F</math>==
==Módulo del desplazamiento transversal==

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

2023-12-13T18:02:04Z

MiguelMoreno: /* Ley de Fourier */

{{ TrabajoED | Visualización de campos escalares y vectoriales en elasticidad (Grupo 16-A) | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC23/24|2023-24]] | Lara Gutiérrez Kreutzer Miguel Moreno Martín Luis Anthony Vera Guapi Emilio Villegas Maroto}}

[[Categoría:Teoría de Campos]]
[[Categoría:TC23/24]]

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.

==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:
[[Archivo:Mallado16a.png|miniaturadeimagen|derecha|Figura 1. Representación del mallado]]

{{matlab|codigo=

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
}}

==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:

[[Archivo:Curvas_nivel_temp1.png|miniaturadeimagen|Figura 2. curvas de nivel del campo de temperaturas ]]

{{matlab|codigo=
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:
[[Archivo:Campo_temp.png|miniaturadeimagen|Figura 3. Campo de temperaturas ]]

{{matlab|codigo=
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.'''

<center> '''INTERPRETACIÓN DEL CAMPO DE TEMPERATURAS''' </center>
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.

=== 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>.

<center><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></center>

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.

[[Archivo:Grad90.png|miniaturadeimagen|350px|Figura 4.Gradiente de Temperatura]]

{{matlab|codigo=

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.
[[Archivo:Grad90_ampliación.png|centre|400px|Figura 5. Ampliación Gradiente de Temperatura]]

<center> '''INTERPRETACIÓN DEL GRADIENTE DEL CAMPO DE TEMPERATURAS''' </center>

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

==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.

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
}}

==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:
<center><math>\vec{u}(x, y)=x/3 \vec{i} sin (\pi( 1/12 \vec{j}·\vec{r_{0}}(x,y))).</math></center>

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

Por lo que:
<center><math>\vec{u}(x, y)=\frac{x}{3}sin (\frac{\pi*y}{12})\vec{i}.</math></center>

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.

[[Archivo:Campodedesplazamientos16a.png|miniaturadeimagen|derecha|425px|Figura 6. Campo de vectores en t=0.]]

{{matlab|codigo=
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
}}

<center> '''INTERPRETACIÓN DEL CAMPO DE VECTORES''' </center>

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.

==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 <math>\vec{r_{d}}(x, y)=\vec{r_{0}}(x, y)+\vec{u}(x, y).</math>.
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.

[[Archivo:Comparacion16a.png|derecha|miniaturadeimagen|400px|Figura 7. Sólido antes de la deformación, después y comparación.]]

{{matlab|codigo=
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,3,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")
}}

==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 será positiva, y si tiene "sumideros", la divergencia será negativa.

Calculamos la divergencia de <math>\vec u </math> :

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

Entonces:

<center> <math>\nabla ·\vec u = \frac{\partial u_1}{\partial x} + \frac{\partial u_2}{\partial y} = \frac{1}{3}sin(\frac{π}{12}y) </math></center>

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

[[Archivo: Diverge.png|miniaturadeimagen|Divergencia 2D]]
[[Archivo:Diverge2.png|miniaturadeimagen|Divergencia 3D]]

{{matlab|codigo=
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 figura se puede apreciar que la divergencia máxima es de 0.333 y la mínima es también la nula, es decir, 0.'''

== Rotacional del sólido ==
Calcular <math>|∇ × \vec{u}|</math> en todos los puntos del sólido y dibujarlo. ¿Qué puntos sufren un mayor rotacional?

=== Cálculo del rotacional ===
<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>.

=== Representación del módulo del rotacional ===
[[Archivo:Ejercicio27.gif|500px|miniatura de imagen|Módulo del rotacional (campo escalar)]]
{{matlab|codigo=
%Limpieza de programas anteriores
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.*My)/36).*cos((pi*My)/12);
%Representación gráfica
surf(Mx,My,rot)
shading flat
axis equal
colorbar
view(2);
axis([-1.5,1.5,-0.5,12.5]);
title('Módulo del rotacional');
xlabel('Eje X');
ylabel('Eje Y');

}}

==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.
<center><math>\vec u(x,y) = \frac{x}{3}\vec{i}sin(\pi(\frac{1}{12}\vec{j}·\vec{r_{0}}(x,y)))</math></center>

Siendo <math>\vec{r_{0}}(x, y)= x \vec{i} + y \vec{j} </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>:

=== 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>

[[Archivo:Tension normal en i.png|400px|miniaturadeimagen|right|Figura 11. Tension Normal en i]]

{{matlab|codigo=
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);
}}

=== 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>

[[Archivo:Tension normal en j.png|400px|miniaturadeimagen|right|Figura 12. Tension Normal en j]]

{{matlab|codigo=
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);
}}

=== 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>

[[Archivo:Tension normal en k1.png|400px|miniaturadeimagen|right|Figura 13. Tension Normal en k]]

{{matlab|codigo=
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);
}}

==Tensiones tangenciales respecto al plano ortogonal a==

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>. \begin{pmatrix}
1\\
0\\
0\end{pmatrix} -sin(frac{π}{12}y)\begin{pmatrix}
1\\
0\\
0\end{pmatrix}| =\frac{πx}{36}cos(frac{π}{12}y)
</math>

==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.

[[Archivo:Von_misesp.jpg|100px|derecha|Richard Edler von Mises]]

Esta viene determinada por la siguiente fórmula:
<center> <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 σ(ρ,θ).</center>

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:

[[Archivo:Tensionvonmises16a.png|400px|miniaturadeimagen|derecha|Tensión Von Mises]]

{{matlab|codigo=
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'''

==Campo de fuerzas <math>\vec F</math>==
==Módulo del desplazamiento transversal==

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

2023-12-13T18:01:08Z

MiguelMoreno:

{{ TrabajoED | Visualización de campos escalares y vectoriales en elasticidad (Grupo 16-A) | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC23/24|2023-24]] | Lara Gutiérrez Kreutzer Miguel Moreno Martín Luis Anthony Vera Guapi Emilio Villegas Maroto}}

[[Categoría:Teoría de Campos]]
[[Categoría:TC23/24]]

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.

==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:
[[Archivo:Mallado16a.png|miniaturadeimagen|derecha|Figura 1. Representación del mallado]]

{{matlab|codigo=

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
}}

==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:

[[Archivo:Curvas_nivel_temp1.png|miniaturadeimagen|Figura 2. curvas de nivel del campo de temperaturas ]]

{{matlab|codigo=
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:
[[Archivo:Campo_temp.png|miniaturadeimagen|Figura 3. Campo de temperaturas ]]

{{matlab|codigo=
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.'''

<center> '''INTERPRETACIÓN DEL CAMPO DE TEMPERATURAS''' </center>
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.

=== 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>.

<center><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></center>

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.

[[Archivo:Grad90.png|miniaturadeimagen|350px|Figura 4.Gradiente de Temperatura]]

{{matlab|codigo=

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.
[[Archivo:Grad90_ampliación.png|centre|400px|Figura 5. Ampliación Gradiente de Temperatura]]

<center> '''INTERPRETACIÓN DEL GRADIENTE DEL CAMPO DE TEMPERATURAS''' </center>

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

==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.

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
}}

==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:
<center><math>\vec{u}(x, y)=x/3 \vec{i} sin (\pi( 1/12 \vec{j}·\vec{r_{0}}(x,y))).</math></center>

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

Por lo que:
<center><math>\vec{u}(x, y)=\frac{x}{3}sin (\frac{\pi*y}{12})\vec{i}.</math></center>

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.

[[Archivo:Campodedesplazamientos16a.png|miniaturadeimagen|derecha|425px|Figura 6. Campo de vectores en t=0.]]

{{matlab|codigo=
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
}}

<center> '''INTERPRETACIÓN DEL CAMPO DE VECTORES''' </center>

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.

==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 <math>\vec{r_{d}}(x, y)=\vec{r_{0}}(x, y)+\vec{u}(x, y).</math>.
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.

[[Archivo:Comparacion16a.png|derecha|miniaturadeimagen|400px|Figura 7. Sólido antes de la deformación, después y comparación.]]

{{matlab|codigo=
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,3,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")
}}

==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 será positiva, y si tiene "sumideros", la divergencia será negativa.

Calculamos la divergencia de <math>\vec u </math> :

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

Entonces:

<center> <math>\nabla ·\vec u = \frac{\partial u_1}{\partial x} + \frac{\partial u_2}{\partial y} = \frac{1}{3}sin(\frac{π}{12}y) </math></center>

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

[[Archivo: Diverge.png|miniaturadeimagen|Divergencia 2D]]
[[Archivo:Diverge2.png|miniaturadeimagen|Divergencia 3D]]

{{matlab|codigo=
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 figura se puede apreciar que la divergencia máxima es de 0.333 y la mínima es también la nula, es decir, 0.'''

== Rotacional del sólido ==
Calcular <math>|∇ × \vec{u}|</math> en todos los puntos del sólido y dibujarlo. ¿Qué puntos sufren un mayor rotacional?

=== Cálculo del rotacional ===
<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>.

=== Representación del módulo del rotacional ===
[[Archivo:Ejercicio27.gif|500px|miniatura de imagen|Módulo del rotacional (campo escalar)]]
{{matlab|codigo=
%Limpieza de programas anteriores
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.*My)/36).*cos((pi*My)/12);
%Representación gráfica
surf(Mx,My,rot)
shading flat
axis equal
colorbar
view(2);
axis([-1.5,1.5,-0.5,12.5]);
title('Módulo del rotacional');
xlabel('Eje X');
ylabel('Eje Y');

}}

==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.
<center><math>\vec u(x,y) = \frac{x}{3}\vec{i}sin(\pi(\frac{1}{12}\vec{j}·\vec{r_{0}}(x,y)))</math></center>

Siendo <math>\vec{r_{0}}(x, y)= x \vec{i} + y \vec{j} </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>:

=== 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>

[[Archivo:Tension normal en i.png|400px|miniaturadeimagen|right|Figura 11. Tension Normal en i]]

{{matlab|codigo=
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);
}}

=== 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>

[[Archivo:Tension normal en j.png|400px|miniaturadeimagen|right|Figura 12. Tension Normal en j]]

{{matlab|codigo=
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);
}}

=== 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>

[[Archivo:Tension normal en k1.png|400px|miniaturadeimagen|right|Figura 13. Tension Normal en k]]

{{matlab|codigo=
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);
}}

==Tensiones tangenciales respecto al plano ortogonal a==

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>. \begin{pmatrix}
1\\
0\\
0\end{pmatrix} -sin(frac{π}{12}y)\begin{pmatrix}
1\\
0\\
0\end{pmatrix}| =\frac{πx}{36}cos(frac{π}{12}y)
</math>

==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.

[[Archivo:Von_misesp.jpg|100px|derecha|Richard Edler von Mises]]

Esta viene determinada por la siguiente fórmula:
<center> <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 σ(ρ,θ).</center>

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:

[[Archivo:Tensionvonmises16a.png|400px|miniaturadeimagen|derecha|Tensión Von Mises]]

{{matlab|codigo=
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'''

==Campo de fuerzas <math>\vec F</math>==
==Módulo del desplazamiento transversal==

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

2023-12-09T10:19:14Z

MiguelMoreno: /* Rotacional del sólido */

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

2023-12-09T10:09:54Z

MiguelMoreno: /* Gradiente de la temperatura */

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

2023-12-09T10:08:08Z

MiguelMoreno: