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

De MateWiki
Saltar a: navegación, buscar
Trabajo realizado por estudiantes
Título Visualización de campos escalares y vectoriales en fluidos. Grupo 2A
Asignatura Teoría de Campos
Curso 2016-17
Autores Luis de la Fuente Puig, Alfonso Martín de Soto Aláez, Álvaro Pumares García
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura




1 -.Introducción

En este trabajo analizaremos el comportamiento de una placa rectangular plana que ocupa la región [math](x,y) ∈ [-2,2]×[0,4][/math] sometida a dos cantidades físicas: la temperatura [math]T(x,y)[/math] que depende de dos variables escalares [math](x,y)[/math] y los desplazamientos [math]\vec u(x,y)[/math] producidos por la acción de una fuerza determinada. Definimos [math]\vec r_0(x,y)[/math] como el vector de posición de los puntos de la placa antes de la deformación por lo que 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(x,y)=\vec r_0(x,y)+\vec u(x,y)[/math] Suponemos 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(x,y)=-\frac{x^3}{20}\vec i+f(y)\vec j[/math], donde [math]f(y)[/math] es una cierta función.


2 -.Representación del sólido

En primer lugar representaremos los puntos interiores de la placa mediante un mallado. Para ello tomaremos los ejes en el rectángulo [math](x,y) ∈ [-3,3]×[0,5][/math] y como paso de muestreo h=1/10 para las variables x e y. Utilizaremos el siguiente código de matlab en el que el mallado lo representamos mediante el comando mesh y en el que los puntos de la placa están distribuidos de esa forma debido al paso de muestreo que hemos usado.

Placa plana rectagular
%Paso de muestreo
h=0.1;
%Región que ocupa el rectángulo
x=[-2:h:2];
y=[0:h:4];
%Creación del mallado
[X,Y]=meshgrid(x,y);
%Representación del mallado 
mesh(X,Y,X.*0)
axis equal;
%Establecemos el límite de los ejes
axis([-3,3,0,5])


3 -.Temperatura del sólido

En este apartado dibujaremos las curvas de nivel que definen la temperatura cuya función viene dada por: [math]T(x,y)=(x+2)^2+(y+2)^2[/math] Para ello hacemos uso del siguiente código de matlab en el que representaremos la temperatura mediante el comando mesh y las curvas de nivel mediante contour en el que especificamos el número de curvas de nivel. También utilizamos el comando colorbar para establecer la escala de colores que nos definen una temperatura mayor (colores claros) o una temperatura menor (colores oscuros).

Temperatura y curvas de nivel
%Paso de muestro
h=0.1;
%Región que ocupa el rectángulo
x=[-2:h:2];
y=[0:h:4];
%Creación del mallado
[X,Y]= meshgrid(x,y);
%Función temperatura
T=(X+2).^2+(Y+2).^2;
%Establecemos el límite de los ejes
axis([-3,3,0,5]);
axis equal;
%Representación de la temperatura y las curvas de nivel
subplot(1,2,1);
mesh(X,Y,T);
subplot(1,2,2);
contour(X,Y,T,50);
colorbar
%Temperatura máxima
Tmax=max(max(T));



Como podemos observar en el gráfico la temperatura máxima se alcanza en el punto (2,4).Introducimos en la ventana de comandos de matlab para que nos dé la temperatura máxima ya que en el código hemos incluido el comando max(max()) y nos proporciona el siguiente resultado: [math]Tmax=52[/math]

3.1 -.Gradiente de temperatura

Ahora vamos a proceder a calcular el gradiente de temperatura (∇T) y dibujarlo como un campo vectorial. El gradiente de un campo escalar f es un campo vectorial que se calcula mediante: [math]∇f=\frac{∂f}{∂x}\vec i+\frac{∂f}{∂y}\vec j+\frac{∂f}{∂z}\vec k[/math]. Por lo que: [math]∇T=2·(x+2)+2·(y+2)[/math]. En el gráfico se puede apreciar que el gradiente de temperatura es ortogonal a las curvas de nivel de la temperatura ya que una de las propiedades del gradiente es:

  • ∇f es ortogonal a las superfices (o curvas si estamos en el plano) de nivel.

En este apartado nos ayudaremos del siguiente código de matlab en el que dibujaremos las curvas de nivel de la temperatura como ya hemos hecho anteriormente y sobre ellas aplicaremos el comando quiver para representar el campo vectorial correspondiente. También hemos incluido el comando gradient como comentario con el que Matlab calcularía el gradiente de un campo escalar sin necesidad de realizarlo a mano por el alumno.

Gradiente de temperatura
%Paso de muestro
h=0.1;
%Región que ocupa el rectángulo
x=[-2:h:2];
y=[0:h:4];
%Creación del mallado
[X,Y]=meshgrid(x,y);
%Gradiente de T
T=(X+2).^2+(Y+2).^2;
contour(x,y,T)
%[X,Y]=gradient(T,1,1);
dx=2.*(X+2);             %Componente X del gradiente
dy=2.*(Y+2);             %Componente Y del gradiente
hold on
quiver(x,y,dx,dy)
axis equal
hold off




4 -.Campo de desplazamientos

Los desplazamientos de los puntos de la placa vienen dados por el vector: [math]\vec u(x,y)=-\frac{x^3}{20}\vec i+f(y)\vec j[/math] Calculamos la función [math]f(y)[/math] a través de las siguientes condiciones que nos proporcionan en el enunciado:

  • Los puntos situados en el eje y = 0 no sufren desplazamiento en la dirección [math]\vec j[/math].
  • [math]∇·\vec u=−\frac{3x^2}{20}+\frac{1}{10}[/math].

Sabiendo que la divergencia de un campo vectorial es [math]∇·\vec u=\frac{∂u_1}{∂x}+\frac{∂u_2}{∂y}+\frac{∂u_3}{∂z}[/math], igualamos la divergencia del vector [math]\vec u[/math] que nos proporcionan en el enunciado con la divergencia dada en este apartado:: [math]-\frac{3x^2}{20}+f'(y)=-\frac{3x^2}{20}+\frac{1}{10}[/math] Por lo que [math]f'(y)=\frac{1}{10}[/math], integramos y nos queda [math]f(y)=\frac{y}{10}+C[/math]. De la primera condición deducimos que C=0 ya que cuando y=0 la componente [math]\vec j[/math] del desplazamiento se anula. Finalmente el vector [math]\vec u[/math] nos queda de la siguiente manera: [math]\vec u=-\frac{x^3}{20}\vec i+\frac{y}{10}\vec j[/math]


4.1 -.Representación del campo de vectores

A continuación dibujamos el campo de vectores en los puntos del mallado del sólido:

Campo de vectores en los puntos del sólido
%Paso de muestro
h=0.1;
%Región que ocupa el rectángulo
x=[-2:h:2];
y=[0:h:4];
%Creación del mallado
[X,Y]=meshgrid(x,y);
%Componente en X de la divergencia
Ux=-X.^3/20;
%Componente en Y de la divergencia
Uy=Y./10;
%Representación del campo de vectores
hold on
mesh(X,Y,0*X);
quiver(X,Y,Ux,Uy);
hold off
%Establecemos el límite de los ejes
axis([-3 3 0 5]);
axis equal;


4.2 -.Posición del sólido antes y después del desplazamiento

En este apartado procederemos a dibujar el sólido antes del desplazamiento cuya posición viene definida por [math]\vec r_0[/math] y se corresponde con la región que ocupa la placa rectangular plana. A continuación representamos la posición de la placa tras haber aplicado la fuerza sumando a la posición inicial el desplazamiento dado por el campo de vectores [math]\vec u[/math] tal como se había explicado en el enunciado. El vector de posición para cada punto inicial de la placa queda definido por la matriz de coordenadas. Los desplazamientos se producen en la dirección del vector unitario [math]\vec i[/math] y [math]\vec j[/math], esto es en la dirección del eje x e y.

Apartado62Ab.jpg 
%Paso de muestro
h=0.1;
%Región que ocupa el rectángulo
x=[-2:h:2];
y=[0:h:4];
%Creación del mallado
[X,Y]=meshgrid(x,y);
%Matriz de los desplazamientos para el eje "x"
Ux=-(X.^3)./20; 
%Matriz de los desplazamientos para el eje "y"
Uy=Y./10;
%Matriz con las coordenadas en x, después del desplazamiento
Rx=Ux+X;
%Matriz con las coordenadas en y, después del desplazamiento
Ry=Uy+Y;
%Representación del sólido antes del desplazamiento 
subplot(1,2,1);
mesh(X,Y,0*X);
%Establecemos el límite de los ejes
axis([-3 3 0 5]);
title('Placa sin deformación');
xlabel('Eje x');
ylabel('Eje y');
%Representación del sólido tras del desplazamiento
subplot(1,2,2);
mesh(Rx,Ry,0*Rx);
%Establecemos el límite de los ejes
axis([-3 3 0 5]);
title('Placa con deformación');
xlabel('Eje x');
ylabel('Eje y');


5 -.Representación de la divergencia de [math]\vec u[/math]

La divergencia de [math]\vec u[/math] dada en el enunciado del problema para poder calcular [math]f(y)[/math] es: [math]∇·\vec u= −\frac{3x^2}{20}+\frac{1}{10}[/math] Dicha función es una parábola simétrica limitada por la placa. Los puntos de corte con el eje [math]x(∇·\vec u=0)[/math] son x=0,816 y x=-0,816 que coincide con los mínimos ya que a partir de estos puntos la divergencia empieza a ser negativa. El máximo se calcula dando un valor de x=0 a la divergencia y nos sale en x,z(0,1/10) maxdivU=0,1. Podemos observar que el valor de la divergencia en la dirección [math]\vec j[/math] no varía, es decir, para los puntos de la placa que tienen la misma coordenada en x la divergencia es la misma.

En este caso para el código del programa utilizamos el comando surf para representar la divergencia en los puntos de la placa y en el eje z representamos los valores de la divergencia:

Divergencia de [math]\vec u[/math] en los puntos de la placa
%Paso de muestro
h=0.1;
%Región que ocupa el rectángulo
x=[-2:h:2];
y=[0:h:4];
%Creación del mallado
[X,Y]=meshgrid(x,y);
divU=-3/20.*X.^2+1/10;
%Representación de la divergencia en los puntos de la placa
surf(X,Y,divU);
%Establecemos el límite de los ejes
axis([-3 3 0 5]);
axis equal;
hold on 
xlabel('Eje x');
ylabel('Eje y');
zlabel('Eje z');


A continuación dibujamos la gráfica de la divergencia de [math]\vec u[/math] proyectada sobre el plana XZ en la que se puede apreciar sus máximos, mínimos, tramos crecientes y decrecientes. En el programa de matlab en este caso bastaría con el comando plot para dibujarla.

Divergencia de [math]\vec u[/math] proyectada en XZ
%Paso de muestro
h=0.1;
%Región que ocupa el rectángulo
x=[-2:h:2];
y=[0:h:4];
divU=-3/20.*x.^2+1/10;
%Gráfica de la divergencia de u
axis equal;
plot(x,divU)
xlabel('Eje x')
ylabel('Eje z')



6 -.Rotacional de [math]\vec u[/math]

En el cálculo vectorial, el rotacional es un operador vectorial sobre campos vectoriales definidos en un abierto de [math]R^3[/math] que muestra la tendencia de un campo vectorial a inducir rotación alrededor de un punto. Definimos el campo vectorial [math]\vec u[/math] en coordenadas cartesianas como: [math]\vec u=u_1\vec i+u_2\vec j+u_3\vec k[/math] El rotacional de este campo vectorial se calcula mediante: [math] \nabla×\vec u=\left| \begin{matrix} \vec i & \vec j & \vec k \\ & & \\ \cfrac{\partial}{\partial x} & \cfrac{\partial}{\partial y} & \cfrac{\partial}{\partial z} \\ & & \\ u_1 & u_2 & u_3 \end{matrix}\right| [/math]


[math] \nabla×\vec u=\left( \frac{\partial u_3}{\partial y} - \frac{\partial u_2}{\partial z}\right)\vec i + \left(\frac{\partial u_1}{\partial z} - \frac{\partial u_3}{\partial x}\right)\vec j + \left(\frac{\partial u_2}{\partial x} - \frac{\partial u_1}{\partial y}\right)\vec k [/math]

Como [math]u_1=-\frac{x^3}{20}[/math] (solo depende de la variable [math]x[/math]) y [math]u_2=\frac{y}{10}[/math] (solo depende de la variable [math]y[/math]) entonces: [math]\nabla×\vec u=\vec 0[/math] Por lo tanto el campo vectorial [math]\vec u[/math] es irrotacional, en un dominio simplemente conexo de donde se deduce que es conservativo


7 -.Tensiones

Definimos [math]\epsilon(\vec u)=(∇\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 u1+2µ\epsilon[/math] donde [math]λ[/math] y [math]µ[/math] 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. Tomamos [math]λ=µ=1[/math].

Primero obtenemos la [math]∇·\vec u·1[/math]: [math][−\frac{3x^2}{20}+\frac{1}{10}]\begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{bmatrix}=\begin{bmatrix} −\frac{3x^2}{20}+\frac{1}{10} & 0 & 0 \\ 0 & −\frac{3x^2}{20}+\frac{1}{10} & 0 \\ 0 & 0 & −\frac{3x^2}{20}+\frac{1}{10} \end{bmatrix}[/math] En segundo lugar calculamos el [math]∇\vec u[/math]: [math]∇\vec u=\begin{bmatrix} -\frac{3x^2}{20} & 0 & 0 \\ 0 & \frac{1}{10} & 0 \\ 0 & 0 & 0 \end{bmatrix}[/math] A continuación hallamos [math](∇\vec u+∇\vec u^t)[/math]: [math](∇\vec u+∇\vec u^t)=\begin{bmatrix} -\frac{3x^2}{10} & 0 & 0 \\ 0 & \frac{1}{5} & 0 \\ 0 & 0 & 0 \end{bmatrix}[/math] Por último operamos el tensor de tensiones ([math]\sigma[/math]): [math]\sigma=\begin{bmatrix} -\frac{9x^2}{20}+\frac{1}{10} & 0 & 0 \\ 0 & -\frac{3x^2}{20}+\frac{3}{10} & 0 \\ 0 & 0 & -\frac{3x^2}{20}+\frac{1}{10} \end{bmatrix}[/math] Dibujamos las tensiones normales en la dirección que marca el eje [math]\vec i[/math], es decir [math]\vec i·σ·\vec i[/math] que resulta ser [math]\sigma_{11}[/math], las tensiones normales en la dirección que marca el eje [math]\vec j[/math], es decir [math]\vec j·σ·\vec j[/math] que coincide con [math]\sigma_{22}[/math] y las correspondientes al eje [math]\vec k[/math], es decir [math]\vec k·σ·\vec k[/math] cuyo resultado es [math]\sigma_{33}[/math].

Para ello empleamos el siguiente código de matlab: 

%Paso de muestreo
h=0.1;
%Región que ocupa el rectángulo
x=[-2:h:2];
y=[0:h:4];
%Creación del mallado
[X,Y]=meshgrid(x,y);
%Tension perpendicular a i
subplot(1,3,1)
ti=1/20.*(-9.*X.^2)+1/10;
surf(X,Y,ti)
axis([-3 3 0 5])
xlabel('Eje x')
ylabel('Eje y')
zlabel('Eje z')
%Tension perpendicular a j
subplot(1,3,2)
tj=1/20.*(-3.*X.^2)+3/10;
surf(X,Y,tj)
view(3)
axis([-3 3 0 5])
xlabel('Eje x')
ylabel('Eje y')
zlabel('Eje z')
%Tension perpendicular a k
subplot(1,3,3)
tk=1/20.*(-3.*X.^2)+1/10;
surf(X,Y,tk)
view(3)
axis([-3 3 0 5])
xlabel('Eje x')
ylabel('Eje y')
zlabel('Eje z')
Tensiones perpendiculares a x, y, z


7.1 -.Tensiones tangenciales

Una vez calculado el tensor de tensiones, vamos a proceder a hallar las tensiones tangenciales dadas por la siguiente fórmula: [math]|\sigma·\vec i-(\vec i·\sigma·\vec i)\vec i)|[/math] Operando con la matriz [math]\sigma[/math] obtenida en el apartado anterior y tomando el vector [math]\vec i[/math] como la dirección del eje x: [math]|\begin{bmatrix} -\frac{9x^2}{20}+\frac{1}{10} & 0 & 0 \\ 0 & -\frac{3x^2}{20}+\frac{3}{10} & 0 \\ 0 & 0 & -\frac{3x^2}{20}+\frac{1}{10} \end{bmatrix}\begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}-[-\frac{9x^2}{20}+\frac{1}{10}]\begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}|=0[/math] Se deduce analíticamente que la tensión tangencial al plano ortogonal a [math]\vec i[/math] es nula.


8 -.Tensión de Von Mises

La tensión de Von Mises se define por la fórmula: [math]\sigma_{VM} = \sqrt{\frac{(\sigma_1-\sigma_2)^2+ (\sigma_2-\sigma_3)^2+ (\sigma_3-\sigma_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 indica cuando un material pasa de comportamiento elástico a plástico.

Utilizamos el siguiente código de matlab en el que hemos hecho uso del comando eig que nos saca los autovalores de la matriz tensor de tensiones: 

%Paso de muestreo
h=0.1;
%Región que ocupa el rectángulo
x=[-2:h:2];
y=[0:h:4];
%Creación del mallado
[X,Y]=meshgrid(x,y);
%Matriz vacía de inicialización
valVM=[]; 
%Bucle para sacar autovalores de la matriz tensor de tensiones
for x=[-2:h:2]
    a=1/20.*(-9.*x.^2)+1/10;
    b=1/20.*(-3.*x.^2)+3/10;
    c=1/20.*(-3.*x.^2)+1/10;
    mat=[a 0 0;0 b 0;0 0 c];                                   
    VM=eig(mat); 
    tVM=sqrt(((VM(1)-VM(2))^2+(VM(2)-VM(3))^2+(VM(3)-VM(1))^2)/2);
    valVM=[valVM,tVM];
end
%Gráfica de las tensiones    
plot(valVM)
ylabel('Valores de Tension de Von Mises')
%Tensión máxima de Von Mises
max(max(valVM))
center

El punto en el que se alcanza mayor valor de la tensión de Von Mises es en x=-2 y x=2 y la tensión máxima es 1.3115.


9 -.Campo de fuerzas sobre la placa

Para obtener el campo de fuerzas [math]\vec F[/math] usamos la ecuación de la elasticidad lineal: [math] \vec F =-\nabla\cdot \sigma=-\frac{\partial \sigma_{ji}}{\partial x_{j}}\vec e_{i}[/math]. Siendo [math]\sigma[/math] la matriz tensor de tensiones podemos ver que la ecuación de la elasticidad lineal solo tendrá valor distinto de 0 cuando i,j sean 1.: [math]F_x=-\frac{\partial \sigma_{11}}{\partial x_{1}}\vec e_{1}=-\frac{\partial (-\frac{9x^2}{20}+\frac{1}{10})}{\partial x}\vec i=\frac {9x}{10}\vec i[/math] Dibujamos el campo de fuerzas [math]\vec F[/math] mediante el siguiente código de matlab: 

%Paso de muestreo
h=0.1;
%Región que ocupa el rectángulo
x=[-2:h:2];
y=[0:h:4];
%Componentes de la fuerza
Fx=9/10.*X;
Fy=0.*Y;
%Fz=0
%Dibujo del campo de fuerzas
quiver(X,Y,Fx,Fy)
%Establecemos el límite de los ejes
axis([-3 3 0 5])
Campo vectorial [math]\vec F[/math]

Si comparamos la placa con deformación del apartado 4.2 podemos observar que se ha contraído y en la gráfica obtenida vemos que la fuerza lleva sentido opuesto a la deformación. Esto se debe al signo negativo de la ecuación de la elasticidad lineal.


10 -.Masa total de la placa

Suponemos que la densidad de la placa es [math]d(x,y)=1+e^{-|x|/(y+1)^2}[/math]. Vamos a proceder a calcular la masa de la placa a través del método del trapecio: 

%Paso de muestreo
h=0.1;
%Región que ocupa el rectángulo
x=-2:h:2;
y=0:h:4;
%Número de puntos
N1=(2-(-2))/h; N2=(4-0)/h;  
%Extremos del intervalo
a=-2; b=2; c=0; d=4;
%Creación del mallado
[X,Y]=meshgrid(x,y);  
%Función densidad
d=1+exp(-abs(X)./((Y+1).*2));
%Masa en los puntos de la malla
w1=ones(N1+1,1);                 
w1(1)=1/2; w1(N1+1)=1/2;
w2=ones(N2+1,1);                 
w2(1)=1/2; w2(N2+1)=1/2;
masa=h*h*w2'*d*w1


La masa total de la placa es 29.2439.

Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Visualización_de_campos_escalares_y_vectoriales_en_elasticidad_(Grupo_2A)&oldid=38402»