Visualizacion de campos escalares y vectoriales en elasticidad (trabajo 7)
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Deformaciones de una placa plana.C |
| Asignatura | Teoría de Campos |
| Curso | 2015-16 |
| Autores | Pablo Goenechea |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
0 INTRODUCCIÓN
El objetivo de este trabajo es estudiar la aplicacion de campos escalares y vectoriales sobre una placa en dos dimensiones.
Un Campo escalar es la reresentación espacial de una magnitud escalar, es decir, un número y en nuestro caso la termperatura. Un Campo Vectorial es la representación espacial de una magnitud vectorial, asociando a cada punto del espacio un vector, en este caso, el desplazamiento producido por una fuerza.
Tambien estudiaremos operadores como la divergencia y el rotacional, las tensiones producidas y la masa total dada la densidad de la placa.
1 MALLADO DEL SÓLIDO
La placa a estudiar ocupa la región comprendida entre dos parábolas: P1 : 18y−81x2−1 = 0, y P2 : 2y +x2 −1 = 0
Para representarlas usaremos un sistema de coordenadas adaptado a la geometria que nos dan: x=uv; y=(u^2−v^2)/2; Hemos representado la placa en (1/3,-1)<(u,v)<(1,1), tomado un paso de muestreo de h=1/20 para u y v.
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 LÍNEAS COORDENADAS Y VECTORES DE LA BASE NATURAL PARA CADA PUNTO
Como la base natural no es constante para todos los puntos; cambiará, junto con las líneas coordenadas, de dirección según cada punto de la placa. Por ello definimos un sistema de coordenadas curvilíneas {0,(u,v,w)} ortogonal y no necesariamente unitario, cuyas componentes (u,v,w) son las derivadas parciales del vector posicion r(x,y,z) respecto de (u,v,w).
for i=1:40 a=u.*i; b=((u.^2)-(i.^2))./2; c=v.*i; d=((i.^2)-(v.^2))./2; end hold on mesh(x,y,0*x) quiver(x,y,v,u); quiver(x,y,u,-v); axis([-1.25,1.25,-1.25,1.25]) hold off view(2)
3 y 4 TEMPERATURA DEL SÓLIDO Y GRADIENTE
La Temperatura del sólido viene dada por la fórmula: T(x,y)=(x-0.5)^2*exp(-y^2). Nos piden dibujar las curvas de nivel y hallar el punto dónde la temperatura es MÁXIMA. Para ello, hallamos los valores de T para cada punto de la malla y creamos una hipermatriz cuya componentes serán (x,y,T). Con el comando contour podemos elegir el número de lineas de nivel que queremos que aparezcan (en nuestro caso 100).
El gradiente nos indica la direccion de máxima pendiente, por lo que teoricamente es ortogonal a las líneas de nivel. Para hallarlo calculamos las derivadas parciales Tx y Ty, que serán las componentes vectoriares del gradiente.
T=(xx.^2-xx+0.25).*exp(-(yy).^2); Tx=(2*xx.^2-xx).*exp(-(yy).^2); Ty=(-2*yy).*(xx.^2-xx+0.25).*exp(-(yy).^2); figure(1) contour(xx,yy,T,100) figure(2) quiver(xx,yy,Tx,Ty) figure(3) quiver(xx,yy,Tx,Ty) hold on contour(xx,yy,T,10,'r') hold off
Cómo vemos, el punto de Máxima temperatura se sitúa a la izquierda de la gráfica, con valores mínimos de "y". También podemos observar cómo el gradiente es ortogonal a las líneas de nivel.
5 y 6 CAMPO DE DESPLAZAMIENTOS (U)
La placa se ve sometida a una fuerza que produce unos desplazamientos U(u,v)=a(b*r0). Debemos calcular analíticamente los desplazamientos U(u,v) en la base (i,j,k). Esto es el desplazamiento de cada punto de la malla, que aplicaremos a la malla para hallar el desplazamiento U*r0=r0+u(r0) y así obtener la malla deformada.
clear all
h=1/20;
uu=1/3:h:1;
vv=-1:h:1;
[u,v]=meshgrid(uu,vv);
x=u.*v ;
y=(1/2).*((u.^2)-(v.^2));
wx=0.25.*u.*v.^2.*sqrt(u.^2+v.^2);
wy=0.25.*u.^2.*v.*sqrt(u.^2+v.^2);
subplot(1,3,2)
quiver(x,y,wx,wy)
title('Campo de deformaciones')
axis([-1.2,1.2,-1,1])
fx=x+wx;
fy=y+wy;
subplot(1,3,3)
surf(fx,fy,0*fx)
title('Placa deformada')
view(2)
subplot(1,3,1)
surf(x,y,0*fx)
title('Placa original')
view(2)
7 DIVERGENCIA
La divergencia es un operador que sirve para calcular los cambios de volúmen del sólido. El resultado es un escalar: si es positivo aumenta el volúmen y si es negativa se reduce. Apoyandonos en los calculos anteriores, calculamos el valor analíticamente y lo sustituimos para cada punto de la malla.
h=1/20; uu=1/3:h:1; vv=-1:h:1; [u,v]=meshgrid(uu,vv); x=u.*v ; y=(1/2).*((u.^2)-(v.^2)); div=0.25.*(u.^2+v.^2)./sqrt(u.^2+v.^2); subplot(1,2,1) surf(x,y,div);colorbar; view(2) subplot(1,2,2) surf(x,y,div);colorbar;
8 ROTACIONAL
El rotacional es un operador que sirve para medir cuanto tiende a rotar un punto inducido a un campo vectorial. Al trabajar en otro sistema de coordenadas al cartesiano, generalizamos la expresion para incluir que los valores de la base dependen de la posición, para ello usamos la notación de Eistein.
rot=(u^3*v-u*v^3)/(4*sqrt(u^2+v^2))
clear all h=1/20; uu=1/3:h:1; vv=-1:h:1; [u,v]=meshgrid(uu,vv); x=u.*v ; y=(1/2).*((u.^2)-(v.^2)); rot=(0.25.*(u.^3.*v)./((u.^2)+(v.^2)).^0.5); subplot(1,2,1) surf(x,y,rot);colorbar; view(2)
9 TENSIONES NORMALES
Se parte de las hipótesis de medio elástico lineal, isótropo y homogéneo. Los valores \(\lambda\) y \(\mu\) se conocen como los Coeficientes de Lamé, dependen de las propiedades elásticas de cada material.
Hemos vuelto a calcular \( \overrightarrow{u} \), pero esta vez dejando que Matlab realizara los cálculos. Las matrices \( \nabla\overrightarrow{u}\) y \(\nabla\overrightarrow{u}^t\) y la matriz del tensor de tensiones calculadas analíticamente son: \( \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}\)
\( \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} \)
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) 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(1) mesh(xx,yy,Gu) view(2) figure(2) mesh(xx,yy,Gv) view(2)
Comparando las gráficas de las tensiones vemos cómo la tensión normal en la direccion de gu se parece tanto al módulo de la divergencia como al rotacional.
10. TENSION DE VON MISES
La tensión de Von Miss es una magnitud escalar que se emplea como indicador de cuando un material inicia un comportamiento plástico (y no elástico puro). La tensión de Von Mises se define por la siguiente fórmula:
En la cual σ1, σ2 y σ3 son autovalores de σ.
sigma1=(-14*uu.^2-6*vv.^2)./(vv.^2+uu.^2); sigma2=(-10*uu.^2-2*vv.^2)./(vv.^2+uu.^2); sigma3=(-6*uu.^2-2*vv.^2)./(vv.^2+uu.^2); Mises=sqrt(((sigma1-sigma2).^2+(sigma2-sigma3).^2+(sigma3-sigma1).^2)./2); surf(xx,yy,Mises); axis([-1,1,-1,1]) max(max(Mises))
Observando la gráfica obtenida vemos cómo la tensión de von mises es opuesta a la tension normal en la direccion de v. Su valor máximo es 6.9282.
11. MASA TOTAL DE LA PLACA
El último apartado nos pide calcular la masa total de la placa a partir de su densidad: d(x,y)=exp(-y^2)*(abs(x)+1). Vamos a usar el método del trapecio para ello. Este método se basa en obtener la densidad en cada punto del mallado, para despues multiplicarla por un vector fila y otro columna, y asi obtener el escalar de la masa. Los valores negativos son posibles e ilógicos por lo que sólo tendremos en cuenta el valor absoluto.
N1=200; N2=200; a=1/3; b=1; c=-1; d=1; h1=(b-a)/N1; h2=(d-c)/N2; u=a:h1:b; v=c:h2:d; [uu,vv]=meshgrid(u,v); xx=uu.*vv; yy=(1/2).*((uu.^2)-(vv.^2)); d=exp(-(yy).^2)*(abs(xx)+1); D=abs(d);
w1=ones(N1+1,1); w1(1)=1/2; w1(N1+1)=1/2; w2=ones(N2+1,1); w2(1)=1/2; w2(N1+1)=1/2; result=h1*h2*w2'*D*w1
El resultado es 341.8297
