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

Trabajo realizado por estudiantes
Título	Visualización de campos escalares y vectoriales en elasticidad. Grupo C13
Asignatura	Teoría de Campos
Curso	2014-15
Autores	Pablo Morales Santón, Gonzalo Sánchez Juez, José Mª Llorente Martiáñez, Alonso Herranz Hudson, Rodrigo Bellot Rodriguez
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

ENUNCIADO DEL TRABAJO:

Consideramos una placa plana (en dimensión 2) que ocupa la región \( [-1/2,1/2] \times [0,4]\). En ella vamos a suponer que tenemos definidas dos cantidades físicas: la temperatura \(T(x,y)\), que depende de las dos variables espaciales \((x,y)\), y los desplazamientos \(\vec u(x,y)\). De esta forma, si definimos \(r_0(x,y)\) el vector de posición de los puntos de la placa antes de la deformación, la posición de cada punto \((x,y)\) de la placa en un instante de tiempo \(t\) viene dada por\[ \vec r (x,y)= \vec r_{0}(x,y)+\vec u(x,y). \] 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\[ \vec u(x,y) = \vec a (\vec b \cdot \vec r_0)^2, \] donde \(\vec a\) y \(\vec b\) son vectores dados. En este trabajo suponemos lo siguente\[ \vec a = \frac{\vec i}{20}, \ \vec b = \vec j \] Por lo tanto, para el caso particular de este trabajo el vector de desplazamiento \(\vec u(x,y)\) será\[ \vec u(x,y) = \vec a (\vec b \cdot \vec r_0)^2=\vec a (y)^2= \frac{y^2}{20}\vec i \]

1 Mallado

Nuestro mallado representa una placa rectangular que ocupa la región [math] [-1/2,1/2] \times [0,4] [/math] , tomando como paso de muestreo h=[math]1/10[/math] para las variables x e y.

Mallado de la placa

% Generamos los vectores
x=-0.5:0.1:0.5;       
y=0:0.1:4;   
% Generamos la malla
[Mx,My]=meshgrid(x,y); 
figure(1) 
% Pintamos la malla
mesh(Mx,My,0*Mx)  
% Limitamos una región en los ejes
axis([-2,2,-1,5])
 % Vista superior
view(2)

2 Temperaturas

La función que define la temperatura de la placa es: [math] T(x,y,z)=(8-y^2+2y)e^{-(x+1)^2} [/math] En la imagen están representadas las curvas de nivel. Las lineas rojas representan los valores máximos, y las azules los mínimos. Observamos que es máxima en el punto [math][-1/2,1][/math],

Curvas de nivel de la temperatura en la placa

% Generamos los vectores
x=-0.5:0.1:0.5;      
y=0:0.1:4;    
% Hacemos la malla
[Mx,My]=meshgrid(x,y);
% Escribimos la función temperatura
T=inline('(8 - y.^2 + 2*y).*exp(-(x + 1).^2)','x','y');
% Aplicamos la función a las matrices de la malla
TT=T(Mx,My);
% Representamos la función
axis equal
p=contour(Mx,My,TT,60);
axis([-0.5,0.5,0,4])

2.1 Gradiente de la temperatura

El gradiente de la temperatura representa la dirección en el que el crecimiento de esta es máximo. La función temperatura es: [math] T(x,y,z)=(8-y^2+2y)e^{-(x+1)^2} [/math] Calculamos su gradiente: [math] \nabla T=\frac{\partial T}{\partial x}\vec i+\frac{\partial T}{\partial y}\vec j+\frac{\partial T}{\partial z}\vec k=(8-y^2+2y)e^{-(x+1)^2}[-2(x+1)]\vec i+e^{-(x+1)^2}(-2y+2)\vec j \\ =(8-y^2+2y)(-2x-2)e^{-(x+1)^2}\vec i + (-2y+2)e^{-(x+1)^2}\vec j [/math] Se observa gráficamente, mediante su representación en MatLab u OCTAVE UPM:

Representación del gradiente sobre la placa

% Generamos los vectores
x=-0.5:0.1:0.5;       
y=0:0.1:4;       
% Hacemos la malla
[Mx,My]=meshgrid(x,y);
% Definimos la función temperatura
T=inline('(8 - y.^2 + 2*y).*exp(-(x + 1).^2)','x','y');
TT=T(Mx,My);
% Definimos las componentes del gradiente de la temperatura
Tx=inline('(8 - y.^2 + 2*y).*exp(-(x + 1).^2).*(-2*(x + 1))','x','y');
Ty=inline('(-2*y + 2).*exp(-(x + 1).^2)','x','y');
TTx=Tx(Mx,My);
TTy=Ty(Mx,My);
% Representamos el campo vectorial generado por el gradiente
quiver(Mx,My,TTx,TTy)
% Unimos en unos ejes el gradiente con las lineas de nivel
hold on 
% Dibujamos las lineas de nivel
contour(Mx,My,TT,15) 
axis equal
axis([-2,2,-1,5]) 
view(2)

Observamos que los vectores del gradiente son ortogonales a las curvas de nivel.

3 Desplazamientos

Siendo el vector desplazamiento el vector [math]\vec u(x,y)[/math] calculado anteriormente: [math] \vec u(x,y)=\frac{y^2}{20}\vec i [/math]

El código MatLab del campo de vectores y su representación se muestran a continuación:

Campo de vectores [math]\vec u(x,y)[/math]

% Generamos los vectores
x=-0.5:0.1:0.5;       
y=0:0.1:4; 
% Hacemos la malla
[Mx,My]=meshgrid(x,y); 
% Introducimos la función del campo
T=inline('(y.^2)/20','y');
TT=T(My);
% Representamos el campo
quiver(Mx,My,TT,TT*0)
axis([-1,1,-1,5])

3.1 Representación del desplazamiento en el sólido

La primera imagen muestra la placa sin aplicarle el desplazamiento, mientras que en la segunda ya sí que se le aplica. En su representación claramente se aprecia el cambio en la forma de la placa.

Sólido antes y después de la deformación

% Intervalos de trabajo
x=[-1/2:0.1:1/2];
y=[0:0.1:4];
% Mallado y z=0
[Mx,My]=meshgrid(x,y);
[m,n]=size(Mx);
z=zeros(m,n);
% Representación del sólido sin deformar
subplot(1,2,1)
mesh(Mx,My,z)
axis equal
view(2)
% Campo u
u=(My.^2)/20;
% Sólido después de la deformación
Mx2=Mx+u;
My2=My;
Mz2=0.*Mx2;
% Representación del sólido deformado
subplot(1,2,2)
mesh(Mx2,My2,Mz2)
axis equal
view(2)

3.2 Divergencia

A la hora de calcular la divergencia, estamos calculando el cambio de área local del sólido elástico cuando le aplicamos el campo vectorial [math] \vec u [/math] en ese punto. Si ésta es positiva, tendremos una fuente (aumento de área), y si es negativa será un sumidero (disminución del área).

[math] \nabla * \vec u = \frac{\delta u_1 }{\delta x} + \frac{\delta u_2 }{\delta y} + \frac{\delta u_3 }{\delta z} = 0 [/math]

Como el campo vectorial solo depende de la componente [math] \vec i [/math] y de la coordenada [math] y [/math], al obtener la divergencia tenemos que la derivada parcial es cero, sacando como conclusión que el incremento de área del sólido elástico es nulo, es decir, que solo sufre desplazamientos horizontales en el sentido de [math] \vec i [/math]. Esto no provoca cambio de área, ya que los puntos con la misma coordenada [math] y [/math] se encuentran siempre a la misma distancia unos de otros, por lo que el área de cada sección infinitesimal será [math] dA= xdy [/math], siendo la [math] x [/math] una constante, y cuya integración da el mismo resultado tanto para el sólido sin deformar como para el sólido deformado, al ser los límites de integración los mismos.

3.3 Rotacional

Representación del rotacional

El cálculo del rotacional se realiza analíticamente y en coordenadas cartesianas, tal y como se muestra a continuación, para su posterior representación gráfica con MatLab o OCTAVE UPM: [math] \nabla \times \vec u = \begin{vmatrix} \vec i & \vec j & \vec k \\ \frac{\delta}{\delta x} & \frac{\delta}{\delta y} & \frac{\delta}{\delta z} \\ \frac{y^2}{20} & 0 & 0 \end{vmatrix} = \frac{-y}{10} \vec k [/math] siendo por lo tanto el módulo del rotacional: [math] \left | \nabla \times \vec u \right | = \frac{y}{10} [/math] Una vez calculado el valor del rotacional, se escribe un código MatLab que permita su representación, para posteriormente dar una interpretación. Para su mejor comprensión, se añaden también las imágenes del sólido antes y después de la deformación.

Sólido antes y después de la deformación

% Intervalos de trabajo
x=[-1/2:0.1:1/2];
y=[0:0.1:4];
% Mallado
[Mx,My]=meshgrid(x,y);
% Definición del rotacional
rot=My/10;
% Dibujo de los gráficos
subplot(2,1,1)
mesh(Mx,My,rot)
subplot(2,1,2)
surf(Mx,My,rot)

El rotacional se puede interpretar como la cantidad de giro que provoca un campo, en este caso el campo de desplazamientos, alrededor del vector normal a la superficie en cada punto, en este caso el vector [math] \vec k [/math]. Es decir, está midiendo cuánto está girando cada cuadrado en los que hemos dividido el sólido, manteniendo siempre el área de cada uno. Tal y como se puede ver en los gráficos, el rotacional es función de [math] y [/math], y es mayor cuanto más aumenta la [math] y [/math]. Por lo tanto, los puntos que tienen mayor rotacional son los puntos cuya [math] y [/math] toma el valor 4.

4 Tensiones

En los apartados siguientes se resolverán todos los enunciados propuestos que están relacionados con tensiones en el sólido sobre el que estamos trabajando.

4.1 Tensiones normales

Primeramente, se define la parte simétrica del tensor gradiente de [math] \vec u [/math] como: [math] \epsilon (\vec u) = (\nabla \vec u + \nabla \vec u ^t)/2 [/math] siendo: [math] \nabla \vec u = \frac{\partial u_i}{\partial x^j} \vec e_i \otimes \vec e_j = \frac{y}{10} \vec i \otimes \vec j [/math] y su transpuesto: [math] \nabla \vec u^t = \frac{y}{10} \vec j \otimes \vec i [/math] Por lo tanto: [math] \epsilon (\vec u) = \frac{y}{20} \vec i \otimes \vec j + \frac{y}{20} \vec j \otimes \vec i [/math] Se define el tensor de tensiones a través de la fórmula: [math] \sigma = \lambda \nabla \cdot \vec u \mathbb{1} + 2 \mu \epsilon [/math] la cual tomando los valores [math] \lambda = \mu \ = 1 [/math] y particularizándola para el caso de nuestro problema es: [math] \sigma = 2 \epsilon = \frac{y}{10} \vec i \otimes \vec j + \frac{y}{10} \vec j \otimes \vec i [/math]

4.1.1 Normales a [math] \vec i[/math]

Las tensiones normales con respecto a [math] \vec i [/math] se definen como: [math] \vec i \cdot \sigma \cdot \vec i [/math] Operando con el tensor [math] \sigma [/math] definido anteriormente se llega al resultado de que: [math] \vec i \cdot \sigma \cdot \vec i = 0 [/math] Por lo tanto, al ser 0 para cualquier punto del sólido, no merece la pena realizar una representación gráfica. Esto nos dice que ningún punto del sólido sufre tensiones normales en la dirección que marca el eje [math] \vec i [/math]. Esta gráfica sería igual que la del módulo de la divergencia, que también toma el valor 0 para cualquier punto del sólido.

4.1.2 Normales a [math] \vec j[/math]

Las tensiones normales con respecto a [math] \vec j [/math] se definen, al igual que en apartado interior con el vector [math] \vec i [/math], como: [math] \vec j \cdot \sigma \cdot \vec j [/math] Operando igualmente con el tensor [math] \sigma [/math] definido anteriormente se llega al resultado de que: [math] \vec j \cdot \sigma \cdot \vec j = 0 [/math] Por lo tanto, al ser 0 para cualquier punto del sólido, tampoco merece la pena realizar una representación gráfica. Esto nos dice que ningún punto del sólido sufre tensiones normales en la dirección que marca el eje [math] \vec j [/math]. Esta gráfica sería, de nuevo, igual que la del módulo de la divergencia, que también toma el valor 0 para cualquier punto del sólido.

4.2 Tensiones tangenciales

4.2.1 Tangenciales a [math] \vec i[/math]

Representación de las tensiones tangenciales respecto al plano ortogonal a i

Las tensiones tangenciales respecto al plano ortogonal a [math] \vec i [/math] vienen definidas por la ecuación [math] |\sigma \cdot \vec i - (\vec i \cdot \sigma \cdot \vec i) \vec i| [/math]. Como hemos comprobado, [math] \vec i \cdot \sigma \cdot \vec i [/math] es nulo y calculando obtenemos que [math] \sigma \cdot \vec i = \frac{y}{10} \vec j [/math]. Escribimos el código en MatLab u OCTAVE UPM para obtener su representación y dar posteriormente una interpretación.

% Intervalos de trabajo
x=[-1/2:0.1:1/2];
y=[0:0.1:4];
% Mallado y z=0
[Mx,My]=meshgrid(x,y);
% Definición del rotacional
sigmai=My/10;
% Dibujo de los gráficos
subplot(2,1,1)
mesh(Mx,My,sigmai)
subplot(2,1,2)
surf(Mx,My,sigmai)
% Las tensiones tangenciales con respecto al plano ortogonal a i son iguales al rotacional

Como vemos, la gráfica de las tensiones tangenciales con respecto al plano ortogonal a [math] \vec i [/math] es igual que la del rotacional. Recordamos que la divergencia era nula, y por tanto distinta de las tensiones tangenciales.

4.2.2 Tangenciales a [math] \vec j [/math]

Representación de las tensiones tangenciales respecto al plano ortogonal a j

Las tensiones tangenciales respecto al plano ortogonal a [math] \vec j [/math] vienen definidas por la ecuación [math] |\sigma \cdot \vec j - (\vec j \cdot \sigma \cdot \vec j) \vec j| [/math]. Como hemos comprobado, [math] \vec j \cdot \sigma \cdot \vec j [/math] es nulo y calculando obtenemos que [math] \sigma \cdot \vec j = \frac{y}{10} \vec i [/math] . Escribimos el código en MatLab u OCTAVE UPM para obtener su representación y dar una interpretación.

% Intervalos de trabajo
x=[-1/2:0.1:1/2];
y=[0:0.1:4];
% Mallado y z=0
[Mx,My]=meshgrid(x,y);
% Definición del rotacional
sigmaj=My/10;
% Dibujo de los gráficos
subplot(2,1,1)
mesh(Mx,My,sigmaj)
subplot(2,1,2)
surf(Mx,My,sigmaj)
% Las tensiones tangenciales con respecto al plano ortogonal a j son iguales al rotacional

Como vemos, la gráfica de las tensiones tangenciales con respecto al plano ortogonal a [math] \vec j [/math] es igual que la del rotacional. Recordamos que la divergencia era nula, y por tanto distinta de las tensiones tangenciales.

4.3 Tensión de Von Mises

La tensión de Von Mises, magnitud escalar utilizada para indicar cuándo un material deja de comportarse como elástico puro y comienza su comportamiento plástico, viene dada por la fórmula: [math] \sigma_{VM} = \sqrt { \frac { \left( \displaystyle\ \sigma_1 - \sigma_2 \right)^2 + \left( \displaystyle\ \sigma_2 - \sigma_3 \right)^2 + \left( \displaystyle\ \sigma_3 - \sigma_1 \right)^2 } {2} } [/math] donde [math] \sigma_1 [/math], [math] \sigma_2 [/math] y [math] \sigma_3 [/math] son los autovalores de [math] \sigma [/math], también conocidos como tensiones principales. Para calcularlos, lo hacemos desde la matriz del tensor de tensiones dado que es la siguiente: [math] \begin{pmatrix} 0 & \frac{y}{10} & 0 \\ \frac{y}{10} & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} [/math] A partir de ella operamos hasta obtener los autovalores de la siguiente forma: [math] \begin{vmatrix} -\lambda & \frac{y}{10} & 0 \\ \frac{y}{10} & -\lambda & 0 \\ 0 & 0 & -\lambda \end{vmatrix} = -\lambda^3 + \frac{\lambda y^2}{100} = 0 [/math] De donde obtenemos los autovalores [math] \lambda_1 = 0 [/math], [math] \lambda_2 = \frac{y}{10} [/math] y [math] \lambda_3 = - \frac{y}{10} [/math]. Estos autovalores, aparte de obtenerlos analíticamente como se muestra, los obtenemos también con la función eig.m de MatLab, tal y como aparece en el código posterior, y como se pide en el enunciado. Una vez obtenidos los autovalores, la tensión de Von Mises sobre la placa queda representada así:

Tensión de Von Mises sobre la placa, máxima en y=4

% Intervalos de trabajo
x=[-1/2:0.1:1/2];
y=[0:0.1:4];
% Mallado y z=0
[Mx,My]=meshgrid(x,y);
[m,n]=size(Mx);
sigma=[0,1/10,0;1/10,0,0;0,0,0];
% Cálculo de los autovalores mediante el uso de MatLab
autovalores=eig(sigma);
t1=My.*autovalores(1);
t2=My.*autovalores(2);
t3=My.*autovalores(3);
% Tensión de Von Mises
tension=sqrt(((t1-t2).^2+(t2-t3).^2+(t3-t1).^2)/2);
%Dibujo de los gráficos
subplot(2,1,1)
mesh(Mx,My,tension)
subplot(2,1,2)
surf(Mx,My,tension)

Como podemos observar en los gráficos, la tensión de Von Mises alcanza su valor máximo en la placa en los puntos cuya componente en [math] y [/math] es igual a 4.

5 Masa

Dada la función de densidad siguiente: [math] d(x,y,z)=xye^{- \frac{1}{x^2}} [/math] Para calcular la masa de la placa debemos integrar la función de densidad respecto a [math]dx[/math] y [math]dy[/math] en toda la superficie de la placa. Pero debido a que la primitiva de esta función de densidad no se puede calcular por los métodos convencionales de integración, utilizaremos métodos numéricos con ayuda de MatLab u OCTAVE UPM.

5.1 Método de Simpson

Método de Simpson

El método de Simpson es un método numérico que se utiliza para el cálculo aproximado de integrales definidas simples. Se basa en la división del intervalo de integración en subintervalos, y aproximando el área debajo de la curva de cada dos subintervalos (siendo esta la definición de integral), por el área de una parábola que pasa por los extremos de la unión de esos dos subintervalos y por el centro del mismo. La fórmula siguiente es la expresión matemática de lo expuesto anteriormente: [math] \int_{a}^{b} f(x) \, dx \approx \frac{b-a}{6}\left[f(a) + 4f\left(\frac{a+b}{2}\right)+f(b)\right][/math] En lo que refiere a su implementación en Matlab u Octave UPM, lo más comodo es crear una función que realice el método de Simpson. De este modo, cada vez que necesitemos calcular una integral definida no sea necesario escribir el código de la función de nuevo, sino que simplemente con hacer una llamada a la función obtendremos su valor. El código MatLab que define la regla de Simpson es el siguiente:

% 'li' es el límite por la izquierda, 'ld' el límite por la derecha, y 'n' el númmero de
% intervalos
function r=simpson(f,li,ld,n)
% Definición de la constante que multiplica el polinomio y del vector de subintervalos
h=(ld-li)/(2*n);
v=linspace(li,ld,n+1);
% Bucle for para sumar todas las áreas que se van obteniendo de las aproximacionse de
% cada intervalo y fórmula de la regla de Simpson
r=0;
for i=1:n
   r=r+(feval(f,v(i))+feval(f,v(i+1))+4*feval(f,(v(i)+h)))*h/3;
end

Para poder usar esta función, habrá que definir otra función Matlab con la función matemática que deseamos integrar, en este caso la densidad definida anteriormente. El código MatLab u OCTAVE UPM es el siguiente:

function d=f(x)
% Fórmula de la densidad proporcionada en el enunciado(después de haber integrado con respecto a x)
d=x*exp(-1/(x^2));
end

Por último, escribimos un programa para el cálculo de la integral definida haciendo uso de las dos funciones. Como la regla de Simpson sirve para el cálculo de integrales simples y esta integral es doble, primero hacemos el cálculo de la integral inmediata con respecto a [math] y [/math] (cuyo valor da 8) y luego usamos el método para la integral con respecto a [math] x [/math]. Como esta última da 0, pero sin embargo la masa no puede tener ese valor, se observa que la función presenta una simetría axial, y que por lo tanto la masa será dos veces la integral desde 0 hasta 0,5. El código Matlab es el siguiente:

% El 10000 es el número de intervalos elegidos
int=simpson('f',0,0.5,10000);
masa=8*2*int

La masa del sólido tratado es por lo tanto: [math] masa = 0.0063965 [/math]

5.2 Otro método numérico

Este es otro método numérico para calcular la masa de la placa algo menos exacto que el de Simpson. Este método se basa en calcular la masa de un [math] dm [/math] y sumar todos estos diferenciales hasta obtener la masa completa de la placa. El método es el siguiente:

% Lado de los cuadrados en los que dividimos el sólido
h=1/1000;
% Creación del mallado del sólido
x=-0.5:h:0.5;
y=0:h:4;
[Mx,My]=meshgrid(x,y);
% Valor de la densidad para cada punto del mallado
d=Mx.*My.*exp(-1./(Mx.^2));
% Masa de cada cuadrado que conforma el mallado
dm=abs(d*h^2);
% Suma de las masas de todos los cuadrados
masa=sum(sum(dm))

La masa del sólido es por lo tanto: [math] masa = 0,0064716 [/math] En consecuencia, se comprueba que los valores obtenidos de la masa por ambos métodos son similares, con lo que ambas aproximaciones se pueden considerar aceptables.