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

2016-12-05T09:25:25Z

Alvaro Pumares Garcia:

{{ TrabajoED | Visualización de campos escalares y vectoriales en fluidos. Grupo 2A | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC16/17|2016-17]] | Luis de la Fuente Puig, Alfonso Martín de Soto Aláez, Álvaro Pumares García}}
[[Categoría:Teoría de Campos]]
[[Categoría:TC16/17]]

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

==-.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.
[[Archivo:Apartado12A.jpg|300px|thumb|right|Placa plana rectagular]]
{{matlab|codigo=
%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])
}}

==-.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).
[[Archivo:Apartado22Ab.jpg|500px|thumb|right|Gradiente de temperatura]]
{{matlab|codigo=
%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>

===-.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.
[[Archivo:Apartado32A.jpg|500px|thumb|right|Temperatura y curvas de nivel]]
{{matlab|codigo=
%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
}}

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

===-.Representación del campo de vectores===

A continuación dibujamos el campo de vectores en los puntos del mallado del sólido:
[[Archivo:Apartado52Ab.jpg|450px|thumb|right|Campo de vectores en los puntos del sólido]]
{{matlab|codigo=
%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;
}}

===-.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 i y j, esto es en la dirección del eje x e y.
[[Archivo:Apartado62Ab.jpg|600px|thumb|right|]]
{{matlab|codigo=
%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');
}}

==-.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= −3x^2/20 + 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 y 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:
[[Archivo:Apartado72Ab.jpg|400px|thumb|right|Divergencia de ū en los puntos de la placa]]
{{matlab|codigo=
%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, puntos de crecimiento y decrecimiento. En el programa de matlab en este caso bastaría conel comando ''plot'' para dibujarla.
[[Archivo:Apartado72Ac.jpg|300px|thumb|right|Divergencia de ū proyectada en XZ]]
{{matlab|codigo=
%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')
}}

==-.Rotacional de <math>\vec u</math>==

En el cálculo vectorial, el rotacional o rotor 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

==-.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 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 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:
{{matlab|codigo=
%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')
}}
[[Archivo:Apartado92Ab.jpg|800px|thumb|center|Divergencia de ū proyectada en XZ]]

===-.Tensiones tangenciales===
jaskdjsd

MateWiki - Contribuciones del usuario [es]

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