Visualización de campos escalares y vectoriales en elasticidad. (Grupo 15-C)
Contenido
1 Introducción
Vamos a realizar el estudio de una placa plana con forma de corona circular centrada en el origen y de radio interior 1 y radio exterior 2, al verse sometida a diversos campos. Para empezar, representamos el mallado sobre el que vamos a trabajar. La visualización de la misma nos sirve para situar la placa en nuestro espacio de trabajo.
2 Espacio de trabajo
En nuestro caso vamos a suponer que tenemos definidas dos cantidades físicas: la temperatura [math]T(\rho,\theta,t)[/math], que depende de las dos coordenadas polares [math](\rho,\theta)[/math] y el tiempo [math]t[/math], y los desplazamientos [math]\vec u(\rho,\theta,t)[/math]. De esta forma, si definimos [math]r_0(\rho,\theta)[/math] el vector de posición de los puntos de la placa en reposo, la posición de cada punto [math](\rho,\theta)[/math] de la placa en un instante de tiempo [math]t[/math] viene dada por:: [math] \vec r (\rho,\theta,t)= \vec r_{0}(\rho,\theta)+\vec u(\rho,\theta,t). [/math] Queda reflejado en Matlab con el siguiente código:
h=0.1 % Intervalos de separación
u=1:h:2; % Intervalo de rho [1,2]
v=0:h:2*pi+h; % Intervalo de theta [0,2*pi]
[uu,vv]=meshgrid(u,v); % Matrices de u y v
figure(1)
xx=uu.*cos(vv); % Parametrización
yy=uu.*sin(vv);
mesh(xx,yy,0*xx) % Dibujo del mallado
axis([-3,3,-3,3]) % Región del dibujo
view(2)
Tras la visualización deseamos estudiar cómo nuestra placa se va a comportar frente a un foco calorífico.
3 Comportamientos ante un foco calorífico
Nuestro foco calorífico se sitúa en el origen de coordenadas, siguiendo la expresión [math]T(\rho,\theta)=-\log(\rho+0.1)[/math]. Representado en Matlab con el siguiente código:
h=(0.1); % Intervalo de separación
u=[1:h:2]; % Intervalo de rho [1,2]
v=[0:h:2*pi+h]; % Intervalo de theta [0,2π]
[uu,vv]=meshgrid(u,v); % Matrices de rho y theta
figure(1)
xx=uu.*cos(vv); % Parametrización
yy=uu.*sin(vv);
f=-log(0.1+sqrt(xx.^2+yy.^2)); % Campo escalar
surf(xx,yy,f) % Dibujar la superficie
axis([-2,2,-2,2]) % Región del dibujo
colorbar
view(2)
Como observamos en la ilustración la temperatura desciende a medida que nos alejamos del origen de coordenadas. Por lo tanto la placa tiene mayor temperatura cuánto menor es su radio (valor de [math]\rho[/math]) , es decir cuanto más cerca del foco nos encontramos. En este caso el foco calorífico no varía en función de [math]\theta[/math].
3.1 Variación de la temperatura
Derivando la función [math]T(\rho,\theta)=-\log(\rho+0.1)[/math] en función de [math]\rho[/math] y [math]\theta[/math], obtenemos una expresión del gradiente de la temperatura que en coordenadas cartesianas nos da como resultado el siguiente campo vectorial::
[math]\nabla T = (\frac{-x}{x^2+y^2+0.1\sqrt{x^2+y^2}},\frac{-y}{x^2+y^2+0.1\sqrt{x^2+y^2}})[/math]
h=0.1; % Intervalo de separación
u=1:h:2; % Intervalo de rho [1,2]
v=0:h:2*pi+h; % Intervalo de theta [0,2*pi]
[uu,vv]=meshgrid(u,v); % Matrices de rho y theta
figure(1)
xx=uu.*cos(vv); % Parametrización
yy=uu.*sin(vv);
f=-log(0.1+sqrt(xx.^2+yy.^2)); % Campo escalar
contour(xx,yy,f) % Dibujar las líneas de nivel
hold on
fx=-xx./((xx.^2+yy.^2)+0.1.*sqrt(xx.^2+yy.^2)); % Derivada parcial respecto de X
fy=-yy./((xx.^2+yy.^2)+0.1.*sqrt(xx.^2+yy.^2)); % Derivada Parcial respecto de Y
quiver(xx,yy,fx,fy) % Dibujar el Campo Vectorial
axis([-2,2,-2,2]) % Región del dibujo
view(2)
colorbar
Si observamos la imagen obtenida con Matlab, podemos ver que las curvas de nivel de la temperatura parten del foco calorífico formando circunferencias concéntricas sin respetar una equidistancia determinada, ya que la función logarítmica no es lineal. Por tanto, estas líneas tienen mayor temperatura cuanto más próximas a él se localicen.
Por otro lado podemos observar que las flechas que representan el gradiente son radiales y que apuntan hacia el origen de coordenadas, debido a que el signo obtenido en las derivadas parciales es negativo. Podemos por tanto decir que el gradiente es menor cuanto más se aproxima al origen de coordenadas.
Gráficamente observamos que [math]\nabla T[/math] es ortogonal a las curvas de nivel de la temperatura.
4 Campo de desplazamiento
Vamos a suponer que sobre la placa se ha aplicado una fuerza que ha provocado una vibración de manera que los desplazamientos en un tiempo [math]t_0[/math] dado vienen dados por:: [math] \vec u(\rho,\theta)=\frac{\sin(\pi \theta/2)}{30\rho}\vec g_{\rho}. [/math] El código de Matlab para dibujar el campo de vectores en los puntos del mallado del sólido es el que mostramos a continuación:
h=0.1 % Intervalo de separación
u=1:h:2; % Intervalo de rho
v=0:h:2*pi+h; % Intervalo de theta
[uu,vv]=meshgrid(u,v);% Matriz de rho y theta
figure(1)
xx=uu.*cos(vv); % Parametrización
yy=uu.*sin(vv);
figure(1)
m=((sin((pi.*vv)./2).*cos(vv))./(30*uu));
n=((sin((pi.*vv)./2).*sin(vv))./(30*uu));
quiver(xx,yy,m,n); % Dibujo de la función
view4.1 Aplicación del desplazamiento
El desplazamiento que realiza la placa consiste en una serie de movimientos transversales en distintas direcciones. Como podemos ver en la siguiente imagen, obtenida con el código Matlab, el desplazamiento es prácticamente inapreciable.
subplot(1,2,1)
h=0.1; % Intervalo de separación
u=1:h:2; % Intervalo de rho
v=0:h:2*pi+h; % Intervalo de theta
[uu,vv]=meshgrid(u,v); % Matriz de rho y theta
figure(1)
xx=uu.*cos(vv); % Parametrización
yy=uu.*sin(vv);
mesh(xx,yy,0*xx) % Dibujo del mallado
axis([-3,3,-3,3]) % Región del dibujo
axis equal
subplot(1,2,2)
h=0.1 % Intervalo de separación
u=1:h:2; % Intervalo de rho
v=0:h:2*pi+h; % Intervalo de theta
[uu,vv]=meshgrid(u,v); % Matriz de rho y theta
figure(1)
xx=uu.*cos(vv); % Parametrización
yy=uu.*sin(vv);
figure(1)
m=((sin((pi.*vv)./2).*cos(vv))./(30*uu));
n=((sin((pi.*vv)./2).*sin(vv))./(30*uu));
mesh(m+xx,n+yy,xx*0) % Dibujo de las funciones
axis([-3,3,-3,3])
axis equal
view
5 Estudio de la divergencia
La divergencia es una medida del cambio de volumen local debido al desplazamiento. En este caso vamos a estudiarla tomando como volumen el de nuestra placa. En este caso, al considerar el campo [math]\vec u(\rho,\theta)=\frac{\sin(\pi \theta/2)}{30\rho}\vec g_{\rho}[/math], la divergencia la obtenemos con la siguiente expresión en coordenadas cilíndricas::
[math] \nabla \cdot \vec u = \frac{1}{\sqrt{|g|}} \frac{\partial}{\partial x^i} \left(\sqrt{|g|} u^i \right)= \frac{1}{\rho} \frac{\partial}{\partial\rho} ({\rho} \frac{\sin(\pi \theta/2)}{30\rho})= 0\ [/math]
En el estudio de la divergencia nos damos cuenta de que todos los puntos tienen la misma, ya que sobre nuestra placa es cero. Esto se debe a que la placa se ve sometida a desplazamientos transversales (como hemos visto en el apartado anterior) sin deformación de volumen, por lo que no habría un cambio del mismo.
6 Estudio del rotacional
En la corona circular de estudio donde los desplazamientos vienen representados por el vector [math] \vec u(\rho,\theta)=\frac{\sin(\pi \theta/2)}{30\rho}\vec g_{\rho} [/math], el rotacional muestra la tendencia de un campo vectorial a inducir rotación alrededor de un punto::
[math] \nabla\times\vec u = \frac{1}{\rho} \begin{vmatrix} \vec g_{\rho} & \vec g_{\theta} & \vec g_{z} \\ \frac{\partial}{\partial \rho} & \frac{\partial}{\partial \theta} & \frac{\partial}{\partial z} \\ \frac{\sin(\pi \theta/2)}{30\rho} & 0 & 0 \end{vmatrix} = \frac{1}{\rho} \frac{\partial}{\partial\theta} (\frac{\sin(\pi \theta/2)}{30\rho})\vec g_{z}= -\frac{\pi}{2}\frac{\cos(\pi \theta/2)}{30\rho^2}\vec g_{z} [/math]
Una vez calculado el rotacional, obtenemos su módulo::
[math]|\nabla × \vec u| =\frac{\pi \cos(\pi \theta/2)}{60\rho^2} [/math]
h=0.1; % Intervalo de separación
u=1:h:2; % Intervalo de rho
v=0:h:2*pi+h; % Intervalo de theta
[uu,vv]=meshgrid(u,v); % Matriz de rho y theta
xx=uu.*cos(vv); % Parametrización
yy=uu.*sin(vv);
figure(1)
f=abs(pi.*cos((pi.*vv)./2))./(60.*uu.^2); % Campo del rotacional
surf(xx,yy,f) % Líneas de nivel del rotacional
hold on
view(2)
h=0.1; % Intervalo de separación
u=1:h:2; % Intervalo de rho
v=0:h:2*pi+h; % Intervalo de theta
[uu,vv]=meshgrid(u,v); % Matriz de rho y theta
figure(1)
xx=uu.*cos(vv); % Parametrización
yy=uu.*sin(vv);
mesh(xx,yy,0*xx) % Dibujo de la malla
axis([-3,3,-3,3]) % Región del dibujo
view(2)
hold off
Matemáticamente, los puntos en los que hay un mayor rotacional son aquellos en los que [math]\cos(\pi\theta/2)=\pm 1 [/math] y [math]\rho=1[/math]. Gráficamente, estas zonas con un mayor valor del rotacional vienen representadas con un color rojizo.
7 Estudio del tensor de tensiones
Definamos [math]\epsilon(\vec u)=(\nabla \vec u + \nabla \vec u^t)/2[/math], la parte simétrica del tensor gradiente de [math]\vec u[/math], que se denomina tensor de deformaciones. En un medio elástico lineal, isótropo y homogéneo los desplazamientos permiten escribir el tensor de tensiones [math]\sigma_{ij}[/math] a través de la fórmula:: [math] \sigma_{ij}=\lambda \nabla \cdot \vec u \delta_{ij} + 2\mu \epsilon_{ij}, [/math] Siendo [math] \nabla \cdot \vec u = 0 [/math] y calculando el gradiente de [math]\vec u[/math] como:: [math]\nabla\vec u=u^i_{.j}\vec g_i \otimes \vec g^j[/math] donde [math]u^i_{.j}= \vec g^i \nabla\vec u \vec g_j=\frac{\partial u^i}{\partial x^j}+\Gamma^i_{kj}u^k.[/math] Podemos calcular el tensor de deformaciones como 1-contravariante 1-covariante:: [math](\epsilon^i_{.j})=(\nabla \vec u + \nabla \vec u^t)/2=\left(\begin{array}{ccc}-\frac{\sin(\pi\theta/2)}{30\rho^2}& \frac{\pi\cos(\pi\theta/2)}{120\rho}&0\\\frac{\pi\cos(\pi\theta/2)}{120\rho}& \frac{\sin(\pi\theta/2)}{30\rho^2}&0\\0&0&0\end{array}\right)[/math] En coordenadas 2-covariantes tendríamos:: [math]\epsilon(\vec u)=(\epsilon_{ij})=G(\epsilon^i_{.j})=\left(\begin{array}{ccc}-\frac{\sin(\pi\theta/2)}{30\rho^2}& \frac{\pi \cos(\pi\theta/2)}{120\rho}&0\\\frac{\pi\rho\cos(\pi\theta/2)}{120}& \frac{\sin(\pi\theta/2)}{30}&0\\0&0&0\end{array}\right)[/math]
Sustituyendo ahora los valores en la fórmula del tensor de tensiones,y tomando [math]\lambda=\mu=1[/math], siendo [math]\lambda[/math] y [math]\mu[/math] los conocidos como coeficientes de Lamé, obtenemos:: [math](\sigma_{ij})=2(\epsilon_{ij})=\left(\begin{array}{ccc}-\frac{\sin(\pi\theta/2)}{15\rho^2}& \frac{\pi \cos(\pi\theta/2)}{60\rho}&0\\\frac{\pi\rho\cos(\pi\theta/2)}{60}& \frac{\sin(\pi\theta/2)}{15}&0\\0&0&0\end{array}\right)[/math]
Las tensiones normales en la dirección [math]\vec g_{\rho}[/math] serían:: [math]\vec g_{\rho} \cdot \sigma \cdot \vec g_\rho =\frac{-sen(\pi \theta/2)}{15\rho^2} [/math]
Las tensiones normales en la dirección [math]\vec g_\theta/\rho[/math] serían:: [math]\vec g_\theta/\rho \cdot \sigma \cdot \vec g_\theta/\rho = \frac{\sin(\pi\theta/2)}{15\rho^2} [/math]
A continuación se muestra el código Matlab empleado para dibujar estas tensiones normales en ambas direcciones:
h=0.1; % Intervalo de separación
u=[1:0.1:2]; % Intervalo de rho [1,2]
v=[0:h:2*pi+h]; % Intervalo de theta [0,2*pi]
[uu,vv]=meshgrid(u,v); % Matrices de rho and theta
figure(1)
subplot(1,2,1)
xx=uu.*cos(vv); % Parametrización
yy=uu.*sin(vv);
f=((-sin(pi*vv./2))./(15*(xx.^2+yy.^2))); % Campo escalar
surf(xx,yy,f) % Dibujar el mallado
axis([-2,2,-2,2]) % Región seleccionada
axis equal
colorbar
title('Tensión normal en dirección g sub rho')
view(2)
subplot(1,2,2)
f=sin(pi.*vv/2)./(15.*uu.^2); % Campo escalar
surf(xx,yy,f) % Dibujar el mallado
axis([-2,2,-2,2]) % Región del dibujo
axis equal
colorbar
title('Tensión normal en dirección g sub theta')
view(2)7.1 Tensiones tangenciales
7.1.1 Respecto al plano ortogonal a [math]\vec g_\rho[/math] y respecto al plano ortogonal a [math]\vec g_\theta/\rho[/math]
Por un lado, vamos a analizar las tensiones tangenciales respecto al plano ortogonal a [math]\vec g_\rho[/math]:: [math]|\sigma \cdot \vec g_\rho-(\vec g_\rho \cdot \sigma \cdot \vec g_\rho) \vec g_\rho| = \frac{\pi\rho cos(\pi \theta/2)}{60} [/math]
Como podemos observar, en este caso, las tensiones son máximas cuando $\rho$ =2 y $\cos(\pi\theta/2)=\pm1$.
Por otro lado, analizamos las tensiones tangenciales respecto al plano ortogonal a [math]\vec g_\theta/\rho[/math]::
[math]|\sigma \cdot \vec g_\theta/\rho-(\vec g_\theta/\rho \cdot \sigma \cdot \vec g_\theta/\rho) \vec g_\theta/\rho|= \frac{\pi cos(\pi \theta/2)}{60 \rho^2}[/math]
En este otro caso, las tensiones son máximas cuando $\rho=1$ y $\cos(\pi\theta/2)=\pm1$, por tanto observamos que coinciden con los puntos en los cuales es máximo el rotacional.
En Matlab podemos dibujar estas tensiones mediante el siguiente código:
h=(0.1); % Intervalo de separación
u=[1:h:2]; % Intervalo de rho [1,2]
v=[0:h:2*pi+h]; % Intervalo de theta [0,2*pi]
[uu,vv]=meshgrid(u,v); % Matrices de rho y theta
figure(1)
subplot(1,2,1)
xx=uu.*cos(vv); % Parametrización
yy=uu.*sin(vv);
f=abs((pi.*uu.*cos(pi*vv./2)))./(60); % Campo escalar
surf(xx,yy,f) % Dibujar el mallado
axis([-2,2,-2,2]) % Región para dibujar
axis equal
colorbar
view(2)
subplot(1,2,2)
f=abs((pi*(cos(pi*vv./2)))./(60*uu.^2)); % Campo escalar
surf(xx,yy,f) % Dibujar el mallado
axis([-2,2,-2,2]) % Región para dibujar
axis equal
colorbar
view(2)
