Visualización de campos escalares y vectoriales en elasticidad. Grupo 18
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Deformaciones de una placa plana en elasticidad. Grupo 18 |
| Asignatura | Teoría de Campos |
| Curso | 2014-15 |
| Autores | Pablo Roman Vegue Sanchez, Alberto Rodriguez Soto, Juan Manuel Rueda Olmedo, Santiago Gomez Fernandez, David Toledo Menendez, Jose Ramon Jimenez Villaseca |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
Contenido
1 Planteamiento del Problema
Se pretende conocer la visualizacion de un campo escalar y vectorial en elasticidad sobre la malla formada por la interseccion de dos parabolas P1= [math]18·y−81^2−1=0[/math] y P2= [math]2·y+x^2−1=0[/math] .Para ello empleamos un sistema de coordenadas apropiado (coordenadas parabolicas): [math]x=u·v[/math]
[math]y=\frac{(u^2−v^2)}{2}[/math] con u,v definidas en (u,v) ∈ [1/3,1] × [−1,1].
Con todo lo anterior implentando en matlab el mallado anterior es:
u=linspace(1/3,1,20); %creamos el vector u en el intervalo [1/3,1]
v=linspace(-1,1,20); %creamos el vector v en el intervalo[-1,1]
[uu,vv]=meshgrid(u,v); %matrices de las coordenadas u y v
xx=uu.*vv;
yy=1/2.*(uu.^2-vv.^2);
mesh(xx,yy,0*xx) %dibujamos el mallado
axis([-1 1 -1 1])%en los sucesivos pasos ajustamos y damos forma a los ejes
xlabel('eje x')
ylabel('eje y')
title('Mallado de los puntos interiores del solido')
grid off
2 Representacion de lineas coordenadas y vectores de la base natural
Para representar los vectores de la base natural en cada punto del mallado hay que tener en cuenta que esta varia en funcion del punto de la placa en que nos encontremos. Por otro lado las lineas coordenadas se obtienen fijando un parametro y haciendo variar el otro. Si calculamos los vectores de la base natural, estos son: [math] \vec{g_u}=v\hat{e_1} +2u \hat{e_2}[/math] : [math] \vec{g_v}=u\hat{e_1} -v^2 \hat{e_2}[/math] : [math]\vec{g}_w=\hat{e_3} [/math]
Para obtenerlo con matlab:
hold on
u=linspace(1/3,1,20);
v=linspace(-1,1,20);
[uu,vv]=meshgrid(u,v);
xx=uu.*vv;
yy=1/2.*(uu.^2-vv.^2);
axis([-1 1 -1 1])
xlabel('eje x')
ylabel('eje y')
title('Vectores de la base natural y lineas coordenadas')
grid off
gux=vv; %componente en la dirección i del vector gu de la base natural
guy=uu; %componente en la dirección j del vector gu de la base natural
quiver(xx,yy,gux,guy)%representamos los vectores gu en cada punto del mallado
gvx=uu; %componente en la dirección i del vector gv de la base natural
gvy=-vv; %componente en la dirección j del vector gv de la base natural
quiver(xx,yy,gvx,gvy)%representamos los vectores gv en cada punto del mallado
for k=1:20; %creamos yy1 correspondiente a las ordenadas de las primeras lineas coordenadas
for p=1:20;
yy1(k,p)=1/2*(((xx(k,p)^2)/(vv(k,p)^2))-(vv(k,p)^2));
end
end
plot(xx,yy1,'k') %representamos las primeras lineas coordenadas
for k=1:20; %creamos yy2 correspondiente a las ordenadas de las segundas lineas coordenadas.
for p=1:20;
yy2(k,p)=1/2*((uu(k,p)^2)-((xx(k,p)^2)/(uu(k,p)^2)));
end
end
plot(xx,yy2,'k')%representamos las segundas lineas coordenadas.
hold off
3 Curvas de nivel de la temperatura
Dado el campo escalar [math]T(x,y)=e^{-y} [/math] que representa la Temperatura, las curvas de nivel obtenidas serian:
u=linspace(1/3,1,20);
v=linspace(-1,1,20);
[uu,vv]=meshgrid(u,v);
xx=uu.*vv;
yy=1/2.*(uu.^2-vv.^2);
f=exp(-yy); %definimos la funcion temperatura
contour(xx,yy,f)%dibujamos las curvas de nivel de la temperatura
colorbar
axis([-1 1 -1 1])
xlabel('eje x')
ylabel('eje y')
title('Curvas de nivel de la temperatura')
Como podemos observar en el grafico adjunto, la leyenda representa con color rojo aquellos en que la temperatura es maxima, se corresponde con la parte inferior de la placa, mientras que el color azul, representa la minima temperatura.
El gradiente de la temperatura es:
hold on
u=linspace(1/3,1,20);
v=linspace(-1,1,20);
[uu,vv]=meshgrid(u,v);
xx=uu.*vv;
yy=1/2.*(uu.^2-vv.^2);
T=inline('exp(-yy)','xx','yy');%funcion temperatura
Tyy=inline('-exp(-yy)','xx','yy');%derivada respecto a y de la temperatura
zz=T(xx,yy);%valores de la funcion temperatura
V=Tyy(xx,yy); %componente en la direccion i del gradiente de la temperatura
U=zeros(size(V)); %componente en la direccion j del gradiente de la temperatura
contour(xx,yy,zz)%dibujamos las lineas de nivel de la temperatura
quiver(xx,yy,U,V)%representamos el campo gradiente de la temperatura
axis([-1 1 -1 1])
xlabel('eje x')
ylabel('eje y')
title('Gradiente y curvas de nivel de la temperatura')
hold off
El campo gradiente de la temperatura se visualiza en la imagen de forma ortogonal a las curvas de nivel.
4 Campo de desplazamientos w(u,v)
Al aplicar una fuerza sobre la placa, esta experimenta un cierto desplazamiento que viene determinado por el vector [math] \vec{w}(u,v) [/math] : Este vector será [math] {w}\vec (u,v)= \vec{a} (\vec{b} \vec{r_o)^{2}} [/math] siendo [math] \vec{a}= \frac{\vec{g}_u}{|\vec{g_u}|}[/math] y [math] \vec{b}=-4 \frac{\vec{g}_u}{|\vec{g_u}|}[/math].
[math]\vec{w}=-4(u^2+v^2)^{1/2}\vec{g_u}=-4v(u^2+v^2)^{1/2} \hat{e_1} + 4u(u^2+v^2)^{1/2} \hat{e_2}[/math]
El desplazamiento provocado por el campo [math]\vec{w}[/math] será:
u=linspace(1/3,1,20);
v=linspace(-1,1,20);
[uu,vv]=meshgrid(u,v);
xx=uu.*vv;
yy=1/2.*(uu.^2-vv.^2);
Wx=4.*vv.*(uu.^2+vv.^2).^(1/2);%componente i del campo u de desplazamientos
Wy=4.*uu.*(uu.^2+vv.^2).^(1/2);%componente j del campo u de desplazamientos
quiver(xx,yy,Wx,Wy)%dibujamos el campo de desplazamintos u
axis([-1 1 -1 1])
xlabel('eje x')
ylabel('eje y')
title('Campo de desplazamientos w')
Una vez aplicado el campo [math]\vec{w}[/math] la placa se deforma adoptando la siguiente forma:
Implementacion en matLAB:
u=linspace(1/3,1,20);
v=linspace(-1,1,20);
[uu,vv]=meshgrid(u,v);
xx=uu.*vv;
yy=1/2.*(uu.^2-vv.^2);
subplot(1,2,1)
mesh(xx,yy,0*xx)%dibujamos el solido sin desplazar
xlabel('eje x')
ylabel('eje y')
title('Solido sin desplazar')
axis([-1 1 -1 1])
grid off
Wx=4.*vv.*(uu.^2+vv.^2).^(1/2);%componente i del campo u de desplazamientos
Wy=4.*uu.*(uu.^2+vv.^2).^(1/2);%componente j del campo u de desplazamientos
xxt=xx+Wx;%solido tras el desplazamiento
yyt=yy+Wy;%solido tras el desplazamiento
subplot(1,2,2)
mesh(xxt,yyt,0*xx)%dibujamos el solido tras el desplazamiento
xlabel('eje x')
ylabel('eje y')
title('Solido tras aplicarle el desplazamiento')
grid off
view(2)
5 Divergencia
La divergencia del campo vectorial nos proporciona el flujo que sale por cada punto por unidad de volumen
[math]\nabla· \vec{w}= \frac{1}{ \sqrt{g} } \frac{\partial \sqrt{g} u^{i} }{\partial u^i}= \frac{1}{u^2+v^2} \frac{\partial}{\partial u}((u^2+v^2)(4(u^2+v^2)^{1/2}) = \frac{12}{(u^2+v^2){1/2} } [/math]
Implementacion en matLAB:
u=linspace(1/3,1,20);
v=linspace(-1,1,20);
[uu,vv]=meshgrid(u,v);
xx=uu.*vv;
yy=1/2.*(uu.^2-vv.^2);
div=12./(uu.^2+vv.^2).^(1/2);
mesh(xx,yy,div)
surf(xx,yy,div)
axis([-1,1,-1,1])
xlabel('eje x')
ylabel('eje y')
title('Divergencia de w')
colorbar
Conforme a la definicion de divergencia, en el grafico podemos observar que aquellos puntos de color rojo son los que tienen mayor divergencia, y como es positiva, actuan como fuentes. En el caso de los puntos azules, tienen menor divergencia. R
6 Rotacional
El rotacional mide la cantidad de giroal rededor del vector normal a la superficie.
[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} \\ 4/(u^2+v^2)^{1/2} &0&0\end{bmatrix}= 4v/(u^2+v^2)^{2}[/math]
Implementacion en matLAB:
u=linspace(1/3,1,20);
v=linspace(-1,1,20);
[uu,vv]=meshgrid(u,v);
xx=uu.*vv;
yy=1/2.*(uu.^2-vv.^2);
rot=4.*vv./(uu.^2+vv.^2).^(3/2);%modulo del rotacional de u
subplot(1,2,1)
mesh(xx,yy,rot)
surf(xx,yy,rot)
axis([-1 1 -1 1])
xlabel('eje x')
ylabel('eje y')
title('Rotacional de w')
view(2)
subplot(1,2,2)
mesh(xx,yy,rot)
surf(xx,yy,rot)
axis([-1 1 -1 1])
xlabel('eje x')
ylabel('eje y')
title('Rotacional de w')
colorbar
7 Tensiones normales en la direcciones gu y gv
Dado є(u) en la placa podemos describir los desplazamientos provocados por las tensiones a traves de la formula:
[math]σ=λ\nabla·\vec{u}1+2μЄ[/math] siendo Є(u): [math]\epsilon (\vec{u})=\frac{ \nabla \vec{u}+ \nabla \vec{u}^t}{2}[/math]
siendo [math]λ=μ=1[/math]
Implementando en matLAB
u=linspace(1/3,1,20);
v=linspace(-1,1,20);
[uu,vv]=meshgrid(u,v);
xx=uu.*vv;
yy=1/2.*(uu.^2-vv.^2);
%Tensiones en la direccion de gu
Ngu=(12*vv.^2+16.*uu.*vv.^2+12.*uu.^2)./(uu.^2+vv.^2).^(3/2);
subplot(1,2,1)
surf(xx,yy,Ngu)
colorbar
axis([-1,1,-1,1])
xlabel('eje x')
ylabel('eje y')
title('Tensiones normales en la direccion gu')
%Tensiones en la direccion de gv
Ngv=(12.*uu.^2-8.*uu.^3-8.*vv.^2.*uu+12.*vv.^2)./(uu.^2+vv.^2).^(3/2);
subplot(1,2,2)
surf(xx,yy,Ngv)
colorbar
axis([-1,1,-1,1])
xlabel('eje x')
ylabel('eje y')
title('Tensiones normales en la direccion gv')
Observando las graficas correspondientes a la divergencia, rotacional y tensiones normales, llegamos a la conclusion de que estas ultimas guardan cierta similitud con la divergencia
8 Tensiones de Von Mises
A partir de la tensión de Von Mises que se define con la siguiente fórmula: [math] \sigma _{VM}= \sqrt{ \frac{( \sigma _1- \sigma _2)^2+(\sigma_2-\sigma_3)^2+(\sigma_3-\sigma_1)^2 }{2} }[/math] donde [math] \sigma_1 \sigma_2 [/math] y [math] \sigma_3 [/math] son los autovalores del tensor σ (tensiones principales), podemos conocer el comportamiento plastico de la placa y asi saber que puntos pueden llegar a colapsar
Implementacion en matLAB:
hold on
u=linspace(1/3,1,20);
v=linspace(-1,1,20);
[uu,vv]=meshgrid(u,v);
xx=uu.*vv;
yy=1/2.*(uu.^2-vv.^2);
x=u.*v;
y=1/2.*(u.^2-v.^2);
for i=1:20
for j=1:20
A=[(12+8.*u(i))./(u(i).^2+v(j).^2).^(1/2),4.*v(j)./(u(i).^2+v(j).^2).^(1/2),0;4.*v(j)./(u(i).^2+v(j).^2).^(1/2),(12)./(u(i).^2+v(j).^2).^(1/2),0;0,0,(12)./(u(i).^2+v(j).^2).^(1/2)];
a=eig(A);
vonmises(i,j)=sqrt(((a(1)-a(2)).^2+(a(2)-a(3)).^2+(a(3)-a(1)).^2)./2);
end
end
subplot(1,2,1)
surf(xx,yy,vonmises)
colorbar
axis([-1,1,-1,1])
xlabel('eje x')
ylabel('eje y')
title('Tensiones de vonmises')
subplot(1,2,2)
surf(xx,yy,vonmises)
colorbar
axis([-1,1,-1,1])
xlabel('eje x')
ylabel('eje y')
title('Tensiones de vonmises')
hold off
A partir de la barra de colores observamos que los puntos de color rojo estan en la cresta de la figura y experimentan una mayor tension de Von Mises.
9 Masa de la placa
A partir de la funcion de densidad calculamos la masa de la placa tomando un elemento notablemente pequeño del mallado, posteriormente se obtiene el area de ese elemento y se realiza para todo el dominio de la placa, efectuando una suma final de todos los valores obtenidos. Esta resulta ser 4.2809e-006.
Implementacion en matLAB:
u=linspace(1/3,1,20);
v=linspace(-1,1,20);
[uu,vv]=meshgrid(u,v);
xx=uu.*vv;
yy=1/2.*(uu.^2-vv.^2);
h=1/100;
f=abs(xx).*exp(1).^(-1./(yy.^2));
a=h^2.*f;
masa=sum(sum(a))
