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

Trabajo realizado por estudiantes
Título	Visualización de campos escalares y vectoriales en elasticidad. Grupo 2-C
Asignatura	Teoría de Campos
Curso	2014-15
Autores	Victoria Manget Sánchez Sacristán, Alejandro Perales Juidías, Manuel Hernández Muñoz, Cristina Martínez Navarro y Alicia Vilalta Duce
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

Trabajo realizado por:

• Victoria Manget Sánchez Sacristán
• Alejandro Perales Juidías
• Manuel Hernández Muñoz
• Cristina Martínez Navarro
• Alicia Vilalta Duce

1 Enunciado

Mallado de los trabajos 1 y 2

Consideramos una placa rectangular plana (en dimensión 2) que ocupa la región [math] [-1/2,1/2] \times [0,4][/math]. En ella vamos a suponer que tenemos definidas dos cantidades físicas: la temperatura [math]T(x,y)[/math], que depende de las dos variables espaciales [math](x,y)[/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]r_0(x,y)[/math] el vector de posición de los puntos de la placa antes de la deformación, 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] 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(x,y) = \vec a \dot(\vec b \cdot \vec r_0)^2, [/math] donde [math]\vec a[/math] y [math]\vec b[/math] son vectores dados. En este trabajo supondremos lo siguiente: [math] \vec a=\frac{\vec j}{10}, \vec b=\vec j. [/math]

Se pide:

1. Dibujar un mallado que represente los puntos interiores del sólido (ver figura 1). Tomar como paso de muestreo h= 1/10 para las variables x e y.
2. La temperatura del sólido viene dada por la función \( T(x,y)=(8-y^2+2y)e^{-x^2}\). Dibujar las curvas de nivel y decidir en qué punto la temperatura es máxima.
3. Calcular \(\nabla T\) y dibujarlo como campo vectorial. Observar gráficamente que \(\nabla T\) es ortogonal a dichas curvas.
4. Consideramos ahora el campo de desplazamientos [math]\vec u(x,y)[/math]. Dibujar el campo de vectores en los puntos del mallado sólido.
5. Dibujar el sólido antes y después del desplazamiento dado por el campo de vectores [math]\vec u(x,y)[/math]. Dibujar ambos en la misma figura usando el comando subplot.
6. Calcular la divergencia de [math]\vec u(x,y)[/math] en todos los puntos del sólido y dibujarla. ¿Qué puntos tienen mayor divergencia? La divergencia es una medida del cambio de volumen local debido al desplazamiento. ¿Se puede apreciar esto en la gráfica?
7. Calcular |\(\nabla T\)x[math]\vec u(x,y)[/math]| en todos los puntos del sólido y dibujarlo. ¿Qué puntos sufren un mayor rotacional?
8. Definamos [math]\epsilon(\vec u)=(\nabla \vec u + \nabla \vec u^t)/2[/math], la parte simétrica del tensor gradiente de [math]\vec u(x,y)[/math]. En un medio elástico lineal, isótropo y homogéneo los desplazamientos permiten escribir el tensor de tensiones [math]\sigma_{ij}[/math] a través de la fórmula

                 [math]\sigma_{i,j}=λ \cdot \nabla \vec u 1 + 2\cdotμ\cdot\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. Tomando [math]λ=μ=1[/math], dibujar las tensiones normales en la dirección que marca el eje [math]\vec i[/math], es decir [math](\vec i \cdot \sigma \cdot \vec i)[/math] y las tensiones normales en la dirección que marca el eje [math]\vec j[/math], es decir [math](\vec j \cdot \sigma \cdot \vec j)[/math]. Comparar las gráficas de las tensiones normales con las del módulo de la divergencia y las del módulo del rotacional.

9. Dibujar las tensiones tangenciales respecto al plano ortogonal a [math]\vec i[/math], es decir [math]|\sigma \cdot \vec i-(\vec i \cdot \sigma \cdot \vec i) \vec i|= 0 [/math]. ¿Dónde son mayores? Comparar las gráficas de las tensiones tangenciales con las del módulo de la divergencia y las del módulo del rotacional.

10. Dibujar las tensiones tangenciales respecto al plano ortogonal a [math]\vec j[/math], es decir [math]|\sigma \cdot \vec j-(\vec j \cdot \sigma \cdot \vec j) \vec j|= 0 [/math]. ¿Dónde son mayores? Comparar las gráficas de las tensiones tangenciales con las del módulo de la divergencia y las del módulo del rotacional.

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

donde [math]\sigma_{1}[/math], [math]\sigma_{2}[/math], [math]\sigma_{3}[/math] son los autovalores de [math]\sigma_{}[/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 (y no elástico puro). Pintar la tensión de Von Mises y señalar en qué punto se alcanza el mayor valor. (Para calcular autovalores con OCTAVE/Matlab usar el comando eig.m).

12. Supongamos que la densidad de la placa es \( d(x,y,z)=xy log{(x+2)}\).

2 Influencia de campos físicos sobre el sólido elástico

2.1 Introducción al método

Un campo físico representa la distribución de una magnitud física que varía en el espacio, en otras palabras, es una aplicación que asocia a cada punto del espacio un tensor (hablamos de campos escalares y vectoriales). Normalmente relacionamos estos tensores con variables físicas (temperatura, velocidad, aceleración, etc.). A continuación se estudiará la influencia que ejercen los campos sobre el sólido elástico, capaz de recuperar su deformación (y volver al estado original) tras ser sometido a tensión.

Partimos de una placa rectangular plana cuyos puntos vienen dados por el vector de posición [math] \vec r (x,y)= \vec r_{0}(x,y)+\vec u(x,y) [/math], donde [math] \vec r_{0}(x,y)[/math] son las coordenadas del punto definidas en [math] [-1/2,1/2] \times [0,4][/math] y [math] \vec u(x,y) [/math] refleja el "cambio de posición" (desplazamiento). Consideramos ahora dos campos físicos aplicados a ésta:

Campo escalar temperatura: [math]T(x,y)[/math]

Campo vectorial desplazamiento: [math]\vec u(x,y)[/math]

que dependen de las variables espaciales (x,y).

x= linspace(-1/2,1/2,10)    %Región de rango x entre los valores dados
y= linspace(0,4,40)         %Región de rango y
[X,Y]= meshgrid(x,y)        %Mallado de la placa 
Z= zeros(40,10)             %Altura de los puntos del intervalo
surf(X,Y,Z)                 %Representación de la superficie
axis equal                  %Caracterización de los ejes   
axis([-1,1,-1,5])           %Limitamos una región en los ejes

2.2 Influencia de la temperatura

Imaginemos que sobre el solido actúa un primer campo que representa la temperatura puntual y no depende del tiempo. Definimos su foco de calor con centro en el origen (0,0) mediante la expresión: \( T(x,y)=(8-y^2+2y)e^{-x^2}\)

Su representación gráfica sobre la placa viene dada por una superficie en tres dimensiones donde los tonos rojizos indican temperaturas más altas y los tonos azules temperaturas más bajas (ver figura 2).

La visualización en dos dimensiones nos permite ver de forma más clara la relación entre temperatura-distancia al foco. Es fácil comprobar que a medida que nos acercamos al foco la temperatura asciende (ver figura 1).

x= linspace(-1/2,1/2,10)                            %Región de rango x entre los valores dados
y= linspace(0,4,40)                                 %Región de rango y
[X,Y]= meshgrid(x,y)                                %Mallado de la placa
F=inline('(8-y.^2+2.*y).*exp(-x.^2)','x','y');      %Función de temperatura
Z1=F(X,Y)                                           %Elevación de los puntos de la función de temperatura
axis([-1,1,-1,5])                                   %Caracterización de los ejes
subplot(1,2,1)                                      %Creación de subventana en el espacio de representación
surf(X,Y,Z1)                                        %Representación superficie en 3D
view([0,0,1])                                       %Vista desde arriba (paso a 2D)
subplot(1,2,2)                                      %Cambio a subventana 2
mesh(X,Y,Z1)                                        %Representación en mallado de la superficie
max(max(F(X,Y)))                                    %Calculo del máximo de la matriz

El valor máximo de la temperatura es 8.9716.

Figura 1 y figura 2

Para obtener la dirección de máximo cambio de temperatura sobre la superficie (hacia dónde tengo que moverme para calentarme lo más rápido posible) nos valemos del gradiente, representado por un campo vectorial superpuesto a las curvas de nivel de la superficie. Se observa ortogonalidad entre vectores y curvas de nivel. Esto se debe a que la manera de recorrer un tramo sobre la superficie se realiza de la forma más corta en esa dirección, la del gradiente, que representa la mínima distancia entre dos puntos de diferente temperatura.

x= linspace(-1/2,1/2,10)                            %Región de rango x entre los valores dados
y= linspace(0,4,40)                                 %Región de rango y
[X,Y]= meshgrid(x,y)                                %Mallado de la placa
T=(8-Y.^2+2.*Y).*exp(-(X.^2));                      %Función de temperatura
[X,Y]=gradient(T,0.1,0.1);                          %Cálculo del gradiente
contour(x,y,T)                                      %Representación de las curvas de nivel
hold on
quiver(x,y,X,Y)                                     %Representación de los vectores del campo gradiente
hold off

Figura 3

2.3 Influencia del desplazamiento

Consideremos ahora [math]\vec u(x,y) = \vec a \dot(\vec b \cdot \vec r_0)^2[/math] un campo vectorial que determine el desplazamiento de los puntos del mallado. Al aplicarlo sobre el conjunto, representa el cambio de posición de cada uno de los puntos del sólido.

x= linspace(-1/2,1/2,10)                            %Región de rango x entre los valores dados
y= linspace(0,4,40)                                 %Región de rango y
[X,Y]= meshgrid(x,y)                                %Mallado de la placa
i=1; j=1;                                           %Creación del comienzo del bucle
Ux= [];                                             %Creación de la matriz Ux
Uy=[];                                              %Creación matriz Uy
for i=1:40                                          %Bucle para obtener matrices Ux, Uy con respecto a los valores i,j
  for j=1:10
    Ux(i,j)=0;
    Uy(i,j)=(Y(i,j).^2)/10;
  end
end
quiver(X,Y,Ux,Uy)                                   %Representación de los vectores de desplazamiento
axis equal

Figura 4

Al aplicar el campo de desplazamiento a la placa rectangular observamos su posición final en comparación a su posición inicial.

Figura 5

2.4 Influencia de la divergencia

La divergencia refleja el incremento local de volumen del sólido sometido al campo vectorial [math]\vec u[/math] , podemos entenderlo como el flujo por unidad de volumen. Gráficamente imaginamos un prisma de base cuadrada y longitud de lado tendiendo a 0. La variación de volumen sufrida en ese punto (a consecuencia del campo) será la divergencia en dicho punto. Resultará mayor cuando el campo aplicado produzca un mayor cambio en las dimensiones de la superficie y máxima en aquellos puntos en los que ejerza mayor influencia (aumente más el volumen). De este modo será mayor la divergencia de un campo que represente un sistema de fuerzas que un campo que simplemente desplace la superficie.

x= linspace(-1/2,1/2,10)          %Región de rango x entre los valores dados
y= linspace(0,4,40)               %Región de rango y
[X,Y]= meshgrid(x,y)              %Mallado de la placa
i=1; j=1;
Ux= [];
Uy=[];
for i=1:40
  for j=1:10
    Ux(i,j)=0;                    %Bucle para obtener matrices Ux, Uy con respecto a los valores i,j
    Uy(i,j)=(Y(i,j).^2)/10; 
  end
end
div= divergence(X,Y,Ux,Uy)        %Divergencia del campo vectorial u de desplazamiento
mesh(div)                         %Representación de la divergencia en mallado
max(max(div))                     %Cálculo del máximo valor de la divergencia

El valor máximo obtenido de la divergencia es 0.7897.

Figura 6

2.5 Influencia del rotacional

El rotacional de un campo vectorial representa la capacidad de girar en torno a un eje. Los puntos que sufren un mayor rotacional serán aquellos en los cuales las derivadas parciales sean máximas ya que el rotacional se halla a partir del producto vectorial de las componentes del campo y los vectores de la base, siendo estas dos fijas.

Procedemos a calcular el rotacional de [math]\vec u(x,y)[/math] y observamos que el resultado es 0.

[math]|\nabla\times\vec u|=\frac{1}{\sqrt{g}} \begin{vmatrix} \vec{g_i} & \vec{g_j} & \vec{g_k} \\ \frac{ \partial}{\partial x} & \frac{\partial }{\partial y} & \frac{\partial }{\partial z}\\ u_1 & u_2 & u_3 \end{vmatrix}= \begin{vmatrix} \vec{g_i} & \vec{g_j} & \vec{g_k} \\ \frac{ \partial}{\partial x} & \frac{\partial }{\partial y} & \frac{\partial }{\partial z}\\ \ 0 & \frac{y^2}{10} & 0 \end{vmatrix}=0[/math]

Por tanto el campo es irrotacional, en un dominio simplemente conexo, es decir, se puede demostrar que el campo vectorial u⃗ es conservativo.

2.6 Tensiones normales de un sólido elástico

Las tensiones sobre un sólido vienen definidas por un tensor que depende de los coeficientes de Lamé (que son dos constantes elásticas que caracterizan por completo el comportamiento elástico de un sólido isótropo para pequeñas deformaciones).

A continuación expondremos la gráfica en la que comparamos las tensiones normales con las del módulo de la divergencia (Figura 7).

Para hallar las tensiones normales en la dirección del eje [math]\vec i[/math] aplicamos [math](\vec i \cdot \sigma \cdot \vec i)[/math] resultando y/5.

Para hallar las tensiones normales en la dirección del eje [math]\vec j[/math] aplicamos [math](\vec j \cdot \sigma \cdot \vec j)[/math] también es 3y/5.

x= linspace(-10/2,1/2,10);          %Región de rango x entre los valores dados
y= linspace(-10,4,40);              %Región de rango y
[X,Y]= meshgrid(x,y);               %Mallado de la placa
sigmax= Y./5                        %Tensión eje x
sigmay=3.*Y./5                      %Tensión eje y
subplot(1,2,1)                      %creación de dos subventanas en el espacio de representación
mesh(sigmax)                        %Representación matriz tensiones en x
hold on             
mesh(sigmay)                        %Representación matriz tensiones en y
hold off
subplot(1,2,2)                      %Cambio a la subventana 2
i=1; j=1;
Ux= [];
Uy=[];
for i=1:40
  for j=1:10
    Ux(i,j)=0;                      %Bucle para obtener matrices Ux, Uy con respecto a los valores i,j
    Uy(i,j)=(Y(i,j).^2)/10; 
  end
end
div= divergence(X,Y,Ux,Uy)        %divergencia del campo vectorial u de desplazamiento
mesh(div)                         %Representación de la divergencia
axis([0,10,0,40,-6,4])            %Caracterización de los ejes

Figura 7

Aplicado a la placa rectangular puesto que tiene un rotacional nulo, por lo explicado anteriormente, dependerá solo de [math]\sigma_{ij}[/math] siendo ésta la parte antisimétrica. Se puede observar analíticamente que las componentes de sigma son las mismas que las de [math]\sigma_{}[/math].

Para hallar las tensiones normales en la dirección del eje [math]\vec i[/math] ,aplicamos [math](\vec i \cdot \sigma \cdot \vec i)[/math] resultando cero pues es la proyección del vector [math](\sigma \cdot \vec i)[/math] al solo depender de la coordenada [math]\vec j[/math] sobre la dirección [math]\vec i[/math] es 0.

Para hallar las tensiones normales en la dirección del eje [math]\vec j[/math] aplicamos [math](\vec j \cdot \sigma \cdot \vec j)[/math] también es cero, puesto que el vector [math](\sigma \cdot \vec j)[/math] solo depende de la componente [math]\vec i[/math] por lo que la proyección sobre el eje [math]\vec j[/math] es nula.

Por ser nulo no obtendremos ninguna comparación respecto al rotacional.

2.7 Tensiones tangenciales de un sólido elástico

Las tensiones tangenciales (componentes del vector tensión rasantes al plano de corte) vienen definidas por la ecuación [math]|\sigma \cdot \vec i-(\vec i \cdot \sigma \cdot \vec i) \vec i|[/math] si son respecto al plano ortogonal a [math]\vec i[/math] y por la ecuación [math]|\sigma \cdot \vec j-(\vec j \cdot \sigma \cdot \vec j) \vec j|[/math] respecto al plano ortogonal a [math]\vec j[/math]. Resolviendo la fórmula descrita a continuación para el cálculo de tensiones tangenciales podemos observar que éstas son nulas.

Esto se debe a que estamos proyectando sobre la dirección [math]\vec i[/math] un vector con componente únicamente en la dirección [math]\vec i[/math], por tanto coincide con la tensión en la misma dirección y al restarlos, ambos se anulan. Análogamente para la proyección sobre la dirección [math]\vec j[/math]. [math] |\sigma \cdot \vec i-(\vec i\cdot\sigma\cdot\vec i)\cdot\vec i| [/math]

[math]|\left( \begin{array}{ll} \frac{y}{5}\ & 0 & 0 \\ 0 & \frac{3y}{5} & 0 \\ 0 & 0 & \frac{y}{5} \end{array} \right)\cdot\left( \begin{array}{l} 1 & 0 & 0 \end{array} \right)-( \left( \begin{array}{l} 1 \\ 0 \\ 0 \end{array} \right)\cdot\left( \begin{array}{ll} \frac{y}{5}\ & 0 & 0 \\ 0 & \frac{3y}{5} & 0 \\ 0 & 0 & \frac{y}{5} \end{array} \right)\cdot\left( \begin{array}{l} 1 & 0 & 0 \end{array} \right))\cdot\left( \begin{array}{l} 1 & 0 & 0 \end{array} \right)|=0 [/math]

[math]|\left( \begin{array}{ll} \frac{y}{5}\ & 0 & 0 \\ 0 & \frac{3y}{5} & 0 \\ 0 & 0 & \frac{y}{5} \end{array} \right)\cdot\left( \begin{array}{l} 0 & 1 & 0 \end{array} \right)-( \left( \begin{array}{l} 0 \\ 1 \\ 0 \end{array} \right)\cdot\left( \begin{array}{ll} \frac{y}{5}\ & 0 & 0 \\ 0 & \frac{3y}{5} & 0 \\ 0 & 0 & \frac{y}{5} \end{array} \right)\cdot\left( \begin{array}{l} 0 & 1 & 0 \end{array} \right))\cdot\left( \begin{array}{l} 0 & 1 & 0 \end{array} \right)|=0 [/math]

2.8 La tensión de Von Mises

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

donde σ1, σ2 y σ3 son los autovalores de σ (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 (y no elástico puro).

sigmatotal=[];                                                          %Vector vacío para la primera interacción.
for y=linspace(0,4,40);                                                 %Bucle para crear el vector sigmatotal
    sigma=[y/5 0 0;0 3*y/5 0;0 0 y/5];                                  %Matriz sigma
    E=eig(sigma);                                                       %Función para sacar los autovalores de la matriz sigma
    sigmavm= sqrt(((E(1)-E(2))^2+(E(2)-E(3))^2+(E(3)-E(1))^2)/2);       %Ecuación de la tensión de Von Mises
    sigmatotal=[sigmatotal,sigmavm]                                     %Creación del vector sigmatotal a partir de las tensiones de las interacciones anteriores
end
plot(sigmatotal)                                                        %Dibujo de la tensión

Figura8

2.9 Cálculo de la masa total

h= 1/10                     %Paso de muestreo
u= -1/2:h:1/2               %Vector u
v= 0:h:4                    %Vector v
[X,Y]=meshgrid(u,v)         %Mallado a partir de los vectores u y v
d=abs(X.*Y.*log(X+2))       %Función de densidad
a=h.^2.*d                   
masa=sum(sum(a))            %Cálculo de la integral

El valor de la masa = 1.6578.