Visualización de campos escalares y vectoriales en elasticidad. (Grupo 28)

De MateWiki
Saltar a: navegación, buscar
Trabajo realizado por estudiantes
Título Visualización de Campos escalares y vectoriales en Placa rectangular plana (en dimensión 2D) (Grupo 28)
Asignatura Teoría de Campos
Curso 2023-24
Autores Carlos De la Torre Pérez
Marcos Herraez
Pablo Jauregui
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

En este artículo se estudiarán los resultados que provocan varios campos escalares y vectoriales sobre un sólido. Para ello, haremos uso de MatLab.
Consideramos una placa rectangular plana que ocupa la región [math](x, y) ∈ [−1, 1] × [0, 12][/math]
En ella vamos a suponer que tenemos definidas dos cantidades físicas:

  • La temperatura
  • Los desplazamientos producidos por la acción de una fuerza determinada.
La temperatura [math]T(x, y)[/math] viene dada por la función:
[math]T(x, y)=3log(1+(x+1)^2)+log(1+(y-2)^2)[/math]

y la posición de cada punto [math](x, y)[/math] de la placa después de la deformación viene dada por:
[math]\vec{r_{d}}(x, y) = \vec{r_{0}}(x, y) + \vec{u}(x, y)[/math]

Vamos a suponer que la fuerza aplicada sobre la placa ha provocado un desplazamiento ondulatorio de los puntos de la misma dado por el vector:
[math]\vec{u}(x, y, t) = \vec{a} sin(πk(\vec{d} · \vec{r_{0}}(x, y) − vt)),[/math]

donde [math]\vec{a}[/math] se conoce como amplitud, [math]k \gt 0[/math] es el número de onda, [math]\vec{d}[/math] es un vector unitario que marca la dirección de propagación y [math]v[/math] es la velocidad de propagación.
Por último, supondremos que los desplazamientos se corresponden a una onda longitudinal en la que la dirección de propagación es la misma que la amplitud. Tomaremos
[math]\vec{a} = \vec{d} = 1/3 \vec{j}, k = 1[/math]


2 Mallado de la placa

Figura 1: Mallado de la placa rectangular

A continuación, se mostrará el mallado de la placa y el código utilizado en MatLab para obtenerlo.
Acto seguido se procederá a una breve explicación del código empleado:

  • En la primera línea del código se utilizan esos comandos para limpiar los programas anteriores y que el nuestro se ejecute de forma correcta.
  • En los dos siguientes comandos discretizamos nuestras variables.
  • En la cuarta y en la quinta línea escribimos los comandos para crear el mallado.
  • Los cuatro siguientes comandos definirán nuestros ejes, mientras que los últimos dos sirven para ponerle un título a nuestra gráfica y visualizar la planta de nuestra placa .
clc,clear,clf 
x1=-1:0.2:1;
y1=0:0.2:12;
[x2,y2]=meshgrid(x1,y1);
mesh(x2,y2,y2*0)
axis equal;
axis([-5,5,-1,13]);
xlabel('Eje X');
ylabel('Eje Y');
title('Placa rectangular')
view(2);


3 Curvas de nivel de la temperatura de la placa

Figura 2: Curvas de nivel de la temperatura y su gradiente

Ahora se expondrán las curvas de nivel de la temperatura mediante las gráficas de la Figura 2.

El gradiente de un campo vectorial es el siguiente:
[math]\nabla T(x,y) =\frac{d∂}{dx} +\frac{d∂}{dy}[/math]

en este caso sería:
[math]\nabla T(x,y) =\frac{6(x+1)}{1+(x+1)^2} +\frac{2(y-2)}{1+(y-2)^2}[/math]

A diferencia del código de MatLab anterior, en este se hace uso de la ecuación de la temperatura y de los comandos [math]quiver[/math] y [math]contour.[/math]
La temperatura máxima alcanzada es de [math]9,44343°C[/math] y se alcanza en el punto [math](1,12)[/math] 

clc,clear,clf
x1=-1:0.2:1;
y1=0:0.2:12;
[x2,y2]=meshgrid(x1,y1);
T= 3.*log(1+(x2+1).^2)+log(1+(y2-2).^2);
axis([-1,1,0,12]);
hold on
U=(6.*(x2+1)./(1+(x2+1).^2));
V=(2.*(y2-2)./(1+(y2-2).^2));
quiver(x1,y1,U,V);
xlabel('Eje X')
ylabel('Eje Y')
title('Gradiente')
contour(x2,y2,T,25)
colorbar
hold off


4 Ley de Fourier

Figura 3: Campo vectorial [math]\vec{Q}[/math]
De acuerdo a la Ley de Fourier, la energía calorífica [math]\vec{Q}[/math] viaja de acuerdo a la fórmula:
[math]\vec{Q} = −κ∇T[/math]

donde [math]κ[/math] es la constante de conductividad térmica de la placa que supondremos [math]κ = 1[/math].
Al ser [math]k=1[/math] se representará la energía calorífica como el negativo del gradiente.
Para representarlo, se tiene el gradiente calculado anteriormente, se procederá al uso de la parte negativa de las dos componentes y con el comando [math]quiver[/math] se representrá lo pedido.
A continuación se calculará [math]\vec{Q}[/math] y se dibujará como campo vectorial. Para ello, se hará uso del siguiente código: 

clc,clear,clf
x1=-1:0.2:1;
y1=0:0.2:12;
[x2,y2]=meshgrid(x1,y1);
T= 3.*log(1.+(x2+1).^2)+log(1.+(y2-2).^2);
hold on
U=(6.*(x2+1)./(1+(x2+1).^2));
V=(2.*(y2-2)./(1+(y2-2).^2));
QU=-U;
QV=-V;
quiver(x1,y1,QU,QV);
axis equal
xlabel('Eje X')
ylabel('Eje Y')
title('Ley de Fourier')
hold off


5 Campo de vectores

Figura 4: Campo de vectores

Ahora se mostrará el campo de vectores en los puntos del mallado del sólido.

Para calcular el campo de vectores, se calculará primero el vector
[math]\vec{u}(x, y, t) = \vec{a} sin(πk(\vec{d} · \vec{r_{0}}(x, y) − vt)),[/math]
siendo
[math]\vec{a} = \vec{d} = 1/3 \vec{j}, k = 1, t = 0[/math]

Con el vector calculado, solo es interesante [math]\vec{j}[/math], ya que [math]\vec{i}[/math] se anulará con el producto escalar. 

clc,clear,clf
x1=-1:0.2:1;
y1=0:0.2:12;
[x2,y2]=meshgrid(x1,y1);
r0x= x2;
r0y= y2;
ux=sin(0*r0x);
uy=(1/3).*sin((1/3).*r0y.*pi);
quiver(x1,y1,ux,uy)
xlabel('Eje X')
ylabel('Eje Y')
title('Campo de vectores')
axis equal
axis([-5,5,-1,13]);


6 Sólido antes y después del desplazamiento

Ahora se mostrará el sólido tanto antes como después del desplazamiento dado por el campo de vectores [math]u[/math]
Se utilizará el comando [math]subplot[/math] en el código empleado.

Figura 5: Antes y después del desplazamiento
clc,clear,clf
x1=-1:0.2:1;
y1=0:0.2:12;
[x2,y2]=meshgrid(x1,y1);
r0x= x2;
r0y= y2;
ux=sin(0*r0x);
uy=(1/3).*sin((1/3).*r0y.*pi);
hold on
subplot(1,2,1)
mesh(r0x,r0y,r0x*0)
title('Situación inicial');
xlabel('Eje X');
ylabel('Eje Y');
axis equal
axis([-5,5,-1,13]);
view(2);
subplot(1,2,2)
mesh(r0x+ux,r0y+uy,r0x*0)
title('Situación final');
xlabel('Eje X');
ylabel('Eje Y');
axis equal
axis([-5,5,-1,13]);
view(2);


7 Divergencia del campo vectorial

La divergencia es una medida del cambio de volumen local debido al desplazamiento, siendo la divergencia

[math]\nabla\cdot\vec{u}=\frac{∂}{∂x_1}(x_1) + \frac{∂}{∂x_2}(x_2) + \frac{∂}{∂x_3}(x_3)[/math]

El vector de posición viene dado por

[math]\vec{r}(x_1,x_2,x_3)= x_1\vec{i} + x_2\vec{j} + x_3\vec{k}[/math]

y en este caso (cuando [math]t=0[/math]) es

[math]\vec{r}(x_1,x_2,x_3)= \frac{1}{3}sen(\pi\cdot\frac{y}{3})\vec{j} [/math]

así que la divergencia es

[math]\nabla\cdot\vec{u}= \frac{\pi}{9}cos(\frac{\pi}{3}y)[/math]

Figura 6: Divergencia

Siendo [math]N[/math] un número natural, los puntos en los que la divergencia de [math]\vec{u}[/math] es máxima, mínima y nula son los siguientes:

  • Máxima: Cuando [math]y= 0+2\pi N[/math]
  • Mínima: Cuando [math]y=3+2\pi N[/math]
  • Nula: Cuando [math]y=\frac{3}{\pi+2\pi N}[/math]

El valor máximo es: [math]0,3491[/math] y el mínimo es: [math]-0,3491[/math] 


clc,clear,clf
x1=-1:0.2:1;
y1=0:0.2:12;
[x2,y2]=meshgrid(x1,y1);
r0x= x2;
r0y= y2;
ux=sin(0*r0x);
uy=(1/3).*sin((1/3).*r0y.*pi);
divu=(pi/9).*cos(pi.*r0y/3);%en t=0
pcolor(x2,y2,divu)
axis equal
axis([-5,5,-1,13]);
colorbar;
divuMAX=max(max(divu))
divuMIN=min(min(divu))


8 Rotacional del sólido

A continuación se calculara y dibujará el rotacional del campo vectorial [math]\vec{u}[/math] cuando [math]t=0[/math]

[math]|\nabla\times\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}\\ 0& \frac{1}{3}sen(\pi\cdot\frac{1}{3}) & 0 \end{vmatrix} = 0[/math]

El rotacional define como rota (gira) una pieza al someterla a un campo vectorial. En este caso, como es nulo, nuestra pieza no se ve afectada en ningún punto.

9 Tensor de tensiones

Definamos
[math]Ԑ(\vec{u}) = \frac{(∇\vec{u} + ∇\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} I + 2 μԐ[/math]

donde [math]λ[/math] y [math]μ[/math] son los conocidos como coeficientes de Lamé que dependen de las propiedades elásticas de cada material y [math]I[/math] el tensor identidad.
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], se dibujarán 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].
Lo primero, se calcula [math]∇\vec{u} + ∇\vec{u}^t[/math]

[math]∇\vec{u} = \begin{pmatrix} 0 & 0 & 0 \\ 0 & \frac{\partial u_{2} }{\partial y} & 0 \\ 0 & 0 & 0 \end{pmatrix} = ∇\vec{u}^t[/math]

Siendo
[math]\frac{\partial u_{2} }{\partial y}=\frac{\pi}{9}cos(\frac{\pi}{3}y)[/math]

entonces
[math]∇\vec{u} + ∇\vec{u}^t = \begin{pmatrix} 0 & 0 & 0 \\ 0 & \frac{2 \pi}{9}cos(\frac{\pi}{3}y) & 0 \\ 0 & 0 & 0 \end{pmatrix}[/math]

Lo siguiente a calcular es [math]\nabla\cdot\vec{u}[/math], para posteriormente multiplicarlo por el tensor identidad [math]I[/math]

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

Por tanto

[math] σ =\nabla\vec{u}\cdot I + Ԑ = \begin{pmatrix} \frac{\pi}{9}cos(\frac{\pi}{3}y) & 0 & 0 \\ 0 & \frac{\pi}{3}cos(\frac{\pi}{3}y) & 0 \\ 0 & 0 & \frac{\pi}{9}cos(\frac{\pi}{3}y) \end{pmatrix} [/math]

Figura 7: Tensiones normales en las direcciones [math]\vec{i},\vec{j}[/math] y [math]\vec{k}[/math]
clc,clear,clf
x1=-1:0.2:1;
y1=0:0.2:12;
[x2,y2]=meshgrid(x1,y1);
r0x= x2;
r0y= y2;
gradu=(pi/9).*cos(pi.*r0y/3);%en t=0
sigmai=(pi/9).*sin((1/3).*r0y.*pi);
sigmaj=(pi/3).*sin((1/3).*r0y.*pi);
sigmak=(pi/9).*sin((1/3).*r0y.*pi);
subplot(1,3,1)
surf(x2,y2,sigmai)
view(2)
xlabel('Eje X')
ylabel('Eje Y')
zlabel('Eje Z')
title('Tensiones normales i')
axis equal
caxis([0 1]);
axis([-5,5,-1,13]);
subplot(1,3,2)
surf(x2,y2,sigmaj)
view(2)
xlabel('Eje X')
ylabel('Eje Y')
zlabel('Eje Z')
title('Tensiones normales j')
axis equal
caxis([0 1]);
axis([-5,5,-1,13]);
subplot(1,3,3)
surf(x2,y2,sigmak)
view(2)
xlabel('Eje X')
ylabel('Eje Y')
zlabel('Eje Z')
title('Tensiones normales k')
axis equal
caxis([0 1]);
axis([-5,5,-1,13]);


10 Tensiones tangenciales

Ahora se calcularán las tensiones tangenciales respecto al plano ortogonal a [math]\vec{i}[/math], es decir [math]|σ·\vec{i} − (\vec{i}·σ·\vec{i})\vec{i}|[/math].

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

[math] σ\cdot\vec{i}= \frac{\pi}{9}cos(\frac{\pi}{3}y)[/math]

[math] (\vec{i}\cdot σ\cdot\vec{i})\cdot \vec{i}= \frac{\pi}{9}cos(\frac{\pi}{3}y)[/math]

Por tanto

[math] |σ·\vec{i} − (\vec{i}·σ·\vec{i})\cdot\vec{i}|= |\frac{\pi}{9}cos(\frac{\pi}{3}y)-\frac{\pi}{9}cos(\frac{\pi}{3}y)|=0 [/math]

En este caso, las tensiones tangenciales respecto al plano ortogonal a [math]\vec{i}[/math] son nulas.

11 Tensión de Von Mises

Figura 8: Tensión de Von Mises

La tensión de Von Mises se define por la fórmula

[math]σ_{VM}=\sqrt{\frac{(σ_{1}-σ_{2})^2+(σ_{2}-σ_{3})^2+(σ_{3}-σ_{1})^2}{2}}[/math]
donde [math]σ_{1}[/math], [math]σ_{2}[/math] y [math]σ_{3}[/math] son los autovalores de [math]σ[/math] (también conocidos como tensiones principales). Se trata de una magnitud escalar que se suele usar como indicador para saber cuando un material inicia un comportamiento plástico.

Gracias a ir analizando el valor en cada punto y almacenando el máximo (dentro del bucle), se puede obtener el punto de valor máximo, este es: [math](0.743590,10.461538,0.854340)[/math]
A diferencia del resto de códigos, en este se hará uso del comando [math]eig.m[/math] 

clc,clear,clf
x1=linspace(-1,1,40);
y1=linspace(0,12,40);
[x2,y2]=meshgrid(x1,y1);
fi=@(x,y) (pi/9).*sin((1/3).*y.*pi);
fj= @(x,y) (pi/3).*sin((1/3).*y.*pi);
fk=@ (x,y) (pi/9).*sin((1/3).*y.*pi);
n=length(x1); m=length(y1);
max=0;
for i=1:n
    for j=1:m
        xn=fi (x2(i,j),y2(i,j));
        yn=fj (x2(i,j),y2(i,j));
        zn=fk (x2(i,j),y2(i,j));
        Mn=[xn,0,0;0,yn,0;0,0,zn];
        autovalores=eig(Mn);
        VM=sqrt(((autovalores(1)-autovalores(2))^2+((autovalores(2)-autovalores(3))^2+((autovalores(3)-autovalores(1))^2))*1/2));
        if VM>max
            max=VM;
            x=x1(i);
            y=y1(i);
            z=VM;
        end
        Z(i,j)=VM;
    end
end
hold on
mesh(x2,y2,Z)
view(3)
shading interp
xlabel('Eje X')
ylabel('Eje Y')
zlabel('Eje Z')
title('Tensión de Von Mises')
axis equal
hold off
fprintf('El punto máximo es (%f,%f,%f)',x,y,z)


12 Campo de fuerzas que actúan sobre la placa

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

[math] \vec F=\frac{\partial^2 \vec u}{\partial t^2}-(\nabla \cdot σ) [/math]

donde [math]\nabla \cdot σ[/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]
A continuación se calculara la velocidad de propagación de las ondas [math]\vec{v}[/math] en términos de las constantes de Lamé, suponiendo que [math]\vec{F}=0[/math]
Se empieza calculando [math]\frac{\partial^2 \vec u}{\partial t^2}[/math], sabiendo que [math]\vec{u}= \frac{1}{3}sen(\frac{\pi y}{3}-vt)\vec{j}[/math]

[math]\frac{\partial \vec u}{\partial t} = \frac{-v}{3}cos(\frac{\pi y}{3}-vt)\vec{j}[/math]

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

[math]σ = λ∇ ·\vec{u} I + 2 μԐ=∇ ·\vec{u} I + (∇\vec{u} + ∇\vec{u}^t) [/math]

En este caso, se empezará primero con [math] ∇\vec{u} + ∇\vec{u}^t [/math]

[math]∇\vec{u} = \begin{pmatrix} 0 & 0 & 0 \\ 0 & \frac{\partial u_{2} }{\partial y} & 0 \\ 0 & 0 & 0 \end{pmatrix} = ∇\vec{u}^t[/math]

[math]\frac{\partial \vec u}{\partial y} = \frac{1}{3}cos(\frac{\pi y}{3}-vt)\cdot\frac{\pi}{3}\vec{j}= \frac{\pi}{9}cos(\frac{\pi y}{3}-vt)[/math]

[math]∇\vec{u} + ∇\vec{u}^t = \begin{pmatrix} 0 & 0 & 0 \\ 0 & \frac{2 \pi}{9}cos(\frac{\pi y}{3}-vt) & 0 \\ 0 & 0 & 0 \end{pmatrix} [/math]

Lo siguiente, será [math]∇ ·\vec{u}\cdot I[/math]

[math]∇ ·\vec{u}= \frac{\pi}{9}cos(\frac{\pi y}{3}-vt)[/math]

[math]∇ ·\vec{u} \cdot I= \begin{pmatrix} \frac{\pi}{9}cos(\frac{\pi y}{3}-vt) & 0 & 0 \\ 0 & \frac{\pi}{9}cos(\frac{\pi y}{3}-vt) & 0 \\ 0 & 0 & \frac{\pi}{9}cos(\frac{\pi y}{3}-vt) \end{pmatrix}[/math]

Por tanto
[math]σ = \begin{pmatrix} \frac{\pi}{9}cos(\frac{\pi y}{3}-vt) & 0 & 0 \\ 0 & \frac{\pi}{3}cos(\frac{\pi y}{3}-vt) & 0 \\ 0 & 0 & \frac{\pi}{9}cos(\frac{\pi y}{3}-vt) \end{pmatrix} [/math]

Se calcula [math]∇ ·σ[/math], sabiendo que es hacer la divergencia de los vectores fila y se iguala en la expresión de [math]\vec{F}=\vec{0}[/math]

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

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

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

13 Módulo de desplazamiento transversal

Fijando ahora el punto [math](1/2, 1)[/math], se calculará el módulo del desplazamiento transversal (dirección [math]\vec{i}[/math]) a lo largo del tiempo en el intervalo [math]t∈ [0, 10].[/math]
Como esta onda es longitudinal, no se desplazará en otra dirección que no sea la de su amplitud, siendo en este caso [math] \vec{a}=\frac{1}{3}\vec{j}[/math]
Se mostrará que no se produce ningún desplazamiento:

[math] \vec{v}(x,y,t)=\frac{1}{3}\vec{j}(sen(\frac{\pi y}{3}-\frac{\pi}{\sqrt{3}}t\vec{j}))[/math]

Para calcular un desplazamiento en una dirección (en este caso [math]\vec{i}[/math]) basta con multiplicar el campo [math]\vec{u}[/math] por dicha dirección escalarmente

[math] \vec{v}(\frac{1}{2},1,t)\cdot \vec{i}= [\frac{1}{3}\vec{j}(sen(\frac{\pi y}{3}-\frac{\pi}{\sqrt{3}}t\vec{j}))]\cdot\vec{i}=\vec{0}[/math]

No se desplaza en la dirección transversal.
Pasaría lo mismo en la dirección [math]\vec{k}[/math], mientras que en la dirección [math]\vec{j}[/math] se produciría todo el desplazamiento.

Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Visualización_de_campos_escalares_y_vectoriales_en_elasticidad._(Grupo_28)&oldid=66854»