Visualización de campos escalares y vectoriales en elasticidad. Grupo C4

Trabajo realizado por estudiantes
Título	Visualización de campos escalares y vectoriales en elasticidad. Grupo C4
Asignatura	Teoría de Campos
Curso	2013-14
Autores	Alberto Rojas Rivero, Eduardo Bonet García, Fernando Sancha Domínguez, Carlos Fernández Bermejo, Ricardo Mazón Cabrera, Emilio Valero Muñoz-Rojas
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

Consideramos una placa circular comprendida entre las circunferencias de radio 1 y radio 2, centradas en el origen. 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 variables espaciales [math](x,y)[/math] , y el tiempo [math]t[/math], y los desplazamientos [math]\vec u(\rho,\theta)[/math].De esta forma, si definimos [math]r_0(\rho,\theta)[/math] el vector de posición de los puntos de la placa antes de la deformación , la posición de cada punto [math](\rho,\theta)[/math] de la placa después de la deformación viene dada por: [math] \vec r (\rho,\theta)= \vec r_{0}(\rho,\theta)+\vec u(\rho,\theta). [/math] Vamos a suponer que la fuerza aplicada sobre la placa ha provocado una desplazamiento de manera que los desplazamientos en un instante [math]t_0[/math] vienen dados por: [math]\vec u(\rho,\theta)=(1-\rho)^2 \vec g_{\theta}. [/math]

1 Mallado del sólido

Nuestro mallado consiste en la representación de los puntos interiores de una placa solida, en forma de corona circular. Tomando como paso de muestreo h=1/10, los comandos en Matlab de la visualizacion de nuestra placa es la siguiente:

Mallado de la 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 de las coordenadas u y v
figure(1)
xx=uu.*cos(vv);        % parametrizacion
yy=uu.*sin(vv);
mesh(xx,yy,0*xx)       % dibuja la malla
axis([-3,3,-3,3])      % selecciona la región
view(2)

2 Temperatura del sólido

La temperatura es un campo escalar representado por [math]T(x,y)=e^{-y}[/math],vamos a tratar de ver como varía la temperatura en la placa que estamos estudiando . Nuestra función depende únicamente de la variable Y, esto quiere decir que la temperatura se mantendrá constante para todo valor de X. Observamos que cuanto mayor es el valor de la Y ,menor es la temperatura , luego los puntos de máxima temperatura son los más alejados del foco( el origen) y con menor valor de la variable Y, como podemos observar en la figura de la derecha.

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
figure(1)
xx=uu.*cos(vv);        % parametrización
yy=uu.*sin(vv);
f=exp(-yy); % Campo escalar
surf(xx,yy,f)          % visualización 
axis([-3,3,-3,3])      % región de visualización
view(2)

3 Cálculo y representación del [math]∇T [/math]

Para estudiar como varía la temperatura al movernos por el plano calculamos el gradiente de T. [math]∇T= -e^{-y} j [/math] En la figura de la derecha podemos observar la representación del campo vectorial [math]∇T [/math] y que este es perpendicular a las curvas de nivel del campo escalar T.

g=-exp(-yy); 
subplot(1,2,1);        
hold on
contour(xx,yy,f,20)
surf(xx,yy,f);
quiver3(xx,yy,f,0*xx,g,0*f);
hold off
subplot(1,2,2)
hold on
quiver(xx,yy,0*xx,Gradiente); colorbar;
contour(xx,yy,f,20)
hold off

4 Campo de desplazamientos

Supongamos la aplicación de una fuerza que provoca una vibración sobre la placa, de manera que, en un t aleatorio los desplazamientos vienen dados por: [math]\vec u(\rho,\theta)=(1-\rho)^2 \vec g_{\theta}. [/math] Para distinguir el campo de vectores en los puntos del mallado del sólido utilizaremos el código matlab siguiente:

u=1:h:2;               % definir intervalo rho [1,2]

     v=0:h:2*pi+h;          % definir intervalo theta [0,2*pi]

     [uu,vv]=meshgrid(u,v); % definir mallado

     xx=uu.*cos(vv);        % parametrizacion

     yy=uu.*sin(vv);

     ux=(xx.*(1-uu).^2);

     uy=(yy.*(1-uu).^2);

     subplot(1,2,1);

     mesh(xx,yy,0.*xx)

     view(2)

     axis([-3,3,-3,3])  

     hold on

     quiver(xx,yy,ux,uy)

     title('campo con vectores')

     subplot(1,2,2)

     plot3(ux,uy,0.*xx)

     title('campo de deformaciones')

5 Sólido antes y despues de desplazamiento

El desplazamiento realizado por la placa viene determinado por una serie de movimientos transversales en direcciones distintas. En las imagenes adjuntas calculadas mediante Matlab, se aprecia que el desplazamiento es practicamente nulo.

h=0.1;                              
ro=1:h:2;                           % Valor de ro según el muestreo del intervalo
teta=0:h:2*pi+h;                    % Valor de teta según el muestreo del intervalo
[roro,tetateta]=meshgrid(ro,teta);  %Matrices de la placa en coordenadas polares
figure(1)
xx=roro.*cos(tetateta);             % parametrización
yy=roro.*sin(tetateta);             
subplot(1,3,1)                      % figura solido momento inicial
mesh(xx,yy,0*xx)                    % mallado
axis([-4,4,-4,4]);                  % ejes 
view(2);                            %visualización
desx=xx+((1-roro).^2).*(-yy);       % movimiento parametrizado 
desy=yy+((1-roro).^2).*(xx);        
subplot(1,3,2)                      % figura solido en movimiento
mesh(desx,desy,yy.*0)               % mallado
axis([-4,4,-4,4]);                  
view(2);                            
subplot(1,3,3);                     % figura solido en movimiento y campo vectorial 
hold on                             
mesh(desx,desy,yy.*0)                
axis([-4,4,-4,4]);                  % ejes
view(2);                            % 
quiver(xx,yy,fx,fy)                 % campo vectorial del gradiente de temperaturas
axis([-4,4,-4,4])                   %  ejes 
view(2)                             
hold off

6 Divergencia de u

La Divergencia mide la diferencia de flujo de un campo vectorial sobre una superficie. Como vemos en la gráfica la divergencia está calculada en todos los puntos de la placa y es siempre mayor que cero, lo que nos indica que no tiene pérdidas de flujo, sino que la diferencia es cero por tanto la placa solamente tendrá fuentes.

h=(0.1);                 
ro=1:h:2;               
teta=0:h:2*pi+h;            
[roro,tetateta]=meshgrid(ro,teta);
figure(1)
xx=roro.*cos(tetateta);        
yy=roro.*sin(tetateta)
f=(1./roro)-4+3.*roro; % El campo divergencia
surf(xx,yy,f) % Dibujmos la malla
axis([-2,2,-2,2])    
colorbar
title('Divergencia del campo del desplazamiento')
view(2)

7 Rotacional de u

A traves del programa de Matlab, que podemos observar a continuacion, hemos calculado el rotacional del campo desplazamiento u en los puntos del solido y los hemos representado, observandose en la grafica cuales son los puntos que sufren mayor rotacional.

h=(0.1);                  
ro=1:h:2;              
teta=0:h:2*pi+h;            
[roro,tetateta]=meshgrid(ro,teta); 
figure(1)
xx=roro.*cos(tetateta);        
yy=roro.*sin(tetateta);
f=(2+4.*roro.^2-6.*roro); 
surf(xx,yy,f) 
axis([-2,2,-2,2])      
colorbar
title('Rotacional del campo del desplazamiento')
view(2)

8 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 = \frac{1}{\sqrt{g}}\cdot \ \frac{\partial(\sqrt{g}u^{i})}{\partial x^{i}}=\frac{1}{\rho}\cdot (\frac{\partial \rho - 2\rho ^{2}+ \rho ^{3}}{\partial \rho }+\frac{\partial 0}{\partial \theta })=\frac{1}{\rho }-4+3\rho [/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.= \left(\begin{array}{ccc}-2+2\rho& 0 \\0& \ 0\end{array}\right)+ (1-2\rho +\rho ^{2}).\left(\begin{array}{ccc}0& 0 \\0& \ 1\rho\end{array}\right) [/math]

Podemos calcular el tensor gradiente como 1-contravariante 1-covariante:

[math]\nabla\vec u=u^i_{.j}=\left(\begin{array}{ccc}-2+2\rho& 0 \\0& \ 1/\rho -2 +\rho\end{array}\right) [/math]

Las coordenadas 1-covariantes 1-contravariante coinciden por lo que la transposicion seá directa y el tensor quedara como [math](\epsilon^i_{.j})=(\nabla \vec u + \nabla \vec u^t)/2[/math] , que corrresponde con la matriz mostrada anteriormente

Sustituyendo ahora los valores en la fórmula del tensor de tensiones,y tomando \(\lambda=\mu=1\), siendo \(\lambda\) y \(\mu\) los conocidos como coeficientes de Lamé, obtenemos:

[math]\left(\begin{array}{ccc}1/\rho-8+7\rho& 0 \\0& \ 3/\rho -8 +5\rho\end{array}\right)[/math]

Las tensiones normales de ambos vectores de la base nos quedaran como:

Las tensiones normales en la dirección \(\vec g_{\rho}\) serían:\[\vec g_{\rho} \cdot \sigma \cdot \vec g_\rho =1/\rho-8+7\rho\]

Las tensiones normales en la dirección \(\vec g_\theta/\rho\) serían:\[\vec g_\theta/\rho \cdot \sigma \cdot \vec g_\theta/\rho = 3/\rho -8 +5\rho\]

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=(1./uu)-8+7.*uu; % Campo escalar
 
    surf(xx,yy,f)              % Dibujar el mallado
 
    axis([-2,2,-2,2])          % Región seleccionada
 
    axis equal
 
    colorbar
 
    title('Tension normal en dirección g sub rho')
 
    view(2)                  
 
    subplot(1,2,2)
 
    f=(3./uu)-8+5.*uu;  % Campo escalar
 
    surf(xx,yy,f)              % Dibujar el mallado
 
    axis([-2,2,-2,2])          % Región del dibujo
 
    axis equal
 
    colorbar
 
    title('Tension normal en dirección g sub theta')
 
    view(2)

Vemos como la tensiones serán mayores a medida que nos alejamos del centro de masas de nuestra placa y simétricas respecto de este ya que dependerá de la distancia de cada uno de los puntos al origen. Además se aprecia que en los bordes interiores de las placas dicha tensión será nula.

8.1 Tensiones tangenciales

Al realizar el calculo de las tensiones tangenciales de nuestra placa para las direcciones [math]\vec g_\rho[/math] y a [math]\vec g_\theta/\rho[/math]:: [math]|\sigma \cdot \vec g_\rho-(\vec g_\rho \cdot \sigma \cdot \vec g_\rho) \vec g_\rho| = 0 [/math][math]|\sigma \cdot \vec g_\theta/\rho-(\vec g_\theta/\rho \cdot \sigma \cdot \vec g_\theta/\rho) \vec g_\theta/\rho|= 0 [/math]

Vemos que ambas tensiones son nulas, algebraicamente hablando esto es debido a que nuestro campo de desplazamientos depende únicamente del valor de la distancia al origen. Esto implica que al estar definido dicho tensor por la divergencia y el tensor unitario, será nula la tensión tangencial siempre que el campo no dependa del angulo, ya que sera igual en todas las direcciones.

9 Tensión de Von Mises

Se define por la fórmula: [math] \sigma _{VM}=\sqrt{\frac{(\sigma 1-\sigma 2)^{2}+(\sigma 2-\sigma 3)^{2}+(\sigma 3-\sigma 1)^{2}}{2}}[/math] Esta tensión viene dada proporcionalmente a la Energía de Distorsión. En ingeniería se utiliza como un indicador de tensiones a las cuales el material falla (teoría de fallo), este cálculo se realiza sobre todo a los materiales dúctiles. En nuestro caso el campo de tensiones de Von Misses viene dado por la fórmula [math] \sigma _{vm}=\frac{1}{\sqrt{2}}\cdot \left ( \frac{-2}{\rho }+2\cdot \rho \right ) [/math].

La representaremos mediante el siguiente programa en MATLAB:

h=(0.1);                 
ro=1:h:2;               
teta=0:h:2*pi+h;            
[roro,tetateta]=meshgrid(ro,teta);
figure(1)
xx=roro.*cos(tetateta);        
yy=roro.*sin(tetateta)
f=((-2./roro)-4+2.*roro)/sqrt(2); % El campo de tensiones
surf(xx,yy,f) % Dibujamos la malla
axis([-2,2,-2,2])    
colorbar
title('Tensiones de Von Misses')
view(2)

10 Masa total de la placa

Si la densidad de la placa es [math] d(x,y,z)=xye^{-1/x^{2}} [/math] y el área del anillo es [math] \pi.R^{2}-\pi.r^{2} [/math] ,siendo R y r los radios exteriores e interiores respectivamente, la masa total de la placa será el producto del área y la densidad. Al tener la funcíon de la densidad, necesitamos integrarla para calcularla en todo el dominio. Para una mejor resolución, efectuaremos el cambio a coordenadas polares:

[math] x=\rho \cdot \cos \theta ; y=\rho \cdot \sin \theta [/math]

Tendremos entonces que resolver la integral: [math]\int_{0 }^{2\pi} \int_{1 }^{2}\rho ^{3}\cos \theta . \sin \theta .\log_{10}(\rho .\cos \theta +2)d\rho d\theta [/math] En matlab sería:

h=1/500;                            % Número de particiones del trapecio
u=1:h:2;                           
v=0:h:2*pi+h;                     
[uu,vv]=meshgrid(ro,teta);   % Malla de la placa en polares
xx=uu.*cos(vv);              % Referencia en cartesianas
yy=uu.*sin(vv);         
d=xx.*yy.*log(xx+2);                 % Función de densidad de la placa
a=h^2*d;                             % Masa en cada punto
mas=sum(sum(a))                      % Suma de todas las masas puntuales)

La masa de la placa es [math] 2,097\cdot 10^{-6} [/math]