Visualización de campos escalares y vectoriales en elasticidad (grupo 9 A)

}} Consideramos una placa plana (en 2 dimensiones) 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 instante [math]t_0[/math] vienen dados por: [math] \vec u(\rho,\theta)=\frac{\sin(\pi \theta)}{20\rho^2}\vec g_{\theta} [/math]

1 Mallado del sólido

Nuestro mallado se ha centrado en la representación de los puntos interiores de una placa sólida, con forma de corona circular. Tomando como paso de muestreo [math]h=\frac{1}{10}[/math] para las variables [math]x[/math] e [math]y[/math], los comandos en Matlab para la visualización de nuestra placa son los siguientes:

Placa mallada con forma de corona circular.

h=0.1;                 % Paso de muestreo.
u=1:h:2;               % Intervalo [1,2].
v=0:h:2*pi+h;          % Intervalo [0,2*pi].
[uu,vv]=meshgrid(u,v); % Matrices.
xx=uu.*cos(vv);        % Parametrización X.
yy=uu.*sin(vv);        % Parametrización Y.
subplot(1,2,1);        % Muestra varias imágenes. 1º Imagen.
mesh(xx,yy,0*xx)       % Mallado.
axis([-3,3,-3,3])      % Selecciona la regíon a dibujar.
view(2)                % Ver imagen desde arriba.
subplot(1,2,2);        % Muestra varias imágenes. 2º Imagen.
mesh(xx,yy,0*xx);      % Visión general.

2 Temperatura del sólido

Teniendo en cuenta que la temperatura es un campo escalar y que en nuestro caso viene definida por [math]T(x,y)=e^{-y}[/math],vamos a tratar de observar su comportamiento a lo largo de nuestra corona circular plana. Nuestra función depende únicamente de la segunda variable, esto quiere decir que la temperatura se mantendrá constante para cada valor [math]y_{0}[/math]. Es decir, cada recta paralela al eje de abcisas tendrá una temperatura [math]T(x,y_{0})=e^{-y_{0}}=cte[/math]. Además, por el signo negativo del exponente, la [math]T[/math] va aumentando a medida que disminuye la coordenada [math]y[/math], ubicándose el foco de calor en la zona más inferior de la placa [math](0,-2)[/math] . Así pues, se observa claramente en la figura como la temperatura va descendiendo a medida que nos vamos alejando del foco.

[math]T[/math] a lo largo de la placa.

f=exp(-yy); 	          % Campo escalar de la Temperatura.
subplot(1,2,1);           % Muestra varias imágenes. 1º Imagen.
surf(xx,yy,f);            % Visualización 2D.
view(2)                   % Ver imagen desde arriba.
axis([-3,3,-3,3])         % Selecciona la regíon a representar en 2D.  
subplot(1,2,2);           % Muestra varias imágenes. 2º Imagen.
surf(xx,yy,f); colorbar;  % Visualización 3D más leyenda en color.

2.1 Variación de la Temperatura ([math]\nabla T[/math])

Para estudiar la variación del campo escalar [math]T[/math], haremos uso del gradiente [math]\nabla T[/math] . Se puede observar, por la propia definición del [math]\nabla[/math], que dicho campo vectorial es perpendicular a las curvas de nivel de nuestra placa plana. A medida que nos acercamos a temperaturas más elevadas, el módulo de [math]\vec\nabla[/math] en cada punto va siendo mayor.

Campo [math]\nabla T[/math].

Grad=-exp(-yy);	                   % Gradiente del campo.
hold on                            % Inicio superposición de gráficos
quiver(xx,yy,0*xx,Grad);           % Representación de los vectores gradientes del campo escalar.
contour(xx,yy,f,20)                % Define 20 líneas de nivel.;
hold off                           % Fin superposición de gráficos.

3 Campo de desplazamientos

Nuestra superficie sufre un desplazamiento en un instante [math]t_{0}[/math] debido a una percusión. Dicho campo viene dado por el siguiente vector [math]\vec u[/math]: [math]\vec u(\rho,\theta)=\frac{\sin(\pi \theta)}{20\rho^2}\vec g_{\theta}=\frac{-\sin(\theta)\sin(\pi \theta)}{20\rho}\vec i + \frac{\cos(\theta)\sin(\pi \theta)}{20\rho}\vec j[/math]

Aplicando este campo a cada punto del solido obtendremos los vectores de desplazamiento de los mismos ante dicha percusión.

Vectores [math]\vec u[/math] de desplazamiento del sólido.

3.1 Resultado de la percusión del solido

En la figura se muestran los vectores de desplazamiento de cada uno de los puntos del sólido según el campo vectorial [math]\vec u(\rho,\theta)[/math].

i=1; j=1;	                                            % Definición de las variables para el bucle.
Ux=[];                                                      % Matrices de vectores componentes x de u(rho,theta).
Uy=[];                                                      % Matrices de vectores componentes y de u(rho,theta).
for i=1:64	                                            % Inicio bucle para calcular el vector de cada punto. De 1 a 64 para v.
  for j=1:11                                                % Inicio bucle para calcular el vector de cada punto. De 1 a 11 para u.
    Ux(i,j)=-sin(vv(i,j))*(sin(pi*vv(i,j)))/(20*uu(i,j));   % Proyección de u(rho,theta) sobre x.
    Uy(i,j)=cos(vv(i,j))*(sin(pi*vv(i,j)))/(20*uu(i,j));    % Proyección de u(rho,theta) sobre y.
    j=j+1;                                                  % Bucle j.
    end                                                     % Fin del bucle j.
  i=i+1;                                                    % Bucle i.
end                                                         % Fin del bucle i.
quiver(xx,yy,Ux,Uy);	                                    % Representación del campo vectorial de desplazamientos.

Podemos apreciar que el vector [math]\vec u[/math] que hace que el solido sufra una percusión, tiene un modulo muy pequeño, lo que dificulta el estudio y visión de la deformación a la que es sometido. Para apreciar mejor el efecto de la percusión en el sólido, procedemos a incrementar el tamaño de la deformación por [math]10[/math]. En la figura adjunta se puede apreciar claramente las distintas fases de la deformación del solido: el sólido en reposo, el sólido deformado debido a [math]\vec u[/math] y el sólido deformado debido a [math]10\cdot\vec u[/math]

Grados de deformación según la variación del módulo del vector [math]\vec u[/math].

subplot(1,3,1); 		                            % Muestra varias imágenes. 1º Imagen: sólido sin deformar.
mesh(xx,yy,0*xx)                                            % Mallado.
axis([-3,3,-3,3])                                           % Selecciona la regíon a dibujar.
view(2)                                                     % Ver imagen desde arriba.
i=1; j=1;	                                            % Definición de las variables para el bucle.
Ux=[];                                                      % Matrices de vectores componentes x de u(rho,theta).
Uy=[];                                                      % Matrices de vectores componentes y de u(rho,theta).
for i=1:64	                                            % Inicio bucle para calcular el vector de cada punto. De 1 a 64 para v.
  for j=1:11                                                % Inicio bucle para calcular el vector de cada punto. De 1 a 11 para u.                                              
    Ux(i,j)=-sin(vv(i,j))*(sin(pi*vv(i,j)))/(20*uu(i,j));   % Proyección de u(rho,theta) sobre x.
    Uy(i,j)=cos(vv(i,j))*(sin(pi*vv(i,j)))/(20*uu(i,j));    % Proyección de u(rho,theta) sobre y.
    Uxx(i,j)=xx(i,j)+Ux(i,j);                               % Proyección vector posición tras la deformación sobre x.
    Uyy(i,j)=yy(i,j)+Uy(i,j);                               % Proyección vector posición tras la deformación sobre y.
    j=j+1;                                                  % Bucle j.
    end                                                     % Fin del bucle j.
  i=i+1;                                                    % Bucle i.
end                                                         % Fin del bucle i.
subplot(1,3,2);                                             % Muestra varias imágenes. 2º Imagen: sólido deformado por u(rho,theta).
mesh(Uxx,Uyy,0*xx);                                         % Mallado.
axis([-3,3,-3,3])                                           % Selecciona la regíon a dibujar.
view(2)                                                     % Ver imagen desde arriba.   
i=1; j=1;	                                            % Definición de las variables para el bucle.
for i=1:64	                                            % Inicio bucle para calcular el vector de cada punto. De 1 a 64 para v.
  for j=1:11                                                % Inicio bucle para calcular el vector de cada punto. De 1 a 11 para u.                                            
    Ux(i,j)=-sin(vv(i,j))*(sin(pi*vv(i,j)))/(2*uu(i,j));    % Proyección de 10u(rho,theta) sobre x.
    Uy(i,j)=cos(vv(i,j))*(sin(pi*vv(i,j)))/(2*uu(i,j));     % Proyección de 10u(rho,theta) sobre y.
    Uxxx(i,j)=xx(i,j)+Ux(i,j);                              % Proyección vector posición tras la deformación ampliada sobre x.
    Uyyy(i,j)=yy(i,j)+Uy(i,j);                              % Proyección vector posición tras la deformación ampliada sobre y.
    j=j+1;                                                  % Bucle j.
    end                                                     % Fin del bucle j.
  i=i+1;                                                    % Bucle i.
end                                                         % Fin del bucle i.
subplot(1,3,3);                                             % Muestra varias imágenes. 3º Imagen: sólido deformado por 10u(rho,theta).
mesh(Uxxx,Uyyy,0*xx);                                       % Mallado.
axis([-3,3,-3,3])                                           % Selecciona la regíon a dibujar.
view(2)                                                     % Ver imagen desde arriba.

3.2 Cambio de volumen local debido al desplazamiento

Para estudiar la variación del volumen después de la percusión, hacemos uso de la divergencia del vector [math]\vec u(\rho,\theta)[/math]. Aplicando este operador al campo vectorial [math]\vec u[/math], obtenemos la función escalar divergencia de [math]\vec u[/math]: [math]\nabla\cdot\vec u=\frac{\pi\cos(\pi \theta)}{20\rho^2}[/math] Aplicando la divergencia observamos las variaciones de volumen en cada punto del sólido.

Variacion del volumen según [math]\nabla\cdot\vec u[/math].

Divergencia=[];                                                % Definición de la matriz Divergencia.
i=1;j=1;	                                               % Definición de las variables para el bucle.
for i= 1:64	                                               % Inicio bucle para calcular el vector de cada punto. De 1 a 64 para v.
  for j = 1:11                                                 % Inicio bucle para calcular el vector de cada punto. De 1 a 11 para u.
    Divergencia(i,j)=pi*(cos(pi.*vv(i,j)))/(20*(uu(i,j).^2));  % Cálculo de la divergencia del campo u.
    j = j+1;                                                   % Bucle j.
    end                                                        % Fin del bucle j.
  i= i+1;                                                      % Bucle i.
end                                                            % Fin del bucle i.
hold on                                                        % Inicio superposición de gráficos.
contour(xx,yy,Divergencia,20); colorbar;                       % Define 20 líneas de nivel más leyenda en color.
hold off                                                       % Fin superposición de gráficos.

Se observa que los puntos con el máximo aumento y disminución de volumen se encuentran en la circunferencia interior de la corona, donde observamos mayor número de líneas de nivel.

3.3 Tendencia a la rotación

El rotacional se define como un operador vectorial que muestra la tendencia de un campo vectorial a inducir rotación alrededor de un punto. En nuestro campo de desplazamiento [math]\vec u[/math], el rotacional será: [math]\frac{\sin(\pi \theta)}{20\rho^2}\vec g_{\theta}=\rho^2\frac{\sin(\pi \theta)}{20\rho^2}\vec g^{\theta}=\frac{\sin(\pi \theta)}{20}\vec g^{\theta}[/math]: [math]\nabla \times \vec u=\frac{1}{\rho}\begin{bmatrix} \vec g_{\rho} & \vec g_{\theta} & \vec g_{z} \\ \frac{\delta}{\delta \rho} & \frac{\delta}{\delta \theta} & \frac{\delta}{\delta z} \\ 0 & \frac{\sin(\pi \theta)}{20} & 0 \end{bmatrix}=\vec {0}\ [/math]:

Por tanto, la tendencia a la rotación en cualquier punto del sólido debido a la acción de la percusión es nula.

3.4 Tensiones debidas a la perturbación

Se define como [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] donde [math]\lambda[/math] y [math]\mu[/math] son los conocidos como coeficientes de Lamé que dependen de las propiedades elásticas de cada material. Tomando [math]\lambda=\mu=1[/math],El resultado de la parte simetrica del gradiente de [math]\vec u[/math]: [math]2\mu \epsilon_{ij}=\begin{bmatrix} 0 & \frac{-sin(\pi \theta)}{20\rho}(\frac{1}{\rho^2}+1) & 0 \\ \frac{-sin(\pi \theta)}{20\rho}(\frac{1}{\rho^2}+1) & \frac{\pi\cos(\pi \theta)}{10\rho^2} & 0 \\ 0 & 0 & 0 \end{bmatrix}[/math]: Calculamos [math]\lambda \nabla \cdot \vec u \delta_{ij}[/math]: [math]\lambda \nabla \cdot \vec u \delta_{ij}=\begin{bmatrix} \frac{\pi\cos(\pi \theta)}{20\rho^2} & 0 & 0 \\ 0 & \frac{\pi\cos(\pi \theta)}{20\rho^2} & 0 \\ 0 & 0 & \frac{\pi\cos(\pi \theta)}{20\rho^2} \end{bmatrix}[/math]: Sumamos ambas y obtenemos [math]\sigma_{ij}[/math] en forma de la siguiente matriz: [math]\sigma_{ij}=\begin{bmatrix} \frac{\pi\cos(\pi \theta)}{20\rho^2} &\frac{-sin(\pi \theta)}{20\rho}(\frac{1}{\rho^2}+1) & 0 \\ \frac{-sin(\pi \theta)}{20\rho}(\frac{1}{\rho^2}+1) & \frac{3\pi\cos(\pi \theta)}{20\rho^2} & 0 \\ 0 & 0 & \frac{\pi\cos(\pi \theta)}{20\rho^2} \end{bmatrix}[/math]:

También calculamos las tensiones normales en la dirección [math]\vec g_{\rho}[/math] según [math]\vec g_{\rho} \cdot \sigma \cdot \vec g_\rho[/math] y en la dirección [math]\vec g_\theta/\rho[/math] según [math]\vec g_\theta/\rho \cdot \sigma \cdot \vec g_\theta/\rho[/math] obteniendo:

[math]\vec g_{\rho} \cdot \sigma \cdot \vec g_\rho = \frac{\pi\cos(\pi\theta)}{20\rho^2}\\ \vec g_\theta/\rho \cdot \sigma \cdot \vec g_\theta/\rho =\frac{3\pi\cos(\pi\theta)}{20\rho^2}[/math]

Tensiones normales en la direccion de [math]\vec g_\rho[/math] y [math]\vec g_\theta /\rho[/math].

rhosigmarho=(pi*cos(pi.*vv))./(20.*(uu.^2));                        % Tensiones normales en la dirección del eje grho.
thetasigmatheta=(3*pi*cos(pi.*vv))./(20.*(uu.^2));                  % Tensiones normales en la dirección del eje gtheta/rho.
subplot(2,2,1);                                                     % Muestra varias imágenes. 1º Imagen. Vista superior tensiones normales según eje grho.
mesh(xx,yy,rhosigmarho); colorbar                                   % Mallado más leyenda en color.
view(2)                                                             % Ver imagen desde arriba.
subplot(2,2,3);                                                     % Muestra varias imágenes. 3º Imagen. Vista superior tensiones normales según eje gtheta/rho.
mesh(xx,yy,thetasigmatheta); colorbar                               % Mallado más leyenda en color.
view(2)                                                             % Ver imagen desde arriba.
subplot(2,2,2);                                                     % Muestra varias imágenes. 2º Imagen. Tensiones normales según eje grho.
mesh(xx,yy,rhosigmarho);                                            % Mallado.
subplot(2,2,4);                                                     % Muestra varias imágenes. 4º Imagen. Tensiones normales según eje gtheta/rho.
mesh(xx,yy,thetasigmatheta);                                        % Mallado.

Tal y como se observa en la figura, las tensiones normales en la dirección [math]\vec g_{\rho}[/math] y en las de dirección [math]\vec g_\theta/\rho[/math] son mayores en las zonas que presentan una mayor deformación como consecuencia de la percusión.

3.4.1 Tensiones tangenciales respecto al plano ortogonal a [math]\vec g_\rho[/math]

Las tensiones tangenciales al plano ortogonal respecto de [math]\vec g_{\rho}[/math] son [math]|\sigma \cdot \vec g_\rho-(\vec g_\rho \cdot \sigma \cdot \vec g_\rho) \vec g_\rho|[/math]. Su representación nos indicará en qué puntos las tensiones son mayores, y dado que el estudio se realiza en un medio elástico lineal, isótropo y homogéneo, en ellos también serán las máximas. Si operamos adecuadamente: [math]|\sigma \cdot \vec g_\rho-(\vec g_\rho \cdot \sigma \cdot \vec g_\rho) \vec g_\rho|=|\frac{\pi\cos(\pi \theta)}{20\rho^2}\vec g_\rho - \frac{sin(\pi \theta)}{20\rho}(\frac{1}{\rho^2}+1)\vec g_\theta - \frac{\pi\cos(\pi \theta)}{20\rho^2}\vec g_\rho|= \frac{sin(\pi \theta)}{20\rho}(\frac{1}{\rho^2}+1)\vec g_\theta[/math]

Distribución de las tensiones normales y tangenciales respecto a [math]\vec g_{\rho}[/math]

rhosigmarho=(pi*cos(pi.*vv))./(20.*(uu.^2)); 	               % Tensiones normales en la dirección del eje grho.
subplot(1,2,1);                                                % Muestra varias imágenes. 1º Imagen. Tensiones normales según eje grho.
mesh(xx,yy,rhosigmarho);                                       % Mallado.
view(2)                                                        % Ver imagen desde arriba.
subplot(1,2,2); 	                                       % Muestra varias imágenes. 1º Imagen. Tensiones tangenciales respecto al plano ortogonal a grho.
tangencialrho=(sin(pi.*vv)).*(1./(20.*uu)+1./(20.*uu.^3))      % Tensiones tangenciales respecto al plano ortogonal a grho.
mesh(xx,yy,tangencialrho);                                     % Mallado.
view(2)                                                        % Ver imagen desde arriba.

Se puede observar la ortogonalidad de las tensiones tangenciales y normales de [math]\vec g_{\rho}[/math]. En los puntos con menor desplazamiento, las tensiones tangenciales alcanzan su valor máximo absoluto.

3.4.2 Tensiones tangenciales respecto al plano ortogonal a [math]\vec g_\theta/\rho[/math]

Las tensiones tangenciales al plano ortogonal respecto de [math]\vec g_\theta/\rho[/math] son [math]|\sigma \cdot \vec g_\theta/\rho-(\vec g_\theta/\rho \cdot \sigma \cdot \vec g_\theta/\rho) \vec g_\theta/\rho|[/math]. Su representación nos indicará en qué puntos las tensiones son mayores, y dado que el estudio se realiza en un medio elástico lineal, isótropo y homogéneo, en ellos también serán las máximas. Si operamos adecuadamente: [math]|\sigma \cdot \vec g_\theta/\rho-(\vec g_\theta/\rho \cdot \sigma \cdot \vec g_\theta/\rho) \vec g_\theta/\rho| = |\frac{-sin(\pi \theta)}{20\rho}(\frac{1}{\rho^2}+1)\vec g_\rho + \frac{3\pi\cos(\pi \theta)}{20\rho^2}\vec g_\theta - \frac{3\pi\cos(\pi \theta)}{20\rho^2}\vec g_\theta|=\frac{sin(\pi \theta)}{20\rho}(\frac{1}{\rho^2}+1)\vec g_\rho[/math].

Distribución de las tensiones normales y tangenciales respecto a [math]\vec g_\theta /\rho[/math]

thetasigmatheta=(3*pi*cos(pi.*vv))./(20.*(uu.^2));                      % Tensiones normales en la dirección del eje gtheta/rho.
subplot(1,2,1);                                                         % Muestra varias imágenes. 1º Imagen. Tensiones normales según eje gtheta/rho.
mesh(xx,yy,thetasigmatheta);                                            % Mallado.
view(2)                                                                 % Ver imagen desde arriba.
subplot(1,2,2); 	                                                % Muestra varias imágenes. 2º Imagen. Tensiones tangenciales respecto al plano ortogonal a gtheta/rho.
tangencialtheta=(sin(pi.*vv)).*(1./(20.*uu)+1./(20.*uu.^3));            % Tensiones tangenciales respecto al plano ortogonal a gtheta/rho.
mesh(xx,yy,tangencialtheta);                                            % Mallado.
view(2)                                                                 % Ver imagen desde arriba.

Se puede observar la ortogonalidad de las tensiones tangenciales y normales de [math]\vec g_\theta /\rho[/math]. En los puntos con menor desplazamiento, las tensiones tangenciales alcanzan su valor máximo absoluto. Evidentemente, queda demostrado que las tensiones tangenciales del sólido se pueden describir mediante una base de vectores ortonormales [ [math]\vec g_\rho[/math] , [math]\vec g_\theta /\rho[/math] ]