Visualización de campos escalares y vectoriales en elasticidad (G4-A)
Contenido
1 Visualización de la placa.
El objetivo de este estudio es la visualización de campos escalares y vectoriales en elasticidad.
Inicialmente consideramos una placa plana que ocupa la región comprendida entre dos parábolas
- P1: 18y − 81x2 −1 = 0
- P2: 2y + x2 −1 = 0.
Para representarla usaremos un sistema de coordenadas adaptado a la geometría que nos dan:
- [math] x=uv [/math]
- [math]y=\frac{(u^2−v^2)}{2}[/math]
sabiendo que las nuevas coordenadas quedan definidas en el intervalo (u,v) ∈ [⅓,1] x [-1,1].
Para representar la placa dibujaremos un mallado que represente los puntos interiores del sólido tomando como
paso de muestreo h = 1/20 para las coordenadas u y v. Ejes en el intervalo (x,y) ∈ [-1; 1] x [-1; 1].
Las líneas de coordenadas están representadas en la figura.
h=0.05;
u=1/3:h:1;
v=-1:h:1;
[uu,vv]=meshgrid(u,v);
xx=uu.*vv;
yy=1/2.*(uu.^2-vv.^2);
mesh(xx,yy,0*xx)
axis([-1,1,-1,1])
xlabel('Eje X')
ylabel('Eje Y')
1.1 Cálculo de la base natural
Al definir un nuevo sistema de coordenadas curvilíneas {0;(u,v,w)} hemos de definir a su vez una base natural que defina dicho espacio vectorial en el que se representan. A esa base la notaremos como base natural. Se trata de tres vectores(en nuestro estudio serán dos ya que la tercera componente de la base es nula) regulares(ortogonales entre sí), que no sean 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] en función de las coordenadas curvilíneas [math](u,v)[/math].
Teniendo en cuenta el cambio de variable:
- [math] x=uv [/math]
- [math]y=\frac{(u^2−v^2)}{2}[/math]
- [math] z=0 [/math]
calculamos los vectores de la base natural.
[math]\overrightarrow{g_u} = v\overrightarrow{i} + u\overrightarrow{j}[/math]
[math]\overrightarrow{g_v} = u\overrightarrow{i} - v\overrightarrow{j}[/math]
[math]\overrightarrow{g_w}= 0[/math]
Así mismo vamos a definir sobre dicha placa dos cantidades físicas: la temperatura T(u,v) que depende de las dos coordenadas curvilíneas (u,v), y los deplazamientos ū(u,v)producidos por la acción de una fuerza determinada.
2 Temperatura de la placa.
Una de las medidas que utilizaremos para el estudio de esta placa será cómo varía la temperatura en su superficie. La temperatura de cualquier material suele estar muy ligada con la elasticidad del propio material de esta forma:[math] \Delta L=\alpha\Delta T[/math]
La temperatura en la placa viene dada por el campo escalar
- [math]T(x,y)=(8-y^2+2y)e^{-(x)^2} [/math]
h=0.05;
u=1/3:h:1;
v=-1:h:1;
[uu,vv]=meshgrid(u,v);
xx= uu.*vv;
yy=1/2.*(uu.^2-vv.^2);
T=(8-yy.^2+2.*yy).*exp(-(xx.^2));
contour(xx,yy,T,30)
colorbar
max(max(T))
xlabel('Eje X');
ylabel('Eje Y');
axis([-1.1,1.1,-0.6,0.6])
La temperatura es una función de x e y, pero al hacer el cambio de variable pasará depender de las variables u y v .
Para saber en que punto de la superficie la temperatura es mas alta tenemos que calcular el máximo de dicha función. Para esto tenemos que hallar sus derivadas parciales, igualarlas a 0 y resolver el sistema.
Nos encontramos un máximo absoluto en (x,y) = (0,1) que esta fuera del dominio, pero como nuestra función está acotada podemos afirmar que el valor máximo se encuentra en el punto de la frontera más cercano al máximo absoluto de la función T, el punto (x,y) = (0,½).
2.1 Estudio del gradiente de temperaturas.
Por lo tanto, y comenzando una estela de gráficas útiles para sacar conclusiones sobre el comportamiento de la placa, estudiaremos ahora la variación de la Temperatura dada por el gradiente de esta.
El [math]\nabla T[/math] es la variación de temperatura por unidad de distancia y lo obtenemos derivando la función de la temperatura T respecto de sus variables (x e y):
[math]\nabla T= \frac{\partial {T_i}}{\partial x^i }= \begin{pmatrix} \frac{\partial{T_1}}{\partial x} & \frac{\partial{T_2}}{\partial y} \\ \frac{\partial{T_1}}{\partial y} & \frac{\partial{T_2}}{\partial x} \end{pmatrix} [/math]
Como podemos observar en la figura, el campo vectorial [math]\nabla T[/math] es ortogonal a las curvas de nivel de la temperatura.
h=0.05;
u=1/3:h:1;n=length(u);
v=-1:h:1;
[uu,vv]=meshgrid(u,v);
xx=uu.*vv;
yy=1/2.*(uu.^2-vv.^2);
T=(8-yy.^2+2.*yy).*exp(-(xx.^2));
dx=-2.*xx.*exp(-(xx).^2).*(8-yy.^2+2.*yy);%Derivada de la T respecto a x
dy=(-2.*yy+2).* exp(-(xx).^2); %Derivada de la T respecto a y
quiver(xx,yy,dx,dy)
hold on
contour(xx,yy,T,30)
hold off
xlabel('Eje X')
ylabel('Eje Y')
3 Efectos del campo de desplazamientos [math]\overrightarrow{u}[/math] sobre la placa.
Continuando el estudio de un sólido aplicaremos un campo vectorial a la superficie anteriormente hallada para observar el comportamiento del sólido.
Partiendo de [math]\overrightarrow{r_0(u,v)}[/math], vector de posición de los puntos de la placa en reposo, la posición de cada punto de la placa vendrá dado por:
[math]\overrightarrow{r(u,v)}=\overrightarrow{r_0(u,v)}+\overrightarrow{u(u,v)}[/math]
3.1 Campo de vectores en el mallado del sólido.
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 desplazmientos:
[math]\overrightarrow{u(x,y)}=\overrightarrow{a}(\overrightarrow{b}\cdot \overrightarrow{r_0})[/math].
Donde [math]\overrightarrow{a} [/math] y [math]\overrightarrow{b} [/math] son vectores dados y los cuales se definen para el estudio del comportamiento de esta placa, serán los siguientes:
[math]\overrightarrow{a}= \frac{\overrightarrow{g_u}}{|\overrightarrow{g_u}|} [/math]
[math]\overrightarrow{b}= -4\frac{\overrightarrow{g_u}}{|\overrightarrow{g_u}|}[/math]
Las componentes del vector ū son halladas mediante la fórmula antes dada.
El campo vectorial del desplazamiento ū está representado en la siguiente figura:
h=0.05;
u=1/3:h:1;n=length(u);
v=-1:h:1;
[uu,vv]=meshgrid(u,v);
aa=((-2*((uu.^3).*vv+(vv.^3).*uu))./(uu.^2+vv.^2))+uu.*vv;
bb=(((-2*(uu.^2).*vv.^2+uu.^4))./(uu.^2+vv.^2))+0.5*(uu.^2-vv.^2);
xx=uu.*vv;
yy=1/2.*(uu.^2-vv.^2);
quiver(xx,yy,aa,bb)
xlabel('Eje X')
ylabel('Eje Y')
3.2 Desplazamientos provocados en la placa.
Aplicando al vector r0 el desplazamiento ū obtenemos el vector r. Mediante la siguiente figura podemos ver el sólido antes del desplazamiento (a la izquierda), y después (a la derecha).
h=0.05;
u=1/3:h:1;n=length(u);
v=-1:h:1;
[uu,vv]=meshgrid(u,v);
xx=uu.*vv;
yy=1/2.*(uu.^2-vv.^2);
%Coordenadas i y j del vector r0
aa=((-2*((uu.^3).*vv+(vv.^3).*uu))./(uu.^2+vv.^2)) +uu.*vv;
bb=(((-2*(uu.^2).*vv.^2+uu.^4))./(uu.^2+vv.^2)) +0.5*(uu.^2-vv.^2);
%Coordenadas i y j del vector r.
subplot(1,2,1)
mesh(xx,yy,0*xx)
xlabel('Eje X'); ylabel('Eje Y')
subplot(1,2,2)
mesh(aa,bb,0*xx)
xlabel('Eje X'); ylabel('Eje Y')
3.3 Estudio de la divergencia de [math]\overrightarrow{u}[/math].
Denotaremos la divergencia de [math]\overrightarrow{u}[/math] como [math]\nabla\cdot\overrightarrow{u} [/math]. Dado un campo [math]\overrightarrow{u}=u^i\overrightarrow{g_i}[/math] se tiene que [math] \nabla\cdot\overrightarrow{u}=\frac{ 1}{ √g}·\frac{\partial }{\partial x^i}\cdot(√g\cdot u^i) [/math] donde [math]g=det(G)[/math] y [math]G[/math] es la matriz de Gram de la base natural {[math]\overrightarrow{g_u},\overrightarrow{g_v},\overrightarrow{g_w}[/math]}.
Nuestro vector desplazamiento en la base natural es: [math]\overrightarrow{u}=u^u\overrightarrow{g_u}+u^v\overrightarrow{g_v}[/math].
Para calcular las componentes contravariantes lo hacemos a partir de [math]u^i=\overrightarrow{u}\overrightarrow{g^i}[/math]. Entonces, operando nos queda:
[math]u^u=\overrightarrow{u}\overrightarrow{g^u}= (\frac{-4uv^2-2u(u^2-v^2)}{u^2+v^2}(v\overrightarrow{i}+u\overrightarrow{j}))\cdot(\frac{1}{u^2+v^2}(v\overrightarrow{i}+u\overrightarrow{j})= \frac{-4uv^2-2u(u^2-v^2)}{u^2+v^2})[/math]
[math]u^v=\overrightarrow{v}\overrightarrow{g^v}=(\frac{-4uv^2-2u(u^2-v^2)}{u^2+v^2}(v\overrightarrow{i}+u\overrightarrow{j}))\cdot(\frac{1}{u^2+v^2}(u\overrightarrow{i}-vu\overrightarrow{j}))= 0[/math]
Con esto ya podemos calcular la divergencia, que quedará
[math] \nabla\cdot\overrightarrow{u}=\frac{1}{√(u^2+v^2)^2}\cdot[\frac{\partial}{\partial u}\cdot (\frac{-4uv^2-2u(u^2-v^2)}{(u^2+v^2)}\cdot √(u^2+v^2)^2)+0]=\frac{-2(3u^3+v^2)}{u^2+v^2} [/math]
Para dibujar la divergencia de [math]\overrightarrow{u}[/math] en todos los puntos del sólido utilizamos el siguiente código MATLAB:
clear all
h=0.05;
u=1/3:h:1;
v=-1:h:1;
[uu,vv]=meshgrid(u,v);
divU= (-2*(3.*uu.^2+(vv.^2)))./(uu.^2+vv.^2);
xx=uu.*vv;
yy=1/2.*(uu.^2-vv.^2);
max(max(divU))
mesh(xx,yy,divU)
xlabel('Eje X'); ylabel('Eje Y')
La divergencia es una medida de variación de volumen local debido al desplazamiento. El resultado de divergencia obtenido es negativo, por lo tanto, sabemos que se trata de un sumidero, es decir, el flujo entrante es mayor al flujo saliente.
En la gráfica se representan las variaciones de volumen y observamos que los puntos de la placa que sufren una mayor cambio de volumen son aquellos situados en los extremos inferiores, correspondiendo a los puntos de mayor divergencia.
3.4 Estudio del rotacional de [math]\overrightarrow{u}[/math].
Denotaremos 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 {[math]\overrightarrow{g_u},\overrightarrow{g_v},\overrightarrow{g_w}.[/math]}
Nuestro vector desplazamiento en la base recíproca queda: [math]\overrightarrow{u}=u_u\overrightarrow{g^u}+u_v\overrightarrow{g^v}.[/math]
Para calcular las componentes convariantes lo hacemos a partir de [math]u_i=\overrightarrow{u}\overrightarrow{g_i}[/math]. Entonces, operando nos queda:
[math]u_u=\overrightarrow{u}\overrightarrow{g_u}=(\frac{-4uv^2-2u(u^2-v^2)}{u^2+v^2}(v\overrightarrow{i}+u\overrightarrow{j}))\cdot(v\overrightarrow{i}+u\overrightarrow{j})= -4u^2-2u^3+2uv^2[/math]
[math]u_v=\overrightarrow{u}\overrightarrow{g_v}= (\frac{-4uv^2-2u(u^2-v^2)}{u^2+v^2}(v\overrightarrow{i}+u\overrightarrow{j}))\cdot(u\overrightarrow{i}-v\overrightarrow{j})= 0[/math]
Entonces el rotacional es [math]\nabla\times\overrightarrow{u}=\frac{ 1}{ √(u^2+v^2)^2} \begin{vmatrix} \overrightarrow{g_u} & \overrightarrow{g_v} & \overrightarrow{g_w} \\ \frac{\partial}{\partial u} & \frac{\partial}{\partial v} & \frac{\partial}{\partial w} \\ -4u^2-2u^3+2uv^2 & 0 & 0 \end{vmatrix}=\frac{1}{u^2+v^2}(4uv\overrightarrow{g_w})=0 [/math]
El rotacional estudia la tendencia del campo a inducir rotación alrededor de un punto. Hemos obtenido que el rotacional es nulo, es decir, es irrotacional. La irrotacionalidad del sólido implica que la divergencia es distinta de cero, como hemos comprobado en el apartado anterior.
4 Tensor de tensiones
Una vez calculados el rotacional y el divergente de [math]\overrightarrow{u}[/math], suponemos realizar el estudio en un medio elástico lineal, isótropo y homogéneo en el cual podemos definir un tensor de tensiones que viene dado por la fórmula:
- [math] \sigma= \lambda \nabla\cdot\overrightarrow{u}1 +2\mu\epsilon [/math]
donde [math]\lambda[/math] y [math]\mu[/math], también llamados los Coeficientes de Lamé, tienen valor la unidad en nuestro caso pero dependen de las propiedades elásticas de cada material. Y [math]\epsilon[/math] viene dado por la ecuación:
[math]\epsilon=(\nabla\overrightarrow{u}+\nabla\overrightarrow{u}^t)/2 [/math]
Procederemos a realizar un gráfico que nos muestre el tensor de tensiones normales en la base natural, [math] \overrightarrow{g_u} [/math] y [math] \overrightarrow{g_v} [/math], suponiendo como se ha ido viendo a lo largo del estudio que la tercera componente de la base natural es inexistente. Para visualizar el tensor de tensiones normales hacemos uso de las ecuaciones:
[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], respectivamente de la base natural.
Hemos vuelto a calcular [math] \overrightarrow{u} [/math], pero esta vez dejando que Matlab realizara los cálculos(no tiene ninguna relevancia para este caso de nuestro estudio).
Las matrices [math] \nabla\overrightarrow{u}[/math] y [math]\nabla\overrightarrow{u}^t[/math]:
[math] \nabla\overrightarrow{u}=\begin{pmatrix}-6u^2-2v^2 & -4uv\\0 & 0\end{pmatrix} \nabla\overrightarrow{u}^t=\begin
{pmatrix}-6u^2-2v^2 &-2uv \\-2uv & 0\end{pmatrix}[/math]
La matriz de tensor de tensores realizada analíticamente:
[math] \sigma= \begin{pmatrix}\ \frac{-6(3u^2+v^2)}{u^2+v^2}&-6uv\\ -2uv & \ \frac{-2(3u^2+v^2)}{u^2+v^2} \end{pmatrix} [/math]
clc
clear all
h=1/20;
u=1/3:h:1;
v=-1:h:1;
[uu,vv]=meshgrid(u,v);
xx=uu.*vv;
yy=1/2*(uu.^2-vv.^2);
r0=[uu.*vv,1/2*(uu.^2-vv.^2)];
gu=[vv,uu];
gv=[uu,-vv];
a=gu./(abs(gu));
b=-4*(gu./(abs(gu)));
uvect=a.*(b.*r0);
divu=[(-2*(3*uu.^2+vv.^2))./(uu.^2+vv.^2)];
%El tensor de tensiones que viene dado por la fórmula arriba dada se ha realizado a parte (analíticamente) debido a problemas con las matrices y sus operaciones
t=[-6.*(3.*uu.^2+vv.^2)./(uu.^2+vv.^2),-6.*uu.*vv;-2*uu.*vv,-2.*(3.*uu.^2+vv.^2)./(uu.^2+vv.^2)];
%una vez hallado el tensor de tensiones en un medio elástico lineal, isótropo y homogéneo
%realizamos una comparación de las gráficas del módulo de la divergencia y el módulo del rotacional
%con las tensiones normales en la dirección que marca g_u y g_v
%mod=[gu./(abs(gu))];
%ten_nor_gu=mod.*t.*mod;
tnoru=[(-6*vv.*(3*uu.^2+vv.^2))./((uu.^2+vv.^2).^2)-((8*uu.^2.*vv.^2)./(uu.^2+vv.^2))-(((2*uu.^2).*(3*uu.^2+vv.^2))./((uu.^2+vv.^2).^2))];
%lo mismo pero ahora en direccion de la g_v
%mod1=[gv/(abs(gv))];
%ten_nor_gv=mod1.*ten_tensiones.*mod1; hubiera sido una opción para realizarlo en matlab pero las matrices no son cuadradas y crean muchos problemas
tnorv=[((-6*uu.^2.*(3*uu.^2+vv.^2))./((uu.^2+vv.^2).^2))+((8*uu.^2.*vv.^2)./(uu.^2+vv.^2))-((2*vv.^2.*(3*uu.^2+vv.^2))./((uu.^2+vv.^2).^2))];
figure(1)
mesh(xx,yy,tnoru)
xlabel('Eje X')
ylabel('Eje Y')
colorbar
figure(2)
mesh(xx,yy,tnorv)
xlabel('Eje X')
ylabel('Eje Y')
colorbar
figure (3)
mesh(xx,yy,divu)
xlabel('Eje X')
ylabel('Eje Y')
colorbar
%la figura del rotacional no hace falta ya que el rotacional es cero y su comparativa no es necesaria.
INTERPRETACIÓN
Podemos observar que las gráficas de la [math]\nabla\cdot\overrightarrow{u}[/math] y de las tensiones normales tienen cierta similitud en cuanto a su forma pero sus ejes son diferentes, una transcurre en el ámbito positivo mientras que la otra esta en espectro negativo.
4.1 Tensión de Von Mises o criterio del fallo elástico
Una vez definida nuestra matriz de tensor de tensiones como [math] \sigma [/math], podremos calcular los autovalores de dicha matriz para hallar la Tensión de Von Mises siguiendo la ecuación:
[math] \sqrt{\frac{(\sigma_1-\sigma_2)^2-(\sigma_2-\sigma_3)^2-(\sigma_3-\sigma_1)^2}{2}} [/math]
La Tensión de Von Mises es una magnitud física proporcional a la energía de distorsión, es decir, una ecuación que actúa como un criterio de fallo elástico sobre un sólido deformable, cuando el tensor de tensiones rebasa un cierto valor y empieza a comportarse de forma plástica. Es por lo tanto una medida de seguridad para materiales dúctiles.
%tomaremos las incógnitas necesarias de apartados anteriores
s1=(-14*uu.^2-6*vv.^2)./(vv.^2+uu.^2);
s2=(-10*uu.^2-2*vv.^2)./(vv.^2+uu.^2);
s3=(-6*uu.^2-2*vv.^2)./(vv.^2+uu.^2);
sVM=sqrt((s1-s2).^2+(s2-s3).^2+(s3-s1).^2)/(sqrt(2))
figure(2)
mesh(xx,yy,sVM)
xlabel('Eje X');ylabel('Eje Y')
colorbar
maxv=(max(max(sVM)))
El valor máximo de la Tensión de Von Mises es 6.9282
5 Masa de la Placa
Para finalizar este estudio realizaremos un último añadido, calcularemos la masa de la placa. Disponemos de la densidad d(x,y) y una superficie anteriormente hallada. Por lo tanto la masa queda definida mediante una integral doble de este modo:
[math] \iint d(x,y)dxdy [/math]
%tomaremos los datos necesarios de apartados anteriores
h1=1/100;
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 -4.4462e-20. Sabemos que sale negativo pero es una transferencia de errores debido a que no se hace una integral como tal y la suma de zonas negativas en la matriz puede ocasionar datos negativos, como ha pasado.
6 Conclusiones Finales
Una vez realizado el estudio hemos podido observar en las diferentes gráficas realizadas que la placa elálastica depende en gran medida del tipo de campo vectorial' que se aplica sobre ella, siendo de gran ayuda para nosotros los ingenieros la Tensión de Von Mises. En la gráfica se aprecia que el criterio de fallo elástico se dispara siguiendo una forma curvilínea cuyo rango es de 6,5, que a nuestro parecer es un cambio drástico. Si miramos las Temperaturas en las unidades que estén su variación no supera la unidad. Lo que implica que con poca variación de temperatura hay una considerable deformación por lo tanto deducimos que es un sólido esta hecho de un material con unas propiedades plásticas elevadas.
