Modelización del comportamiento de una placa bidimensional

De MateWiki
Saltar a: navegación, buscar
Trabajo realizado por estudiantes
Título Modelización del comportamiento de una placa bidimensional. Grupo 1-C
Asignatura Teoría de Campos
Curso 2022-23
Autores Mario Andres Silva Peñaloza, Matias Rodriguez Obon, Javier Nogales Sanchez, Sofia Benito Gomez, Miguel Angel Abad Robles
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura


En este trabajo empezaremos considerando una placa rectangular plana (bidimensional). Ésta ocupa la región definida entre 0 y 10 en el eje X y entre 0 y 2 en el eje Y, con su región en Z manteniéndose indefinida dado que sus cualidades físicas se realizan independientemente de la sección transversal de la placa. De antemano, se definirán dos cualidades físicas pertenecientes a esta placa, siendo su temperatura (T) y sus desplazamientos producidos por una fuerza determinada ([math]\vec u [/math]).

La temperatura viene dada por la siguiente función que depende de la posicion en la que se encuentra cada punto de la placa respecto del eje X (x) y del eje Y (y).

[math] T(x, y) = (x − 3)^2 + (10(y-1/2))^2 [/math]

Los desplazamientos ([math]\vec u [/math]) vendrán definidos por los vectores [math]\vec a, \vec d, \vec r_0(x,y) [/math] , que son la amplitud, el numero de onda, y un vector unitario que marca la dirección de propagación respectivamente. También dependerá de una serie de escalares, que son la velocidad de propagación v y el tiempo transcurrido t. Para los primeros 10 apartados consideraremos el tiempo t=0, y tomaremos en cuenta que [math]\vec a=\frac{2}{5}\vec j[/math], k=1, y [math]\vec d=\vec i[/math].

[math] \vec u=\vec a sen(k(\vec d ·\vec r_0(x,y)-vt)) [/math]

siendo

[math] \vec r_0(x,y)= x\vec i + y\vec j [/math]


1 Representacion del solido en 2D

Para representar el sólido, que es una placa rectangular plana con dimensiones (x, y) ∈ [0, 10] × [0, 2], se han hecho dos vectores x e y, con paso [math] h=\frac {2}{10} [/math], después se crea una cuadricula 2D con meshgrid y con el comando mesh se grafica el mallado en 2D.

Representación del sólido
h=2/10;
x= 0:h:10;
y = 0:h:2;
[XX,YY]=meshgrid(x,y);
mesh(XX,YY,0*XX);
view(2)
title("Mallado")
axis equal
axis([-0.5,10.5,-0.5,2.5]);
hold on
xlabel('Eje X')
ylabel('Eje Y')
plot([0,10],[2,2],'k',[10,10],[2,0],'k',[0,10],[0,0],'k',[0,0],[2,0],'k');
hold off


2 Curvas de nivel de la temperatura

La temperatura viene dada por la función [math] T(x, y) = (x − 3)^2 + (10(y-1/2))^2 [/math] , para obtener las curvas de nivel usamos el comando contour

Curvas de nivel
clear all
clc
h=2/10;
x= 0:h:10;
y = 0:h:2;
[XX,YY]=meshgrid(x,y);mesh(XX,YY,0*XX);
z=(XX-3).^2+(10*(YY-1/2)).^2;
contour(XX,YY,z)
title("Curvas de nivel de la temperatura")
axis([-0.5,10.5,-0.5,2.5])
hold on
colorbar
plot([0,10],[2,2],'k',[10,10],[2,0],'k',[0,10],[0,0],'k',[0,0],[2,0],'k');
axis equal
axis([-0.5,10.5,-0.5,2.5]);
xlabel('Eje X')
ylabel('Eje Y')
%valor maximo
[m,n]=size(z)
maximus=0
for a=1:1:m
    for b=1:1:n
        if z(a,b)>maximus
            maximus=z(a,b);
            ex=a;
            ye=b;
        else
        end
    end
end
maximus
fprintf('el maximo de temperatura es %f y se encuentra en las coordenadas x=%d; y=%d \n',maximus,XX(ex,ye),YY(ex,ye))
hold off

Con el bucle obtenemos que la temperatura máxima es 274 y se encuentra en las coordenadas X = 10 e Y = 2

3 Gradiente de la temperatura

Para dibujar el campo vectorial, obtenemos el gradiente de nuestra función de temperatura.

[math]\nabla T =\frac{\partial T}{\partial x}\vec i + \frac{\partial T}{\partial y} \vec j =(2x-6)\vec i + (20y-10) \vec j [/math]

Usamos la grafica del apartado anterior y graficamos el campo vectorial con quiver, se puede ver que los vectores son ortogonales a las curvas de nivel.

Curvas de nivel y gradiente
h=2/10;
x= 0:h:10;
y = 0:h:2;
[XX,YY]=meshgrid(x,y);
mesh(XX,YY,0*XX);
i= (2.*XX-6)
j=(200.*YY-100)
z=(XX-3).^2+(10*(YY-0.5)).^2
quiver(XX,YY,i,j)
hold on
title("Curvas de nivel y gradiente")
contour(XX,YY,z,7)
axis equal
axis([-0.5,10.5,-0.5,2.5]);
xlabel('Eje X')
ylabel('Eje Y')
hold off


4 Campo de vectores del desplazamiento

Se puede visualizar el desplazamiento que van a tener las partículas de la placa al ser sometida a los desplazamientos descritos por el vector [math] \vec u [/math], al ser t = 0,[math] \vec a =\frac {2}{5} \vec j, k = 1 , \vec d = \vec i [/math]. El vector [math] \vec u [/math] seria [math] \vec u=\frac {2}{5}sen(x)\vec j [/math] generando una gráfica utilizando Matlab con el siguiente código.

Campo de vectores del desplazamiento
h=2/10;
x= 0:h:10;
y = 0:h:2;
[XX,YY]=meshgrid(x,y);
ui=XX*0 
uj=2/5*sin(XX) 
quiver(XX,YY,ui,uj)
view(2)
title("Campo de vectores del desplazamiento")
axis equal


Como se puede ver en la gráfica resultante, el campo vectorial describe un movimiento ondulatorio. Esto llevará a la creación de curvas en lo que antes era una placa completamente planas. Sin embargo, se puede ver que el campo incrementa y disminuye gradualmente, lo cual demuestra que las curvas serán fluidas; no serán pliegues bruscos y marcados.

5 Solido antes y después del desplazamiento

Al ser sometido a una fuerza que provoca un desplazamiento la placa es deformada. El vector [math]\vec u [/math] que describe el desplazamiento de los vectores es sinusoidal, genera el campo vectorial visto en el apartado anterior y por tanto es predecible que provoque que la placa se curve. Podemos observar la posición de los puntos de la placa al generar dos gráficas utilizando Matlab para ver el antes y el después.

Solido antes y después del desplazamiento
clear all
clc
%problema 5
h=2/10;
x= 0:h:10;
y = 0:h:2;
[XX,YY]=meshgrid(x,y);
figure(5)
subplot(2,1,1)
mesh(XX,YY,0*XX);
hold on
view(2)
title("Solido antes del desplazamiento")
axis([-0.5,10.5,-0.5,2.5]);
axis equal
xlabel('Eje X')
ylabel('Eje Y')
plot([0,10],[2,2],'k',[10,10],[2,0],'k',[0,10],[0,0],'k',[0,0],[2,0],'k');
hold off
subplot(2,1,2)
ui=XX;
uj=YY+2/5*sin(XX);
mesh(XX,YY+uj,0*XX)
view(2)
title("Solido despues del desplazamiento")
axis([-0.5,10.5,-0.5,2.5]);
axis equal
xlabel('Eje X')
ylabel('Eje Y')


6 Divergencia de [math]\vec u [/math]

La placa puede sufrir deformaciones debido a la fuerza aplicada a ella, y también cambios de volúmenes. Éstos se pueden determinar analíticamente al calcular la divergencia del vector [math]\vec u [/math], [math]\nabla · \vec u [/math].

Dado el vector [math] \vec u=\frac {2}{5}sen(x)\vec j [/math] su divergencia sera [math]\nabla ·\vec u = \frac{\partial u_1}{\partial x} + \frac{\partial u_2}{\partial y} = 0. [/math]

Se puede observar que la divergencia es 0, lo cual demuestra que la placa no sufre cambios de volumen locales. Esto se debe a que estamos tratando con un solido rígido, con una masa y una densidad constantes, por lo que el volumen también lo será. Esto hace que ilustrar una gráfica mostrando el cambio de la gradiente en cada punto sea redundante, ya que generaría un plano nulo.

7 Rotacional de [math]\vec u [/math]

Calculando el rotacional del vector que describe el desplazamiento de las partículas que forman la placa proporciona más información respecto a los puntos que sufrirán mayores alteraciones físicas. Se hace de la siguiente manera

Sabiendo que [math] \vec u=\frac {2}{5}sen(x)\vec j [/math], |∇× [math]\vec u [/math] | seria igual a: [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} \\ & & \\ 0 & \frac{2}{5}sen(x) & 0 \end{matrix}\right| =\left( \frac{\partial 0}{\partial y} - \frac{ \partial 0}{\partial z}\right)\vec i + \left(\frac{\partial 0}{\partial z} - \frac{\partial 0}{\partial x}\right)\vec j + \left(\frac{ \partial \frac{2}{5}sen(x)}{\partial x} - \frac{ \partial 0}{\partial y}\right)\vec k =0\vec i + 0\vec j + \left(\frac{2}{5}cos(x)\right)\vec k =\frac{2}{5}cos(x)\vec k [/math]

|∇× [math]\vec u [/math] | = [math] \frac {2}{5}cos(x) [/math]


Este resultado puede ser interpretado más fácilmente generando una gráfica y un bucle que nos indicará el punto que sufre el mayor rotacional utilizando Matlab. El código utilizado es el siguiente.

Rotacional en 2-D y 3-D
h=2/10;
x= 0:h:10;
y = 0:h:2;
[XX,YY]=meshgrid(x,y);
rot=2/5*cos(XX)
subplot(1,2,1)
surf(XX,YY,rot)
view(2)
axis equal
axis([-0.5,10.5,-0.5,2.5]);
title("Rotacional 2-D en t=0")
colorbar
xlabel('Eje X')
ylabel('Eje Y')
subplot(1,2,2)
surf(XX,YY,rot)
view(3)
title("Rotacional 3-D en t=0")
colorbar
xlabel('Eje X')
ylabel('Eje Y')
zlabel('Eje Z')


8 Tensor de tensiones

Consideramos que la placa esta siendo tratada con un medio elástico lineal, isótopo, y homogéneo, por lo que podemos describir un tensor de tensores ([math]\sigma [/math]). λ y µ son los coeficientes de Lamé, escalares que dependen de las propiedades elásticas del material en cuestión. Por otra parte tenemos la divergencia de [math]\vec U [/math] y ϵ que representa la parte simétrica del tensor gradiente de [math]\vec U [/math]. Estos se operan de la siguiente manera para obtener el tensor de tensiones.

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

Particularizando para nuestro caso, podemos expandir la función y realizar los cálculos recordando que estamos considerando t=0, y tomando en cuenta que λ=µ=1.

[math]\sigma = 2\epsilon [/math], siendo [math]\epsilon(\vec u) = \frac{1}{2} (∇\vec u + ∇\vec u ^t)[/math], operando nos queda [math]\sigma = \epsilon [/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 & 0 & 0 \\ & \\ \frac{2}{5}cos(x) & 0 & 0 \\ & \\ 0 & 0 & 0 \end{matrix} \right ) [/math]

Su traspuesta seria:


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

Y el resultado de [math]\sigma = \epsilon [/math]


[math]\sigma_{ij} = \left ( \begin{matrix} 0 & \frac{2}{5}cos(x) & 0 \\ & \\ \frac{2}{5}cos(x) & 0 & 0 \\ & \\ 0 & 0 & 0 \end{matrix} \right ) [/math]

Este resultado puede ser aplicado para calcular las tensiones normales en la dirección de cada uno de los ejes ([math]\vec i [/math], [math]\vec j [/math], [math]\vec k [/math]). Se opera de la siguiente manera.

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

Usamos: [math]\vec i \cdot \sigma_{ij} \cdot \vec i = 0[/math]

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

Usamos: [math]\vec j \cdot \sigma_{ij} \cdot \vec j = 0[/math]

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

Usamos: [math]\vec k \cdot \sigma_{ij} \cdot \vec k = 0[/math]

Al resultar todas nulas, podemos deducir que no hay tensiones normales en ninguna dirección provocados por el desplazamiento de los puntos de la placa.

9 Tensión tangencial al plano ortogonal [math]\vec i [/math]

Aprovechando que hemos obtenido el tensor de tensores en el apartado anterior, podemos utilizarlo para calcular la tensión tangencial respecto al plano ortogonal a [math]\vec i [/math] tras realizar los siguientes cálculos.

La tensión tangencial respecto a [math]\vec i [/math] es [math]|\sigma \cdot \vec i-(\vec i \cdot \sigma \cdot \vec i)\vec i |[/math], como [math] (\vec i \cdot \sigma \cdot \vec i)=0[/math], nos quedaría solo [math]\sigma \cdot \vec i [/math]


[math]\sigma \cdot \vec i = \left ( \begin{matrix} 0 & \frac{2}{5}cos(x) & 0 \\ & \\ \frac{2}{5}cos(x) & 0 & 0 \\ & \\ 0 & 0 & 0 \end{matrix} \right ) \left ( \begin {matrix} 1 \\ & \\ 0 \\ & \\ 0 \end{matrix} \right ) = \frac{2}{5}cos(x) [/math]


Al generar un resultado no nulo podemos dibujarlo para facilitar su visualización e interpretación con el siguiente código en Matlab.

Tensiones tangenciales respecto al plano ortogonal i
h=2/10;
x= 0:h:10;
y = 0:h:2;
[XX,YY]=meshgrid(x,y);
figure(7)
tension= 2/5*cos(XX);
quiver(XX,YY,tension,tension*0)
view(2)
axis equal
axis([-0.5,10.5,-0.5,2.5]);
title("Tensiones tangenciales respecto al plano ortogonal i")
xlabel('Eje X')
ylabel('Eje Y')


Podemos observar que oscilan entre unos valores máximos y mínimos, cambiando con un frente de onda paralelo al Eje Y de manera gradual y progresiva. Esto era predecible por el hecho que la tensión tangencial respecto al plano ortogonal venia descrito por una función trigonométrica dependiente de x.

10 Tensión de Von Mises

La tension de Von Mises es una magnitud escalar que sirve como indicadora del tipo de comportamiento de un material, si es completamente elástico o si tiene un comportamiento plástico. Se calcula aplicando la siguiente formula que utiliza los autovalores de la matriz del tensor de tensores.

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

Empezamos calculando estos autovalores igualando el polinomio característico de [math] \vec \sigma [/math] a 0, y procedemos insertando estos autovalores en la formula. Se obtienen los siguientes resultados.


[math]\sigma_{VM} = \left ( \begin{matrix} -\lambda & \frac{2}{5}cos(x) & 0 \\ & \\ \frac{2}{5}cos(x) & -\lambda & 0 \\ & \\ 0 & 0 & -\lambda \end{matrix} \right ) = -\lambda (\lambda ^2 - \frac{4}{25}cos(x)^2) [/math] Obtenemos los siguientes autovalores [math] \lambda_{1} = 0, \lambda_{2} = \frac{2}{5}cos(x) [/math] y [math] \lambda_{3} = -\frac{2}{5}cos(x) [/math]

Desarrollando la formula de Von Mises nos queda

[math]\sigma_{VM} = \frac{2 \sqrt3}{5}cos(x) [/math]

Podemos generar una grafica utilizando Matlab que nos mostrará como se distribuye esta tensión a lo largo de la placa. Esto va a generar unos resultados visiblemente parecidos a los generados en el apartado anterior pero con un significado ligeramente distinto y con una dilatación dado que el coseno de x esta siendo multiplicado por una constante diferente.

Tension de Von Mises
clc
clear all
h=2/10;
x= 0:h:10;
y = 0:h:2;
[XX,YY]=meshgrid(x,y);
figure(7)
VM= 2*sqrt(3)/5*cos(XX);
surf(XX,YY,VM)
view(2)
axis equal
axis([-0.5,10.5,-0.5,2.5]);
title("Tension de Von Mises")
colorbar
xlabel('Eje X')
ylabel('Eje Y')

%valor maximo
[m,n]=size(VM)
maximus=0
for a=1:1:m
    for b=1:1:n
        if VM(a,b)>maximus
            maximus=VM(a,b);
            ex=a;
            ye=b;
        else
        end
    end
end
maximus
fprintf('el maximo es %f y se encuentra en las coordenadas x=%d; y=%d \n',maximus,XX(ex,ye),YY(ex,ye))


Se puede ver que la descripción que hicimos anteriormente era acertada, y se puede observar que el primer máximo se encuentra donde x es igual a 0 y se repite periódicamente.

11 Velocidad de propagación

En muchos momentos se ha hecho referencia a la fuerza aplicada a la placa que ha causado los desplazamientos descritos por [math]\vec u [/math]. Se puede aproximar al aplicar la siguiente ecuación e la elasticidad lineal

[math] \vec F=\frac{\partial^2 u}{\partial t^2}-\nabla · \sigma [/math]

Si dejamos de asumir que t=0 y que los coeficientes de Lamé no son igual a 1, podemos utilizarla para calcular la velocidad de propagación de las ondas en función de estos coeficientes. Para la realización de estos cálculos tenemos que tomar en cuenta que la definición de [math]\vec u [/math] cambia por lo que también cambiara nuestro tensor de tensores (se realiza el mismo procedimiento que en el octavo apartado con [math]\vec u [/math] dependiente de t). [math]\vec u [/math] y [math]\vec sigma [/math] nos quedan de la siguiente manera.

[math] \vec u=\frac {2}{5}sen(x-vt)\vec j [/math]

[math]\sigma= µ\left ( \begin{matrix} 0 & \frac{2}{5}cos(x-vt) & 0 \\ & \\ \frac{2}{5}cos(x-vt) & 0 & 0 \\ & \\ 0 & 0 & 0 \end{matrix} \right ) [/math]

Podemos comenzar a realizar los caulculos para obtener [math] \vec F [/math] utilizando nuestro nuevo [math] \vec u [/math] y [math] sigma [/math]. Estaremos definiendo [math]\nabla · sigma [/math] como el campo vectorial obtenido tras hacer la divergencia de las filas de la matriz [math] sigma [/math]. Tras operar obtenemos el siguiente resultado.

[math] \vec F=-v^2\frac {2}{5}sen(x-vt)\vec j-\frac {2}{5}µsen(x-vt)\vec j [/math]

Si igualamos [math] \vec F=\vec 0 [/math] podemos despejar v que nos queda en función de µ, siendo [math]v=\sqrt{µ}[/math].

Estos resultados serían diferentes si tomamos [math]\vec a=\frac{2}{5}\vec i[/math]. [math]\vec u , \sigma, y \vec F [/math] obtendrían los siguientes valores.

[math] \vec u=\frac {2}{5}sen(x-vt)\vec i[/math]

[math]\sigma= -λ\frac {2}{5}sen(x-vt)\vec j \left ( \begin{matrix} 1 & 0 & 0 \\ & \\ 0 & 1 & 0 \\ & \\ 0 & 0 & 1 \end{matrix} \right ) [/math] [math] +µ\left ( \begin{matrix} \frac{4}{5}cos(x-vt) & 0 & 0 \\ & \\ 0 & 0 & 0 \\ & \\ 0 & 0 & 0 \end{matrix} \right ) [/math]

[math] \vec F=-v^2\frac {2}{5}sen(x-vt)\vec i+\frac {2}{5}(λ+2µ)sen(x-vt)\vec i [/math]

[math]v=\sqrt{λ+2µ}[/math]

Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Modelización_del_comportamiento_de_una_placa_bidimensional&oldid=54281»