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

Trabajo realizado por estudiantes
Título	Visualización de campos escalares y vectoriales en elasticidad. Grupo 5-C
Asignatura	Teoría de Campos
Curso	2014-15
Autores	Manuel Morales, Jorge Villa, Sergio Rodríguez, Diego Pontiveros, Manuel Jugo, Lourdes Sánchez-Ocaña
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Enunciado

Consideramos una placa plana (en dimensión 2) que ocupa la región comprendida entre las parábolas [math]P_1=18y-81x^2-1=0[/math] y [math]P_2=2y+x^2-1=0[/math].Consideramos una placa para representarla usaremos un sistema de coordenadas adaptado a la geometría que nos dan: \begin{cases} x=u·v; \\ y=1/2(u^2-v^2); \end{cases}

con u y v definidas en (u, v) ∈ [1/3, 1] × [−1, 1].

En ella vamos a suponer que tenemos definidas dos cantidades físicas: la temperatura [math]T(x,y)[/math] y los desplazamientos [math]\vec u(u,v)[/math].De esta forma, si definimos [math]r_0(u,v)[/math] el vector de posición de los puntos de la placa en reposo, la posición de cada punto [math](u,v)[/math] de la placa viene dada por: [math] \vec r (u,v)= \vec r_{0}(u,v)+\vec u(u,v). [/math] Vamos a suponer que sobre la placa se ha aplicado una fuerza que ha provocado una vibración de manera que los desplazamientos vienen dados por: [math]\vec u(u,v)=\vec a · (\vec b · \vec r_0). [/math] donde [math]\vec a[/math] y [math]\vec b[/math] son vectores dados.

En este trabajo supondremos lo siguiente:

1. [math]\vec a=\frac{\vec g_u}{|\vec g_u|}[/math] ;
2. [math]\vec b=-4\frac{\vec g_u}{|\vec g_u|}[/math].

2 Mallado de los puntos interiores del sólido

Para la representación del mallado, hemos utilizado el comando "mesh", que nos permite definir los límites de la placa en las variables x e y, utilizando un muestreo de 1/20. La malla se creará mediante una matriz en la que al última componente sea cero, ya que estamos hablando de una malla plana.

Mallado de los puntos interiores del sólido

h=1/20;
u=1/3:h:1;
v=-1:h:1;
[uu,vv]=meshgrid(u,v);
xx=uu.*vv;
yy=0.5*(uu.^2-vv.^2);
mesh(xx,yy,0*xx)
view(2)

2.1 Cálculo de la base natural

Es necesario definir un nuevo sistema de coordenadas curvilíneas {0;(u,v,w)} respecto a una base natural, compuesta por tres vectores regulares ortogonales entre sí, que no son necesariamente unitarios. Las componentes de la base natural {[math]\overrightarrow{g_u},\overrightarrow{g_v},\overrightarrow{g_w}[/math]} son las derivadas parciales del vector posición, [math]\overrightarrow{r}=x\overrightarrow{i}+y\overrightarrow{j}+z\overrightarrow{k}[/math] respecto de las coordenadas curvilíneas [math](u,v)[/math].

Al realizar el cambio de variable,

\begin{cases} x=uv\\ y=\frac{(u^2−v^2)}{2}\\ z=0 \end{cases}

calculamos los vectores de la base natural.

[math]\overrightarrow{g_u} = \frac{∂r}{∂u} = v\overrightarrow{i} + u\overrightarrow{j}[/math]

[math]\overrightarrow{g_v} = \frac{∂r}{∂v} = u\overrightarrow{i} - v\overrightarrow{j}[/math]

[math]\overrightarrow{g_w}= 0[/math]

Siendo [math]r(x, y, z)=x\overrightarrow{i} + y\overrightarrow{j} + z\overrightarrow{k}[/math]

3 Influencia de un foco de calor

La temperatura viene dada por la función [math]T(x,y)=(8-y^2+2y)e^(-x^2)[/math]

clear all

h=1/20;
u=1/3:h:1;
v=-1:h:1;
[uu,vv]=meshgrid(u,v);
xx=uu.*vv;
yy=0.5*(uu.^2-vv.^2);
figure(1)
mesh(xx,yy,0*xx)


T=(8-yy.^2+2*yy).*exp(-(xx).^2);
A=(8-yy.^2+2*yy).*(-2*xx).*exp(-(xx).^2);
B=exp(-(xx).^2).*(-2*yy+2);
figure(2)
contour(xx,yy,T,20)
figure(3)
quiver(xx,yy,A,B)

figure(4)
quiver(xx,yy,A,B)
hold on
contour(xx,yy,T,20,'r') %dibujar con 10 curvas de nivel
hold off

Curvas de nivel de la temperatura

Vector gradiente de la temperatura

Gradiente y curvas de nivel de la temperatura superpuestos

Se observa en los gráficos que la temperatura será mayor en el centro de la placa debido a la distribución dada en el enunciado.

4 Campo de desplazamiento

En este apartado, veremos, gracias a una gráfica que se produce antes y después del desplazamiento.

clear all

h=1/20;
u=1/3:h:1;
v=-1:h:1;
[uu,vv]=meshgrid(u,v);
xx=uu.*vv;
yy=0.5*(uu.^2-vv.^2);
figure(1)
subplot(1,2,1)
mesh(xx,yy,0*xx)
view(2)
%axis([1/3,1,-1,1])


fx=inline('-2*xx.*yy','xx','yy'); 
fy=inline('-2*(xx.^2)','xx','yy');                
U=fx(uu,vv);                             
V=fy(uu,vv);  


A=xx+U;
B=yy+V;
subplot(1,2,2)
mesh(A,B,0*A)
view(2)
figure(2) 
quiver(xx,yy,U,V); 
view(2)

Campo de vectores desplazamiento sobre el mallado

Sólido antes y después de aplicar el vector ū

Puede apreciarse en las gráficas cómo el campo vectorial tiene un efecto mayor en los puntos exteriores de la malla, ya que el campo vectorial tiene un mayor módulo.

5 Divergencia

La divergencia es una muestra del incremento de volumen que se da en un punto cualquiera debido a la acción de una fuerza sobre el mismo. La variación de volumen sufrida en ese punto (a consecuencia de la fuerza) 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 cree un incremento de volumen que la de un campo que simplemente desplace la superficie.

h=1/20;
u=1/3:h:1;
v=-1:h:1;
[uu,vv]=meshgrid(u,v);
xx=uu.*vv;
yy=0.5*(uu.^2-vv.^2);
figure(1)
mesh(xx,yy,0*xx)
view(2)

div=-2*(3*uu.^2+vv.^2)./(uu.^2+vv.^2);
figure(2)
surf(xx,yy,div)
view(2)

Representación gráfica de la divergencia

Gráfica 3-D de la divergencia

Comparando la gráfica del desplazamiento entre el antes y el después, se puede apreciar una mayor variación de volumen local en los puntos que se encuentren a una mayor distancia de los puntos del centro y cuando su ordenada es menor. En la gráfica de la divergencia se puede observar como es ésta máximo en los puntos anteriormente descritos.

6 Rotacional

El rotacional es un operador vectorial que indica la velocidad y dirección de giro que produce un campo.

Se define el rotacional de [math]\overrightarrow{u}[/math] como [math]\nabla \times \overrightarrow{u}. [/math] Dado un campo [math]\overrightarrow{u}=u_i\overrightarrow{g^i}[/math] se tiene que [math] \nabla\times\overrightarrow{u}=\frac{ 1}{ √g} \begin{vmatrix} \overrightarrow{g_u} & \overrightarrow{g_v} & \overrightarrow{g_w} \\ \frac{\partial}{\partial u} & \frac{\partial}{\partial v} & \frac{\partial}{\partial w} \\ u_u & u_v & u_w \end{vmatrix}[/math] donde [math]g=det(G)[/math] y [math]G[/math] es la matriz de Gram de la base natural.

En este caso se define el campo a partir de los vectores [math] \vec{a}= \frac{\vec{g}_u}{|\vec{g_u}|}[/math] y [math] \vec{b}=-4 \frac{\vec{g}_u}{|\vec{g_u}|}[/math] de forma que el campo es [math] \vec {u}(u,v)= \vec{a} (\vec{b} \vec{r_o}) [/math], sabiendo que [math]\overrightarrow{r_o}[/math] es el vector posición.

Al calcular la base natural se demuestra el primer vector de la base natural es [math]\vec{g_u}=v \hat{e_1} +u \hat{e_2} [/math]. A partir de este y el vector posición se calcula el campo de desplazamientos: [math]\vec{u}= \frac{\vec{g}_u}{|\vec{g_u}|}(-4 \frac{\vec{g}_u}{|\vec{g_u}|} \vec{r_o})= \frac{\vec{g}_u}{|\vec{g_u}|^2}(-4 \vec{g_u} \vec{r_o})=\frac{ \vec{g_u}}{u^2+v^2}(-4(uv^2+\frac{u}{2}(u^2-v^2)))=-4(\frac{uv^2+u^3}{2})\frac{\vec{g_u}}{u^2+v^2}=-2u\vec{g_u}=-2uv \hat{e_1} -2u^2 \hat{e_2}[/math]

En este caso se demuestra que el rotacional es nulo [math]\nabla \times\vec{u}= \frac{1}{ \sqrt{g} } \begin{bmatrix} g_{u} & g_{v} & g_{w} \\ \frac{\partial}{\partial u} & \frac{\partial}{\partial v} & \frac{\partial}{\partial w} \\ u_{u} &u_{v}&u_{w}\end{bmatrix} = \frac{1}{ u^2+v^2 } \begin{bmatrix} g_{u} & g_{v} & g_{w} \\ \frac{\partial}{\partial u} & \frac{\partial}{\partial v} & \frac{\partial}{\partial w} \\ -2u(u^2+v^2) &0&0\end{bmatrix}=\frac{u^2+v^2}{ u^2+v^2 } \begin{bmatrix} g_{u} & g_{v} & g_{w} \\ \frac{\partial}{\partial u} & \frac{\partial}{\partial v} & \frac{\partial}{\partial w} \\-2u&0&0\end{bmatrix}= -2u ( \frac{\partial g_{v} }{\partial w} - \frac{\partial g_{w} }{\partial v} )=0[/math]

7 Tensiones en la dirección de los vectores de la base natural.

Calculados la divergencia y el rotacional de [math]\overrightarrow{u}[/math], pasamos a definir el tensor de tensiones:

[math] \sigma= \lambda \nabla\cdot\overrightarrow{u}1 +2\mu\epsilon [/math]
Se parte de las hipótesis de medio elástico lineal, isótropo y homogéneo. Los valores [math]\lambda[/math] y [math]\mu[/math] se conocen como los Coeficientes de Lamé, dependen de las propiedades elásticas de cada material y se supone su valor la unidad. Por último [math]\epsilon[/math] viene dado por la ecuación:

[math]\epsilon=(\nabla\overrightarrow{u}+\nabla\overrightarrow{u}^t)/2 [/math]
El siguiente paso es realizar un gráfico que nos muestre el tensor de tensiones normales en las direcciones de la base natural, [math] \overrightarrow{g_u} [/math] y [math] \overrightarrow{g_v} [/math]. Para visualizarlo hacemos uso de las ecuaciones e función de la base natural: [math] \frac{\overrightarrow{g_u}}{|\overrightarrow{g_u}|}\sigma\frac{\overrightarrow{g_u}}{|\overrightarrow{g_u}|} , \frac{\overrightarrow{g_v}}{|\overrightarrow{g_v}|}\sigma\frac{\overrightarrow{g_v}}{|\overrightarrow{g_v}|} [/math]

Hemos vuelto a calcular [math] \overrightarrow{u} [/math], pero esta vez dejando que Matlab realizara los cálculos. Las matrices [math] \nabla\overrightarrow{u}[/math] y [math]\nabla\overrightarrow{u}^t[/math] y la matriz del tensor de tensiones calculadas analíticamente son:
[math] \nabla\overrightarrow{u}=\begin{pmatrix}-6u^2-2v^2 & -4uv & 0\\0 & 0 & 0\\0 & 0 & 0\end{pmatrix} \nabla\overrightarrow{u}^t=\begin {pmatrix}-6u^2-2v^2 & 0 & 0 \\-4uv & 0 & 0\\0 & 0 & 0\end{pmatrix}[/math]

[math] \sigma= \begin{pmatrix}\ \frac{-6(3u^2+v^2)}{u^2+v^2}-12u^2-4v^2&-4uv&0\\ -4uv & \frac{-6(3u^2+v^2)}{u^2+v^2} & 0\\0 & 0 & \frac{-6(3u^2+v^2)}{u^2+v^2} \end{pmatrix} [/math]

Y ya por último, con la fórmula anteriormente expuesta, calcularemos el tensor de tensiones:

[math] \sigma ^ij=\frac{-6u^2-2v^2}{u^2+v^2} \begin{bmatrix}1&0&0 \\0&1&0\\0&0&1 \end{bmatrix} +2 \begin{bmatrix} -6u^2-2v^2 & -2uv & 0 \\-2uv & 0 & 0\\0 & 0 & 0\end{bmatrix}=\begin{bmatrix}\frac{-6(3u^2+v^2)}{u^2+v^2}-12u^2-4v^2&-4uv&0\\ -4uv & \frac{-6(3u^2+v^2)}{u^2+v^2} & 0\\0 & 0 & \frac{-6(3u^2+v^2)}{u^2+v^2} \end{bmatrix} [/math] El calculo de las tensiones en una dirección se calcula mediante un producto tensorial que proyecta las tensiones sobre la dirección indicada.

clear all
h=1/20;
u=1/3:h:1;
v=-1:h:1;
[uu,vv]=meshgrid(u,v);
xx=uu.*vv;
yy=0.5*(uu.^2-vv.^2);
figure(1)
mesh(xx,yy,0*xx)
Gu=-(vv.*((4.*uu.^2.*vv)./(uu.^2+vv.^2).^(1/2)+(vv.*((6.*uu.^2+2.*vv.^2)./(uu.^2+vv.^2)+12.*uu.^2+4.*vv.^2))./(uu.^2+vv.^2).^(1/2)))./(uu.^2+vv.^2).^(1/2)-(uu.*((uu.*(6.*uu.^2+2.*vv.^2))./(uu.^2+vv.^2).^(3/2)+(4.*uu.*vv.^2)./(uu.^2+vv.^2).^(1/2)))./(uu.^2+vv.^2).^(1/2);
Gv=(uu.*((4.*uu.*vv.^2)./(uu.^2+vv.^2).^(1/2)-(uu.*((6.*uu.^2+2.*vv.^2)./(uu.^2+vv.^2)+12.*uu.^2+4.*vv.^2))./(uu.^2+vv.^2).^(1/2))./(uu.^2+vv.^2).^(1/2)-(vv.*((vv.*(6.*uu.^2+2.*vv.^2))./(uu.^2+vv.^2).^(3/2)-(4.*uu.^2.*vv)./(uu.^2+vv.^2).^(1/2)))./(uu.^2+vv.^2).^(1/2));
figure(2)
mesh(xx,yy,Gu)
view(2)

figure(3)
mesh(xx,yy,Gv)
view(2)

Tensión en la dirección g_u

Tensión en la dirección g_v

Divergencia del vector desplazamiento

El módulo del rotacional del vector desplazamiento es 0, por lo que no se ha dibujado la gráfica. Se puede apreciar que las tensiones en la dirección [math]\overrightarrow{g_v}[/math], son similares a las de la divergencia, mientras que las de la dirección [math]\overrightarrow{g_u}[/math] no tienen una clara relación.

8 Tensión de Von Mises

La tensión de Von Mises (o el Esfuerzo) es un índice obtenido de la combinación de los Esfuerzos Principales en un momento dado para determinar en qué puntos ocurre el esfuerzo en el eje X, Y y Z y provoca el fallo. Este método de cálculo se utiliza para medir el esfuerzo y las distribuciones de tensión dentro de un material dúctil. Para el calculo es necesario calcular los autovalores de la matriz de tensiones en cada punto del sólido.

La matriz del tensor de tensiones se define como [math] \sigma [/math], puede calcularse fácilmente a partir de las tensiones principales (autovalores) del tensor tensión en un punto del sólido deformable, mediante la expresión,
[math] \sqrt{\frac{(\sigma_1-\sigma_2)^2-(\sigma_2-\sigma_3)^2-(\sigma_3-\sigma_1)^2}{2}} [/math]

%vonmises
clear all
h=1/20;
u=1/3:h:1;
v=-1:h:1;
[uu,vv]=meshgrid(u,v);
xx=uu.*vv;
yy=0.5*(uu.^2-vv.^2);
mesh(xx,yy,0*xx);
%syms xx yy
lambdatotal=zeros(41,14);
for i=1:41
    for j=1:14
        x=uu(i,j);
        y=vv(i,j);
        A=[-12*x.^2-4*y.^2+(-6*x.^2-2*y.^2)./(x.^2+y.^2) (-4*x.*y) 0;-4*x.*y (-6*x.^2-2*y.^2)./(x.^2+y.^2) 0;0 0 (-6*x.^2-2*y.^2)./(x.^2+y.^2)];
        lambda=eig(A);
        lambdavm=sqrt((((lambda(1)-lambda(2))^2+(lambda(2)-lambda(3))^2+(lambda(3)-lambda(1))^2)/2));
        lambdatotal(i,j)=lambdavm;
    end
end
surf(xx,yy,lambdatotal)
colorbar

Representación 3-D de la tensión de Von Mises

De la gráfica se observa que las tensiones son mayores en los extremos que en el centro y también mayores cuando aumenta su ordenada. Se puede explicar este fenómeno ya que es en los extremos donde se produce una mayor deformación que hace necesaria una mayor tensión en dichos puntos.

9 Masa total del sólido

Por último se calculará la masa de la placa. Se da como dato la función densidad f(x,y) y la malla hallada anteriormente. Por lo tanto la masa queda definida mediante la integral doble de la función en el dominio x perteneciente a [-0.9833,0.9833] e y a [-0.4444,0.4835]
[math] \iint f(x,y)dxdy [/math]

%tomaremos los datos necesarios de apartados anteriores
h1=1/20;
f=xx.*yy.*exp(-1./xx.^2);
c=(h1.^2).*f;
masa=(sum(sum(c)))

El valor de la masa de la placa es de 1.9111e-020. Sim embargo, si utilizamos la funcion de matlab quad2d y aplicamos los límites de integración para x e y, nos sale un valor de la masa -4.8789e-019. Por tanto, decidimos hacer la integral doble de la función densidad en x e y, y nos sale que la masa es nula.

[math] \iint x·y·exp(-1/x^2)dxdy=0 [/math]