Grupo B9-T2: Deformaciones elásticas de una placa plana en 2D

De MateWiki
Saltar a: navegación, buscar
Trabajo realizado por estudiantes
Título Grupo B9-T2:

Deformaciones elásticas de una placa plana en 2D.

Asignatura Teoría de Campos
Curso 2021-22
Autores Álvaro Parrilla

Alfonso Esplá

Teresa Perera

Pablo Sotos

Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

Enunciado del Trabajo 2:
Consideramos una placa rectangular plana (en dimensión 2) que ocupa la región (x, y) ∈ [0, 10] × [−1, 1].

En ella vamos a suponer que tenemos definidas dos cantidades físicas:

La temperatura T(x, y), que viene dada por: [math]T(ρ, θ) = (ρ − 2)^2 cos(θ) [/math] ; y los desplazamientos [math]\vec u(x,y) [/math] producidos por la acción de una fuerza determinada.

De esta forma, si definimos [math]\vec r_{0}(x,y) [/math] el vector de posición de los puntos de la placa antes de la deformación, la posición del punto (x,y) de la placa después de la deformación viene dada por:

[math]\vec r(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 de los puntos de la misma dado por el vector de desplazamientos

[math] \vec u = \frac {y^2}{80} \vec i + xf(y)\vec j [/math]

donde f(y) es una función que no conocemos.


1 Representación del sólido 2D

Para la representación del sólido delimitado por (x, y) ∈ [0, 10] × [−1, 1], hemos dibujado un mallado como paso de muestreo de [math] h=\frac {1}{2} [/math] para las variables x e y.


Representación del sólido
%Mallado (x,y) = (0,10)x(-1,1)
% Región del sólido
h=1/10;
x=[0:h:10];
y=[-1:h:1];
z=[0:0];
% Matriz de la región
[Mx,My,Mz]= meshgrid(x,y,z);
% Mallado
figure(1)
mesh(Mx,My,Mz);
% Gráfico ajustado
axis([-0.5 10.5 -1.5 1.5]);
view(0,90);
title('Sólido 2D')
xlabel('Eje X')
ylabel('Eje Y')


2 Líneas de nivel de la temperatura definida como un campo escalar y temperatura máxima

La distribución de temperatura en el sólido viene dada por una función escalar en coordenadas cilíndricas (polares ya que estamos en 2D): [math] T(ρ, θ) = (ρ − 2)^2 cos(θ) [/math]


Está distribución se puede resolver de dos formas:


Forma 1: Para que nos sea más fácil cálcularla y aplicarla en el sólido, convertiremos las coordenadas cilíndricas en cartesianas. [math] T(x,y) = x \sqrt{x^2+y^2}\ - 4x + \frac {4x}{\sqrt{x^2+y^2}} [/math]


Una vez convertida, ya podríamos dibujarlo en la gráfica.


Curvas de nivel de la temperatura (forma 1 y 2)
% Temperatura en cartesianas
% Definimos parámetros y creamos las matrices
h=1/10;
x=0:h:10;
y=-1:h:1;
[Mx,My]=meshgrid(x,y);
% Definimos la función de temperatura
T= Mx.* sqrt(Mx.^2+My.^2)-4*Mx+4*(Mx./sqrt(Mx.^2+My.^2));
% Dibujamos las líneas de nivel
figure(1)
contour(Mx,My,T)
axis([-0.5 10.5 -1.5 1.5])
xlabel('Eje X')
ylabel('EjeY')
colorbar
title('Curvas de nivel de la temperatura')


Forma 2: La segunda forma de resolverlo es un poco más tediosa, pero se obtiene lo mismo. Primero parametrizamos la región de interés en cartesianas, luego convertimos estas coordenadas en cilíndricas, aplicamos la función escalar de temperatura: [math] T(ρ, θ) = (ρ − 2)^2 cos(θ) [/math] , y por último dibujamos las curvas de nivel "contour(Mx,My,T)" 


%Curvas de nivel de la temperatura en cilíndricas
% Creo los vectores
h= 1/10;
x= 0:h:10;
y= -1:h:1;
%Creo las matrices
[Mx,My]=meshgrid(x,y);
%Convierto cada matriz a coordenadas cilíndricas
Mp = sqrt(Mx.^2+My.^2);
Mtt= atan(My./Mx);
%Aplico la fórmula de la temperatura
T=(Mp-2).^2.*cos(Mtt);
%Creo las curvas de nivel
figure(1)
contour(Mx,My,T);
axis([-0.5 10.5 -1.5 1.5])
xlabel('Eje X')
ylabel('EjeY')
colorbar
title('Curvas de nivel de la temperatura')


3 Cálculo y representación del Gradiente de la temperatura

En primer lugar, obtendremos la función temperatura en coordenadas cartesianas, proceso que ya se llevó a cabo en el ejercicio anterior. Posteriormente se realiza el gradiente de dicha función.


[math]\nabla T =\frac{\partial T}{\partial x}\vec i + \frac{\partial T}{\partial y} \vec j = (\sqrt{x^2+y^2}+\frac{4y^2}{(x^2+y^2)\sqrt{x^2+y^2}}-4+\frac{x^2}{\sqrt{x^2+y^2}})\vec i + (\frac{x^3y+xy^3-4xy}{{(x^2+y^2)}\sqrt{x^2+y^2}}) \vec j [/math]


Como se puede observar en el dibujo, el gradiente es ortogonal a las curvas de nivel del campo T, y por tanto, ortogonal a la superficie del campo T en todos sus puntos.


Gradiente de temperatura
%Gradiente de la temperatura en coordenadas cartesianas
%Parámetros y matrices
h=1/10;
x=0:h:10;
y=-1:h:1;
[Mx,My]=meshgrid(x,y);
figure(1)
% Definimos la función de temperatura
T= Mx.* sqrt(Mx.^2+My.^2)-4*Mx+4*(Mx./sqrt(Mx.^2+My.^2));
% Dibujamos las líneas de nivel
 contour(Mx,My,T)
 hold on
 xlabel('Eje X')
 ylabel('EjeY')
 axis([-0.5,10.5,-1.5,1.5])
 title('Gradiente de la temperatura')
 view(2)
 colorbar
%Componentes del campo vectorial
U=(sqrt(Mx.^2+My.^2))+(Mx.^2./(sqrt(Mx.^2+My.^2)).*(Mx.^2+My.^2))+(4.*My.^2./(sqrt(Mx.^2+My.^2)).*(Mx.^2+My.^2))-4;
V=((My.*Mx.^3)+(Mx.*My.^3)-(4.*Mx.*My))./((sqrt(Mx.^2+My.^2)).*(Mx.^2+My.^2));
%Representación del campo
quiver(Mx,My,U,V)
view(2)
hold off


4 Campo de desplazamientos

Sea [math]\vec u(x,y) = \frac{y^2}{80} \vec i + x f(y) \vec j [/math] calcularemos su rotacional, y aplicando las condiciones proporcionadas por el enunciado, obtendremos f(y).


[math]\nabla \times \vec u = \frac{\partial x f(y)}{\partial x} - \frac{\partial \frac{y^2}{80}}{\partial y} = f(y) -\frac{y}{40} = 0 [/math]


Por tanto, f(y)= [math] \frac{y}{10} [/math]


Una vez calculado f(y), podremos sustituirlo en la expresión del campo de desplazamiento, obteniendo que [math]\vec u(x,y) = \frac{y^2}{80} \vec i + \frac{xy}{40} \vec j [/math]


5 Representación del campo vectorial en los puntos del sólido

A continuación se representa mediante un gráfico realizado en Octave, el campo de vectores en los puntos del mallado del sólido.


Campo de desplazamientos
u=0:0.1:10;
v=-1:0.1:1;
[Mx,My]=meshgrid(u,v);
%Definimos el vector desplazamiento  
Ux= (My.^2)/80;                
Uy= Mx.*My /40;
%Representamos el campo vectorial
quiver(Mx,My,Ux,Uy);           
axis([-0.5,10.5,-1.5,1.5]);
axis equal
xlabel('Eje X')
ylabel('Eje Y')
title ('Campo de desplazamientos')


6 Situación inicial y final del sólido

En este apartado, procederemos a mostrar gráficamente el antes y el después del desplazamiento de la placa, junto con una superposición de las imágenes con la función plot3 para que se pueda apreciar visualmente la deformación


centro 

u=0:0.1:10;
v=-1:0.1:1;
[Mx,My]=meshgrid(u,v);
Ux= My.^2/80 ;
Uy= Mx.*My /40 ;
figure(1)
%ANTES
subplot(1,3,1)
mesh(Mx,My,0*Mx);
axis([-0.5,10.5,-1.5,1.5]);
axis equal
title('Antes del desplazamiento');
view(2)
xlabel('Eje X')
ylabel('Eje Y') 
%DESPUÉS
subplot(1,3,2)
mesh(Ux+Mx ,Uy+My,0*Ux);
axis([-0.5,10.5,-1.5,1.5]);
axis equal 
title('Despues del desplazamiento')
view(2)
xlabel('Eje X')
ylabel('Eje Y') 
%COMPARACIÓN
subplot(1,3,3)
plot3(Mx,My,Mx*0);
hold on 
plot3(Ux+Mx,Uy+My,0*Ux);
hold off
view(2)
axis([-0.5,10.5,-1.5,1.5]);
axis equal
title ('Comparacion')
xlabel('Eje X')
ylabel('Eje Y')


7 Divergencia máxima, mínima y nula de [math]\vec u [/math]

Ahora procedemos a dibujar la divergencia del vector [math]\vec u [/math], [math]\nabla · \vec u [/math] , y determinar analíticamente los puntos en los que es máxima, mínima y nula.


Procedemos a calcularla analíticamente: dado el vector [math]\vec u(x,y) = \frac{y^2}{80} \vec i + \frac{xy}{40}\vec j [/math] , la divergencia: [math]\nabla ·\vec u = \frac{\partial u_1}{\partial x} + \frac{\partial u_2}{\partial y} = \frac{x}{40}. [/math]


centro 

u=0:0.1:10;
v=-1:0.1:1;
[Mx,My]=meshgrid(u,v);
%Representamos la divergencia
div= Mx./40;           
surf(Mx,My, div);
axis ([-0.5,10.5,-1.5,1.5])
axis equal
view(2)
xlabel('X')
ylabel('Y')
colorbar
title ('Divergencia')


La divergencia es máxima en x=10, y su valor es 0,25. Por otro lado, será mínima (y nula) en x=0


8 Módulo del rotacional de [math]\vec u [/math]

Calculamos el |∇× [math]\vec u [/math] | en todos los puntos del sólido.


Sabiendo que:


[math] f(y) =\frac{y}{40} [/math]


[math] \vec{u} (x,y) = \frac {y^2}{80} \vec i + x f(y) \vec j = \frac {y^2}{80}\vec i + \frac{xy}{40} \vec j [/math]


[math] \nabla \times \vec u =\left| \begin{matrix}\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{matrix}\right| =\left| \begin{matrix} \vec i & \vec j & \vec k \\ & & \\ \frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \\ & & \\ \frac{y^2}{80} & \frac{xy}{40} & 0 \end{matrix}\right| =\left( \frac{\partial 0}{\partial y} - \frac{ \partial \frac{xy}{40}}{\partial z}\right)\vec i + \left(\frac{\partial \frac{y^2}{80}}{\partial z} - \frac{\partial 0}{\partial x}\right)\vec j + \left(\frac{ \partial \frac{xy}{40}}{\partial x} - \frac{ \partial \frac{y^2}{80}}{\partial y}\right)\vec k =\left(0 - 0\right)\vec i + \left(0 - 0\right)\vec j + \left(\frac{y}{40}-\frac{2y}{80}\right)\vec k =0 [/math]


El rotacional es nulo, por lo que su módulo también lo es y no existe rotación en la placa.


9 Cálculo del tensor de tensiones y representación de este en la base física cartesiana

La fórmula que aplicaremos para el cálculo del tensor de tensiones es:


[math]\sigma = λ \nabla · \vec u I + 2µ \epsilon [/math]


Donde λ , μ son los coeficiente de Lamé, los cuales dependen de las cualidades elásticas del material. En este caso ambos coeficientes valen 1.


Primero, para calcular la [math]\epsilon[/math], obtendremos el gradiente de [math]\vec u [/math] y su traspuesto:


[math] \vec{u} (x,y) =\frac {y^2}{80}\vec i + \frac{xy}{40} \vec j [/math]


[math]\nabla \vec u =\left ( \begin{matrix} \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{matrix} \right ) = \left ( \begin{matrix} 0 & \frac{y}{40} & 0 \\ & \\ \frac{y}{40} & \frac{x}{40} & 0 \\ & \\ 0 & 0 & 0 \end{matrix} \right ) [/math]


Al ser una matriz simétrica, su traspuesto es la misma matriz:


[math]\nabla \vec u ^t= \left ( \begin{matrix} 0 & \frac{y}{40} & 0 \\ & \\ \frac{y}{40} & \frac{x}{40} & 0 \\ & \\ 0 & 0 & 0 \end{matrix} \right ) [/math]


Por lo que la matriz simétrica es ella misma, aplicando la fórmula:


[math]\epsilon(\vec u) = \frac{1}{2} (∇\vec u + ∇\vec u ^t)[/math]


Obtenemos:


[math]\epsilon(\vec u)= \left ( \begin{matrix} 0 & \frac{y}{40} & 0 \\ & \\ \frac{y}{40} & \frac{x}{40} & 0 \\ & \\ 0 & 0 & 0 \end{matrix} \right ) [/math]


Con esto ya tendríamos todos los datos necesarios (la [math] \nabla · \vec u [/math] ya la hemos obtenido en apartados anteriores), solo nos falta sustituir en la fórmula:


[math]\sigma_{ij} = \left( \begin{matrix} \frac{x}{40} & 0 & 0 \\ & \\ 0 & \frac{x}{40} & 0 \\ & \\ 0 & 0 & \frac{x}{40} \end{matrix} \right ) + 2 · \left ( \begin{matrix} 0 & \frac{y}{40} & 0 \\ & \\ \frac{y}{40} & \frac{x}{40} & 0 \\ & \\ 0 & 0 & 0 \end{matrix} \right ) = \left ( \begin{matrix} \frac{x}{40} & \frac{y}{20} & 0 \\ & \\ \frac{y}{20} & \frac{3x}{40} & 0 \\ & \\ 0 & 0 & \frac{x}{40} \end{matrix} \right ) [/math]


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

Aplicamos: [math]\vec i \cdot \sigma_{ij} \cdot \vec i = \frac{x}{40}[/math]


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

Aplicamos: [math]\vec j \cdot \sigma_{ij} \cdot \vec j = \frac{3x}{40}[/math]


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

Aplicamos: [math]\vec k \cdot \sigma_{ij} \cdot \vec k = \frac{x}{40}[/math]


Tensiones normales direccionales
h=1/10;
x=0:h:10;
 y=-1:h:1;
 [xx,yy]=meshgrid(x,y);
 %Tensiones normales en las direcciones de los ejes coordenados 
 tni=xx./40;
 tnj=(3.*xx)./40;
 tnk=xx./40;
  %Representación tensiones normales
  figure
 subplot(1,3,1)
 surf(xx,yy,tni)
 axis image
 view(2)
 colorbar
 xlabel('x')
 ylabel('y')
 zlabel('z')
 title('Tensión normal dirección eje i')
 subplot(1,3,2)
 surf(xx,yy,tnj)
 axis image
 view(2)
 colorbar
 xlabel('x')
 ylabel('y')
 zlabel('z')
 title('Tensión normal dirección eje j')
 subplot(1,3,3)
 surf(xx,yy,tnk)
 axis image
 view(2)
 colorbar
 xlabel('x')
 ylabel('y')
 zlabel('z')
 title('Tensión normal dirección eje k')


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

Para poder obtener las tensiones tangenciales al plano ortogonal a [math]\vec i [/math], aplicaremos la fórmula:


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


sustituyendo los datos del apartado anterior.


[math]\left | \sigma \cdot \vec i - (\vec i \cdot \sigma \cdot \vec i ) \vec i \right | = \left | \left ( \begin{matrix} \frac{x}{40} & \frac{y}{20} & 0 \\ & \\ \frac{y}{20} & \frac{3x}{40} & 0 \\ & \\ 0 & 0 & \frac{x}{40} \end{matrix} \right ) \left ( \begin {matrix} 1 \\ & \\ 0 \\ & \\ 0 \end{matrix} \right )- (\frac{x}{40}) \left ( \begin {matrix} 1 \\ & \\ 0 \\ & \\ 0 \end{matrix} \right ) \right | = \frac {y}{20} [/math]


11 Cálculo y representación de la tensión de Von Mises [math](σ_{VM})[/math]

En este apartado habremos de calcular la tensión de Von Mises, usada para saber cuando un material inicia un comportamiento plástico.


Su expresión es [math] \sigma_{MV} [/math]=[math]\sqrt{\frac{(σ_1-σ_2)^2+(σ_2-σ_3)^2+(σ_3-σ_1)^2}{2}}[/math] donde [math] \sigma_{1} [/math], [math] \sigma_{2} [/math] y [math] \sigma_{3} [/math] son los autovalores de σ, cuya matriz se calculó en el apartado anterior:


[math]\sigma_{ij} = \left( \begin{matrix} \frac{x}{40} & 0 & 0 \\ & \\ 0 & \frac{x}{40} & 0 \\ & \\ 0 & 0 & \frac{x}{40} \end{matrix} \right ) + 2 · \left ( \begin{matrix} 0 & \frac{y}{40} & 0 \\ & \\ \frac{y}{40} & \frac{x}{40} & 0 \\ & \\ 0 & 0 & 0 \end{matrix} \right ) = \left ( \begin{matrix} \frac{x}{40} & \frac{y}{20} & 0 \\ & \\ \frac{y}{20} & \frac{3x}{40} & 0 \\ & \\ 0 & 0 & \frac{x}{40} \end{matrix} \right ) [/math]


A continuación, se calcularan con Matlab dichos autovalores, obteniendo el valor y la representación de la tensión de Von Mises.


Tensión de Von Mises
%Definimos las variables y las matrices
h=1/10;
x=0:h:10;
y=-1:h:1;
[xx,yy]=meshgrid(x,y);
%Matriz tensión
A=inline('x./40','x','y');
B=inline('y./20','x','y');
C=inline('3.*x./40','x','y');
tension=[];
for i=1:length(y)
     for j=1:length(x)
         D=A(xx(i,j),yy(i,j));
         E=B(xx(i,j),yy(i,j));
         F=B(xx(i,j),yy(i,j));
         G=C(xx(i,j),yy(i,j));
         H=A(xx(i,j),yy(i,j));
         tension=[D E 0; F G 0; 0 0 H];       
%Se obtienen los autovalores con la función eig.m
         autoval=eig(tension);
         tensionVM=sqrt(((autoval(1)-autoval(2))^2+(autoval(2)-autoval(3))^2+(autoval(3)-autoval(1))^2)/2);
         TENSION(i,j)=tensionVM;
     end
end
%Gráfica de la tensión de Von Mises
figure
surf(xx,yy,TENSION)
title('Tensión de Von Mises')
view(2)
colorbar
xlabel('x')
ylabel('y')
%Cálculo del máximo valor de la tensión de Von Mises
M=TENSION(1,1:length(TENSION));
maxTENSION=max(M);
fprintf('El valor máximo de la tensión de Von Mises es %f\n',maxTENSION)


El valor máximo de la tensión en la placa será de 0.5074.


12 Cálculo y representación del campo de fuerzas [math]\vec F [/math] que actúa sobre la placa

El campo de fuerzas [math]\vec F [/math], es el causante del desplazamiento observado. Este se aproxima con la ecuación de elasticidad lineal:


[math]\vec {F}=-\nabla \cdot σ [/math]


La divergencia del campo tensorial "σ", con signo negativo, es la fuerza aplicada en la placa.


12.1 Representación del campo

Para su cálculo, una vez definido el campo tensorial: [math] \sigma_{i,j}[/math] = \begin{pmatrix} \frac{x}{40} & \frac{y}{20} & 0\\ \frac{y}{20} & \frac{3x}{40} & 0\\ 0 & 0 & \frac{x}{40} \end{pmatrix}


Para calcular la divergencia del campo, sumamos las derivadas de cada fila de la matriz de componentes del tensor σ, respecto a "x" (la primera columna), "y" (la segunda columna) y "z" (la tercera columna). Así obtendremos un campo vectorial teniendo como componentes: [math] F_1,F_2,F_3[/math] correspondientes al vector [math] \vec F[/math]


[math] F_1 = - \frac{\partial σ_{j1}}{\partial x} = -\frac{3}{40}[/math] --> correspondiente a [math]\vec i[/math]


[math] F_2 = - \frac{\partial σ_{j2}}{\partial y} = 0[/math] --> correspondiente a [math]\vec j[/math]


[math] F_3 = - \frac{\partial σ_{j3}}{\partial z} = 0[/math] --> correspondiente a [math]\vec k[/math]


Finalmente obtenemos:


[math] \vec F = -\frac{3}{40} \vec i [/math]


12.2 Interpretación de la gráfica

Al dibujar la gráfica de este campo vectorial, observamos que la fuerza que actúa para deformar la placa tiene la dirección, pero el sentido contrario al vector [math]\vec i [/math]. Esta fuerza hace que el extremo izquierdo de la placa no se vea perjudicado (x=0), pero a medida que te hacercas al extremo x=10, existe una deformación progresiva del sólido, siendo máxima en este extremo.


Representación del campo de fuerzas [math]\vec F [/math]
%Representamos el mallado
h=1/10;
x=0:h:10;
y=-1:h:1;
[Mx,My]=meshgrid(x,y);
% Componentes del campo vectorial F
Fx= -3/40+0*Mx;
Fy= 0*My;
%Representamos el campo de Fuerzas F
figure(1)
quiver(Mx,My,Fx,Fy)
axis([-0.5 10.5 -1.5 1.5])
xlabel('Eje X')
ylabel('Eje Y')
title('Representación del campo de fuerzas')
view(0,90)
Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Grupo_B9-T2:_Deformaciones_elásticas_de_una_placa_plana_en_2D&oldid=47666»