Elasticidad en Campos escalares y vectoriales (Grupo C3)

1 Introducción

Consideramos una placa plana (en dimensión 2) que ocupa el anillo comprendido entre las circunferencias centradas en el origen de radios 2 y 1. En ella 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] 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]

La representación de los puntos interiores de la placa sólida nos queda en forma de corona circular como podemos observar:

Mallado de la placa con paso de muestreo h=1/10

Dicha representación se ha obtenido a través de Matlab mediante el código siguiente:

u=1:0.1:2;               % mallado del intervalo [1,2] con tamaño h= 0.1
v=0:0.1:2*pi+0.1;            % mallado del intervalo [0,2*pi] con tamaño h=0.1
[uu,vv]=meshgrid(u,v); % matrices de coordenadas 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 gráfico
view(2)                % punto de vista

En los siguientes apartados se utilizará los mallados y parametrizaciones de este primer código(líneas de 1 a 6).

2 Influencia de un foco de calor

La temperatura proviene de un foco de calor muy concentrado en el origen y viene dada por un campo escalar:[math]T(\rho,\theta)=-\log(\rho+0.1)[/math] Como se puede observar la temperatura varía únicamente con el parámetro [math]\rho[/math], por lo tanto para un mismo valor de éste la temperatura es constante sea cual sea el valor de [math]\theta[/math] comprendido entre 0 y [math]2Π[/math]. Vemos que cuanto mayor es el valor de [math]\rho[/math] la temperatura va disminuyendo debido al sentido negativo de la fórmula; luego los puntos de máxima temperatura serán aquellos que estén más cerca al origen pues están mas cerca al foco.

Temperatura a lo largo del sólido

%campo escalar T
figure(2)
T=-log(0.1+uu); % Campo escalar
subplot(1,2,1),surf(xx,yy,T); % mallado 
axis([-3,3,-3,3])      % región
view(2)                % punto de vista
subplot(1,2,2),surf(xx,yy,T)

A partir de la definición de gradiente se puede observar mediante las gráficas la dirección en la cual la temperatura varía mas rápido. El ritmo de variación de la temperatura se establece a partir del módulo del gradiente y por lo tanto será menor en los puntos más próximos al foco.

[math]\nabla T[/math] =-[math]\frac{\partial T}{\partial \rho}\vec{g}^{\rho}[/math]= -[math]\frac{1}{0,1+\rho}\vec{g}^{\rho}\Rightarrow |\nabla T|=\frac{1}{0.1+\rho}[/math]

Como se observa en la figura, el gradiente es perpendicular a la superficie y por lo tanto es perpendicular a las curvas de nivel.

Campo de temperatura y curvas de nivel

%gradiente de T y curvas de nivel
figure(3)
Tx=-1./(0.1+uu).*cos(vv);                 % componente x del gradiente 
Ty=-1./(0.1+uu).*sin(yy);                 % componente y del gradiente
subplot(1,3,1),quiver(xx,yy,Tx,Ty)     % dibujo del gradiente
axis([-3,3,-3,3])      % Región del grafico
subplot(1,3,2),contour(xx,yy,T,10)
axis([-3,3,-3,3])
subplot(1,3,3)
hold on
quiver(xx,yy,Tx,Ty)
contour(xx,yy,T,10)
axis([-3,3,-3,3])
hold off
view(2)% punto de vista

3 Deformaciones que sufre el sólido

Las deformaciones que sufre el sólido viene dado por un campo vectorial [math] \vec u(\rho,\theta)[/math] que nos indica el desplazamiento que sufre cada punto del sólido debido a la vibración producida al aplicar una fuerza sobre éste::

[math] \vec u(\rho,\theta)=\frac{\sin(\pi \theta/2)}{30\rho}\vec g_{\rho}. [/math]

El campo vectorial en los puntos de la placa viene representado en la siguiente imagen:

Campo de desplazamiento de cada punto del sólido

Se puede apreciar que el campo no es uniforme en todos los puntos de la placa plana ya que tanto la dirección como el módulo varía de unos puntos a otros del anillo circular.

Se ha obtenido con el siguiente código:

%campo vectorial u
figure(4)
ux=sin(pi.*vv./2)./(30.*uu).*cos(vv);      % componente x de u 
uy=sin(pi.*vv./2)./(30.*uu).*sin(vv);      % componente y de u
quiver(xx,yy,ux,uy)     % dibujo del gradiente
axis([-3,3,-3,3])      % Región del grafico

Para poder apreciar mejor este desplazamiento vamos a mostrar dos imágenes: en la primera se muestra el sólido antes del desplazamiento y en la segunda una vez desplazado cada punto del sólido. Como podemos apreciar la deformación del sólido no es muy notable debido a que el módulo del desplazamiento es muy pequeño (muy próximo a 1/30)

Mallado de la placa antes y despues del desplazamiento

El código de matlab es el siguiente:

%desplazamiento del sólido
figure(5)
xd=xx+ux; %componente x del sólido desplazado
yd=yy+uy; %componente y del sólido desplazado
subplot(1,2,1), mesh(xx,yy,0*xx)
axis([-2.2,2.2,-2.2,2.2])      % región del gráfico
view(2)                % punto de vista
subplot(1,2,2), mesh(xd,yd,0*xd)
axis([-2.2,2.2,-2.2,2.2])      % región del gráfico
view(2)                % punto de vista

3.1 Cambio de volumen local del sólido

A partir del operador divergencia podemos establecer el cambio de volumen local como consecuencia del campo vectorial [math]\vec u[/math]. Conocido el campo vectorial [math]\vec u[/math] obtenemos su divergencia:: [math]\nabla \cdot u=\frac{1}{\sqrt{g}}\frac{\partial}{\partial x^i}(\sqrt{g}u^i)=\frac{1}{\rho}\left(\frac{\partial}{\partial\rho}\left(\frac{\sin(\pi \theta/2}{30}\right)\right)=0[/math]

Por tanto podemos decir que el cambio de volumen local es nulo y el cuerpo no se expande ni se comprime en ninguna dirección siendo un sólido incompresible.

3.2 Velocidad de giro

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}} \left|\begin{array}{ccc}\vec g_{\rho} & \vec g_{\theta} & \vec g_{z}\\\frac{\partial}{\partial \rho}& \frac{\partial}{\partial \theta}&\frac{\partial}{\partial z}\\u_{\rho}& u_{\theta}& u_z\end{array}\right|= \frac{1}{\rho} \left|\begin{array}{ccc}\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{array}\right|=-\frac{1}{30\rho^2}\left(\frac{\pi}{2}\cos(\pi\theta/2)\right) \vec g_{z} [/math]

Entonces [math] |\nabla \times\vec u|=\frac{1}{30\rho^2}\left(\frac{\pi}{2}|\cos(\pi\theta/2)|\right) [/math] y los puntos donde inducirá mayor rotación es cuando [math]\cos(\pi\theta/2)=\pm 1 [/math] y [math]\rho=1[/math].

El código de matlab es:

Módulo del rotacional

%módulo del rotacional
figure(6)
rot=abs(1./(30*uu.^2)*(pi/2).*cos(pi*vv/2));
surf(xx,yy,rot)   % mallado
axis([-3,3,-3,3])      % región
view(2)                % punto de vista

4 Tensiones normales

El tensor de deformaciones es la parte simétrica del tensor gradiente de [math]\vec u: \epsilon(\vec u)=\frac{\nabla \vec u +\nabla \vec u^t}{2}[/math]. En un medio elástico lineal, isótropo y homogéneo los desplazamientos permiten escribir el tensor de tensiones [math] \sigma_{ij}[/math] como: [math]\sigma_{ij}=\lambda \nabla \cdot \vec u\delta_{ij}+2\mu \epsilon_{ij}[/math] donde [math] \lambda[/math] y [math]\mu[/math] son los coeficientes de Lamé que dependen de las propiedades elásticas de cada material.

Calculamos el gradiente de [math]\vec u[/math]: [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]

Como [math]u^{\rho}=\frac{\sin(\pi \theta/2)}{30\rho}, u^{\theta}=u^z=0[/math] y teniendo en cuenta los símbolos de Christoffel para coordenadas cilíndricas se obtiene

[math]u^{\rho}_{.\rho}=\frac{\partial u^{\rho}}{\partial \rho}=-\frac{\sin(\pi\theta/2)}{30\rho^2}, \qquad u^{\rho}_{.\theta}=\frac{\partial u^{\rho}}{\partial \theta}=\frac{\pi\sin(\pi\theta/2)}{60\rho}, \qquad u^{\rho}_{.z}=0 [/math]

[math]u^{\theta}_{.\rho}=0, \qquad u^{\theta}_{.\theta}=\Gamma^{\theta}_{\rho\theta}u^{\theta}=\frac{\sin(\pi\theta/2)}{30\rho^2},\qquad u^{\theta}_z=0,[/math]

[math]u^z_{.\rho}=u^z_{.\theta}=u^z_{.z}=0.[/math]

Entonces el tensor de deformaciones en coordenadas contra-covas es

[math](\epsilon^i_{.j})=\left(\begin{array}{ccc}-\frac{\sin(\pi\theta/2)}{30\rho^2}& \frac{\pi\sin(\pi\theta/2)}{120\rho}&0\\\frac{\pi\sin(\pi\theta/2)}{120\rho}& \frac{\sin(\pi\theta/2)}{30\rho^2}&0\\0&0&0\end{array}\right)[/math]

y en coordenadas 2-covas

[math](\epsilon_{ij})=G(\epsilon^i_{.j})=\left(\begin{array}{ccc}-\frac{\sin(\pi\theta/2)}{30\rho^2}& \frac{\pi \rho\sin(\pi\theta/2)}{120}&0\\\frac{\pi\sin(\pi\theta/2)}{120\rho}& \frac{\sin(\pi\theta/2)}{30}&0\\0&0&0\end{array}\right)[/math]

El tensor de tensiones tomando [math] \lambda=\mu=1[/math] y teniendo en cuenta que [math]\nabla \cdot \vec u =0 [/math] es

[math](\sigma_{ij})=2(\epsilon_{ij})=\left(\begin{array}{ccc}-\frac{\sin(\pi\theta/2)}{15\rho^2}& \frac{\pi \rho\sin(\pi\theta/2)}{60}&0\\\frac{\pi\sin(\pi\theta/2)}{60\rho}& \frac{\sin(\pi\theta/2)}{15}&0\\0&0&0\end{array}\right)[/math]

Por tanto las tensiones normales en la dirección $\vec g_{\rho}$ son $\vec g_{\rho}\cdot \sigma\cdot \vec g_{\rho}=\sigma_{11}=-\frac{\sin(\pi\theta/2)}{15\rho^2}$.

Las tensiones normales en la dirección $\vec g_{\theta}/\rho$ son $\vec g_{\theta}/\rho \cdot \sigma\cdot\vec g_{\theta}/\rho=\frac{\sigma_{22}}{\rho^2}=\frac{\sin(\pi\theta/2)}{15\rho^2} $ç

Tensiones normales

El código en matlab es:

%Tensiones normales
figure(7)
Ep=-sin(pi.*vv/2)./(15.*uu.^2); % Campo escalar
Ett=sin(pi.*vv/2)./(15.*uu.^2); % Campo escalar
subplot(1,2,1);surf(xx,yy,Ep)          % mallado
axis([-3,3,-3,3])      % región
view(2)                % punto de vista
subplot(1,2,2);surf(xx,yy,Ett)          % mallado2
axis([-3,3,-3,3])      % región
view(2)                % punto de vista

5 Tensiones tangenciales

Las tensiones tangenciales respecto al plano ortogonal a $\vec g_\rho$ vienen dadas por: $|\sigma\cdot \vec g_\rho-(\vec g_\rho\cdot\sigma\cdot\vec g_\rho)\vec g_\rho|$.

$\sigma\cdot\vec g_\rho=-\frac{\sin(\pi\theta/2)}{15\rho^2}\vec g^\rho+\rho\frac{\cos(\pi\theta/2)}{60}\vec g^\theta=-\frac{\sin(\pi\theta/2)}{15\rho^2}\vec g_\rho+\frac{\cos(\pi\theta/2)}{60\rho}\vec g_\theta$

$(\vec g_\rho\cdot\sigma\cdot\vec g_\rho)\vec g_\rho=-\frac{\sin(\pi\theta/2)}{15\rho^2}\vec g_\rho$

Entonces $|\sigma\cdot \vec g_\rho-(\vec g_\rho\cdot\sigma\cdot\vec g_\rho)\vec g_\rho|=\pi\left|\frac{\cos(\pi\theta/2)}{60}\vec g_\theta\right|=\pi\frac{|\cos(\pi\theta/2)|}{60}$

Las tensiones solo dependen de $\theta$ y son máximas cuando $\cos(\pi\theta/2)=\pm1$ que son los puntos de menor deformación al aplicarle el campo vectorial $\vec u$.

Las tensiones tangenciales respecto al plano ortogonal a $\vec g_\theta/\rho$ vienen dadas por: $|\sigma\cdot \vec g_\theta/\rho-(\vec g_\theta/\rho\cdot\sigma\cdot\vec g_\theta/\rho)\vec g_\theta/\rho|$.

$\sigma\cdot\vec g_\theta/\rho=\frac{\pi\cos(\pi\theta/2)}{60\rho}\vec g^\rho+\frac{\sin(\pi\theta/2)}{15\rho}\vec g^\theta=\frac{\pi\cos(\pi\theta/2)}{60\rho}\vec g_\rho+\frac{\sin(\pi\theta/2)}{15\rho^3}\vec g_\theta$

$(\vec g_\theta/\rho\cdot\sigma\cdot\vec g_\theta/\rho)\vec g_\theta/\rho=\frac{\sin(\pi\theta/2)}{15\rho^3}\vec g_\theta$

Entonces $|\sigma\cdot \vec g_\theta/\rho-(\vec g_\theta/\rho\cdot\sigma\cdot\vec g_\theta/\rho)\vec g_\theta/\rho|=\left|\frac{\pi\cos(\pi\theta/2)}{60\rho}\vec g_\rho\right|=\frac{\pi|\cos(\pi\theta/2)|}{60\rho}$

Las tensiones son máximas cuando $\cos(\pi\theta/2)=\pm1$ y $\rho=1$ que coinciden con los puntos en los que el rotacional es máximo.

Tensiones tangenciales

El código en matlab es:

%Tensiones tangenciales
figure(8)
subplot(1,2,1),ttp=abs(pi.*cos(pi.*vv./2))./(60.*uu); % Campo escalar
surf(xx,yy,ttp)          % mallado
axis([-3,3,-3,3])      % región
subplot(1,2,2),ttp=abs(pi.*cos(pi.*vv./2))./(60.*uu); % Campo escalar
surf(xx,yy,ttp)          % mallado
axis([-3,3,-3,3])      % región
view(2)                % punto de vista