Visualización de campos escalares y vectoriales en elasticidad. Grupo 17. Trabajo 5
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Deformaciones de una placa plana. Grupo 17-C |
| Asignatura | Teoría de Campos |
| Curso | 2014-15 |
| Autores | Claudia Jalón Manzano, Nacho Mollá Carcaño, Jorge Javier Rodriguez Anzules, Pedro Torrecilla, Pablo Revuelta Aragón |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
Contenido
- 1 Introducción
- 2 Líneas coordenadas y vectores de la base natural
- 3 Influencia de un foco de calor
- 4 Gradiente de T y curvas de nivel
- 5 Campo de desplazamientos
- 6 Visualización gráfica del desplazamiento producido en el sólido
- 7 Divergencia
- 8 Rotacional
- 9 Tensor de tensiones
- 10 Tensiones de Von Mises
- 11 Masa total del sólido
1 Introducción
Consideramos una placa plana (en dimensión 2) que ocupa la región comprendida entre las parábolas P1= [math]18y−81x^2−1=0[/math] y P2= [math]2y+x^2−1=0[/math]. Para representarla usaremos un sistema de coordenadas adaptado a la geometría que nos da:
[math]x=uv[/math]
[math]y=\frac{(u^2−v^2)}{2}[/math] con u,v definidas en (u,v) ∈ [1/3,1] × [−1,1].
La representación de los puntos interiores de la placa sólida nos queda como podemos observar:
%mallado del sólido
h=1/20; %paso de muestreo
u=1/3:h:1; %intervalo[1/3,1]
v=-1:h:1; %intervalo [-1,1]
[U,V]=meshgrid(u,v); %matrices de las coordenadas u y v
figure(1)
xx=U.*V; %paramerización
yy=1/2*(U.^2-V.^2);
mesh(xx,yy,0*xx) %dibujar la malla
axis([-1,1,-1,1]) %selecciona la región
2 Líneas coordenadas y vectores de la base natural
Las líneas coordenadas y los vectores de la base natural irán cambiado en direccion según el punto de la placa, ya que la base natural en estas coordenadas no es constante. La base natural será la siguiente: [math] \vec{g_u}=v\hat{e_1} +u \hat{e_2}[/math] : [math] \vec{g_v}=u\hat{e_1} -v \hat{e_2}[/math] : [math]\vec{g}_w=\hat{e_3} [/math]
La representación la podemos observar en este gráfico:
%lineas coordenadas y vectores de la base natural
figure(2)
guu=V; %componente i de gu
guv=U; %componente j de gu
gvu=U; %componente i de gv
gvv=-V; %componente j de gv
hold on
mesh(xx,yy,0*xx); %dibujar la malla
quiver(xx,yy,guu,guv) %dibujar el campo vectorial de los vectores gu
quiver(xx,yy,gvu,gvv) %dibujar el campo vectorial de los vectores gv
hold off
3 Influencia de un foco de calor
La temperatura proviene de un foco de calor dada por el campo escalar [math]T(x,y)=(8-y^2+2y)e^{-(x)^2} [/math]
%temperatura del sólido T
figure(3)
T=(8-yy.^2+2.*yy).*exp(-(xx).^2); %función de la temperatura (campo esalar)
surf(xx,yy,T) %visualización campo escalar
colorbar %escala de colores
axis([-1,1,-1,1])
La temperatura (T) es una función continua que depende de x e y, que a su vez dependen de u y v . Si estudiamos los extremos de la función temperatura, observamos que obtenemos un máximo absoluto en (x,y)=(0,1). Sin embargo, este punto no pertenece al dominio, dado que nuestra función está acotada. Por ello, debemos buscar el máximo relativo en el contorno, que al ser una función suave se encontrará en el punto de la frontera más próximo al máximo absoluto de la función T, es decir, la temperatura es máxima en (x,y)=(0,0.5).
4 Gradiente de T y curvas de nivel
El gradiente de una función escalar expresa la dirección en la cual el campo crece más rápido. Por otra parte, las curvas de nivel expresan los puntos que se encuentran a la misma temperatura. Por tanto, si tenemos dos curvas de nivel y estamos en un punto de la menor, el gradiente será el vector de mínima distancia a la otra curva de nivel, siendo por tanto perpendicular a ambas.
[math]\nabla T(x,y)= \frac{\partial T}{\partial x} \hat{e_1}+\frac{\partial T}{\partial y}\hat{e_2}=[-2xe^{-x^2}(8-y^2+2y)]\hat{e_1}+[(-2y+2)e^{-x^2}]\hat{e_2}[/math]
%gradiente de T y curvas de nivel
figure(4)
Tx=-2.*xx.*exp(-(xx).^2).*(8-yy.^2+2.*yy); %derivada de la temperatura respecto a x
Ty=(-2.*yy+2).* exp(-(xx).^2); %derivada de la temperatura respecto a y
subplot(1,3,1),quiver(xx,yy,Tx,Ty) %representación del gradiente
axis([-1,1,-1,1]);
subplot(1,3,2),contour(xx,yy,T,10) %respresentación de las curvas de nivel
axis([-1,1,-1,1]);
subplot(1,3,3)
hold on
quiver(xx,yy,Tx,Ty) %representación del gradiente y de las curvas de nivel
contour(xx,yy,T,10)
axis([-1,1,-1,1]);
hold off
5 Campo de desplazamientos
Una fuerza determinada aplicada sobre nuestro sólido ha provocado un desplazamiento del mismo que viene dado por [math] \vec{u}(x,y) [/math] . Este vector será [math] \vec {u}(u,v)= \vec{a} (\vec{b} \vec{r_o}) [/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]. Como hemos hallado anteriormente [math]\vec{g_u}=v \hat{e_1} +u \hat{e_2} [/math]. Tomaremos: [math]\vec{r_o}= x_1 \hat{i} + x_2 \hat{j} + x_3 \hat{k}= uv{e_1}+ \frac{1}{2}(u^2-v^2){e_2}+w{e_3}[/math] Con todo esto: [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] La representación del campo de desplazamiento [math]\vec{u}[/math] será la siguiente:
%campo vectorial de desplazamientos u
figure(5)
Ux=-2.*U.*V; %componente i del campo de desplazamientos
Uy=-2.*U.^2; %componente j del campo de desplazamientos
quiver(xx,yy,Ux,Uy) %representacion del campo vectorial de desplazamientos
axis([-1,1,-1,1])
6 Visualización gráfica del desplazamiento producido en el sólido
Lo que se trata de mostrar en este apartado es el desplazamiento producido en el sólido provocado por la fuerza de la que hemos hablado anteriormente. Para ello debemos considerar [math] \vec{r_o} [/math] como el vector posición inicial de los puntos del sólido y [math] \vec{d}=\vec{r_o}+\vec{u}[/math] como el vector de los puntos del sólido al final del desplazamiento. La visualización de este desplazamiento será la siguiente:
%desplazamiento del sólido
figure(6)
xd=xx+Ux; %componente i del sólido desplazado
yd=yy+Uy; %componente j del sólido desplazado
subplot(1,2,1), mesh(xx,yy,0*xx) %dibujo del sólido antes de desplazarse
axis([-1,1,-2,1]) %región del gráfico
subplot(1,2,2), mesh(xd,yd,0*xd) %dibujo del sólido despues de desplazarse
axis([-1,1,-2,1]) %región del gráfico
7 Divergencia
La divergencia nos sirve para medir la diferencia entre el flujo saliente y el flujo entrante de un campo vectorial sobre la superficie que rodea a un volumen de control, por tanto, si el campo tiene "fuentes" la divergencia será positiva, y si tiene "sumideros", la divergencia será negativa.
[math]\nabla· \vec{u}= \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)(-2u)) =\frac{1}{u^2+v^2} \frac{\partial}{\partial u}(-2u^3-2uv^2)= \frac{1}{u^2+v^2}(-6u^2-2v^2)= \frac{-6u^2-2v^2}{u^2+v^2} [/math]
Por lo que es trivial que [math]\nabla ·\vec{u} \leq 0[/math] , por lo que será un sumidero.
%divergencia
figure(7)
div=(-6*U.^2-2*V.^2)./(U.^2+V.^2); %campo escalar Divergencia
surf(xx,yy,div) %dibujar divergencia
axis([-1,1,-1,1]) %selecciona la región
8 Rotacional
A partir del operador rotacional se puede establecer la tendencia del campo vectorial a inducir rotación en un punto alrededor de un vector. El rotacional obtenido en nuestro caso es: [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{4uv}{u^2+v^2} g_{w} [/math]
%rotacional
figure(8)
rot=abs((4.*U.*V)./(U.^2+V.^2); %modulo del rotacional
surf(xx,yy,rot) %dibujar rotacional
axis([-1,1,-1,1]) %selecciona la región
9 Tensor de tensiones
Definiendo la parte simétrica del tensor gradiente de u, є(u). En un medio elástico lineal, isótropo y homogéneo los desplazamientos permiten escribir el tensor de tensiones a través de la siguiente fórmula :
[math]σ=λ\nabla·\vec{u}1+2μЄ[/math] siendo Є(u): [math]\epsilon (\vec{u})=\frac{ \nabla \vec{u}+ \nabla \vec{u}^t}{2}[/math] con [math]λ=μ=1[/math].
Procederemos a calcular el gradiente de [math]\vec{u}[/math] con la fórmula de las derivadas parciales covariantes: [math]u^{i}, _j = \frac{\partial u^i}{\partial x^j}+\Gamma^i_{jk}u^k[/math]
[math]u^{1}, _1 = \frac{\partial u^u}{\partial u} + \Gamma ^1_{11} u^u=-2- \frac{2u^2}{u^2+v^2}[/math] [math]u^{1}, _2 = \frac{\partial u^u}{\partial v} + \Gamma ^1_{12} u^u= \frac{-2uv}{u^2+v^2} [/math] [math]u^{1}, _{3} = \frac{\partial u^u}{\partial w} + \Gamma ^1_{13} u^u=0 [/math]
[math]u^2,_1= \frac{\partial u^v }{\partial u} + \Gamma ^{2}_{11}u^u= \frac{2uv}{u^2+v^2} [/math] [math]u^2,_2= \frac{\partial u^v }{\partial v} + \Gamma ^{2}_{12}u^u= \frac{-2u^2}{u^2+v^2} [/math] [math]u^2,_3= \frac{\partial u^v }{\partial w} + \Gamma ^{2}_{13}u^u=0 [/math]
[math]u^3,_1= \frac{\partial u^w }{\partial u} + \Gamma ^{3}_{11}u^u=0 [/math] [math]u^3,_2= \frac{\partial u^w }{\partial v} + \Gamma ^{3}_{12}u^u=0 [/math] [math]u^3,_3= \frac{\partial u^w }{\partial w} + \Gamma ^{3}_{13}u^u=0 [/math]
Una vez calculadas las componentes del gradiente, sacaremos su parte simétrica Є(u) haciendo la semisuma del gradiente y su transpuesta.
[math]\nabla\vec{u}=\begin{bmatrix}-2- \frac{2u^2}{u^2+v^2} & \frac{-2uv}{u^2+v^2} & 0 \\\frac{2uv}{u^2+v^2} & \frac{-2u^2}{u^2+v^2} & 0\\0&0&0 \end{bmatrix}[/math], [math] \nabla \vec{u}^t=\begin{bmatrix}-2- \frac{2u^2}{u^2+v^2} & \frac{2uv}{u^2+v^2} & 0 \\\frac{-2uv}{u^2+v^2} & \frac{-2u^2}{u^2+v^2} & 0\\0&0&0 \end{bmatrix} [/math], [math] Є(\vec{u})=\begin{bmatrix}\frac{-8u^2-4v^2}{2(u^2+v^2)} & 0 & 0 \\\ 0 & \frac{-4u^2}{2(u^2+v^2)} & 0\\0&0&0 \end{bmatrix} [/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} \frac{-8u^2-4v^2}{2(u^2+v^2)}&0&0 \\0&\frac{-4u^2}{2(u^2+v^2)}&0\\0&0&0 \end{bmatrix}=\begin{bmatrix}\frac{-14u^2-6v^2}{u^2+v^2} & 0 & 0 \\\ 0 & \frac{-10u^2-2v^2}{u^2+v^2} & 0\\0 &0& \frac{-6u^2-2v^2}{u^2+v^2} \end{bmatrix} [/math]
Ahora ya tenemos el tensor de tensiones, a continuación con matlab dibujaremos las tensiones normales en las direciones [math] \vec{g}_u [/math] y [math] \vec{g}_v [/math].
%tensor de tensiones
figure(8)
Eu=(-14*U.^2-6*V.^2)./(U.^2+V.^2); %Campo escalar
Ev=(-10*U.^2-2*V.^2)./(U.^2+V.^2); %Campo escalar
subplot(1,2,1);surf(xx,yy,Eu) %mallado
axis([-1,1,-1,1]) %región
view(2) %punto de vista
subplot(1,2,2);surf(xx,yy,Ev) %mallado2
axis([-1,1,-1,1]) %región
El sentido físico de la divergencia es cuánto se están estirando o contrayendo los puntos por acción de las fuerzas y el rotacional es la rotación inducida en los puntos por lo mismo, sumando ambos nos daría la deformación de la placa, lo cual está muy relacionado con el tensor de tensiones, de ahí la similitud que hay.
10 Tensiones de Von Mises
La tensión de Von Mises es una magnitud escalar que indica el comportamiento plástico del sólido se calcula con la siguiente fórmula, donde [math] \sigma_1 \sigma_2 [/math] y [math] \sigma_3 [/math] son los autovalores del tensor σ (tensiones principales): [math] \sigma _{VM}= \sqrt{ \frac{( \sigma _1- \sigma _2)^2+(\sigma_2-\sigma_3)^2+(\sigma_3-\sigma_1)^2 }{2} }[/math]
%tensiones de Von Misses
figure(9)
e1=(-14*U.^2-6*V.^2)./(V.^2+U.^2); %autovalor 1
e2=(-10*U.^2-2*V.^2)./(V.^2+U.^2); %autovalor 2
e3=(-6*U.^2-2*V.^2)./(V.^2+U.^2); %autovalor 3
RVM=sqrt(((e1-e2).^2+(e2-e3).^2+(e3-e1).^2)./2); %fórmula de Von Misses (campo escalar)
surf(xx,yy,RVM) %reseprentación de la fórmula en el mallado
axis([-1,1,-1,1]) %selecciona la región
Como podemos observar la tensión de Von Mises es mayor en el eje de nuestro sólido esto indica que es la parte más propensa a sufrir una deformación plástica, esto puede ser por la naturaleza de las fuerzas o porque es la seción más pequeña.
11 Masa total del sólido
La densidad de la placa viene dada por la función:
Calculamos la masa con el siguiente programa. Utilizamos la defición de integral para calcularla, es decir, la suma de los diferenciales de masa dividiendo cada eje entre 100. Si ampliaramos esta división cada vez se acercaría más al valor real de la masa, siendo ya suficientemente exacta con este grado de división.
El resultado del programa es:
%masa del sólido
hh=1/100; %paso de muestreo
u=1/3:hh:1; %intervalo [1/3,1]
v=-1:hh:0; %intervalo [-1,1]
[U,V]=meshgrid(u,v); %matrices de las coordenadas u y v
X=U.*V; %paramerización
Y=1/2*(U.^2-V.^2);
D=X.*Y.*exp(-1./X.^2); %función de densidad
a=(hh.^2).*D;
M1=abs(sum(sum(a))) %masa
M1 = 6.5985e^-5
hh=1/100; %paso de muestreo
u=1/3:hh:1; %intervalo [1/3,1]
v=-1:hh:0; %intervalo [-1,1]
[U,V]=meshgrid(u,v); %matrices de las coordenadas u y v
X=U.*V; %paramerización
Y=1/2*(U.^2-V.^2);
D=X.*Y.*exp(-1./X.^2); %función de densidad
a=(hh.^2).*D;
M2=abs(sum(sum(a))) %masa
M2 = 6.5985e^-5
M1+M2=0.088921
