Deformaciones elásticas de un anillo circular (Grupo A2)

Considerando una placa en forma de anillo circular realizaremos una serie de visualizaciones en 2 dimensiones de diversos campos actuando sobre ella. En primer lugar estudiaremos el efecto de la temperatura dada por [math]T(\rho,\theta,t)[/math] para después plantear unas vibraciones [math]\vec u(\rho,\theta,t)[/math] de las que analizaremos cómo actúan y las tensiones que producen. Utilizaremos, además de las cartesianas, coordenadas polares ρ y θ y definiremos [math]r_0(\rho,\theta)[/math] como 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]

1 Visualización del sólido (mallado)

El anillo circular está comprendido entre las circunferencias centradas en el origen de radios 1 y 2. Para definirlo creamos una malla delimitada por el dominio de las coordenadas polares utilizadas:ρ y θ (ρ є [1,2], θ є [0,2Π)), y planteamos las expresiones de la parametrización.

Sólido en forma de anillo circular

1.1 Código MATLAB

h=0.1;     % Paso de muestreo
%Intervalos en los que se definen las polares
r=1:h:2;                            
t=0:h:2*pi+h;          
[rr,tt]=meshgrid(r,t);   %Creación del mallado        
%Parametrización
xx=rr.*cos(tt);                 
yy=rr.*sin(tt);
mesh(xx,yy,0*xx)           
axis([-3,3,-3,3])

2 Efecto térmico

2.1 Distribución de temperaturas en el anillo

En el origen de coordenadas (centro de las circunferencias que forman la placa) se sitúa un foco de calor. La temperatura en cada punto del sólido debido a este aumento térmico se puede hallar mediante la expresión en coordenadas polares:[math]T(\rho,\theta)=-\log(\rho+0.1)[/math] En el gráfico se muestra la distribución térmica en el anillo. Se trata de una distribución con centro de simetría el foco de calor (origen de coordenadas) e independiente del ángulo θ. De esta forma se observan circunferencias concéntricas formadas por puntos que se encuentran a la misma temperatura (curvas de nivel).

Distribución de temperaturas del anillo.

2.1.1 Código MATLAB

h=0.1;                               
r=1:h:2;                            
t=0:h:2*pi+h;          
[rr,tt]=meshgrid(r,t);     
xx=rr.*cos(tt);                 
yy=rr.*sin(tt);
mesh(xx,yy,0*xx
axis([-3,3,-3,3])    
f=-log(rr+0.1);   % Expresión del campo escalar térmico
surf(xx,yy,f)          
axis([-3,3,-3,3])

2.2 Gradiente térmico

Al calcular el gradiente del campo escalar T, [math]\nabla T[/math] resulta un campo vectorial de expresión:[math] \nabla T=\frac{\cos(\theta)}{0.1-ρ}\vec i + \frac{\sin(\theta)}{0.1-ρ}\vec j [/math] La representación del campo vectorial gradiente indica, en cada punto del sólido, la dirección en la cual la temperatura aumenta más rápido, y su módulo nos dirá cuan rápido aumenta la temperatura en esa dirección.Es decir, obtenemos la dirección en la cual el campo T varía más rápidamente. Estas direcciones coincidirán con el vector que pase de forma más rápida de una curva de nivel a la siguiente, por tanto, el gradiente será perpendicular a dichas curvas de nivel, como podemos observar en el gráfico. La dirección de los vectores del campo es aquella que une el punto en el que están situados con el foco térmico, con sentido hacia éste, pues es el punto con temperatura máxima y cuanto más cerca se encuentren de él más alta será la temperatura. Resulta un campo vectorial cuyos vectores están orientados hacia el origen de coordenadas (centro de simetría del campo), punto en el que dicho campo es nulo.

Gradiente y curvas de nivel del campo de temperaturas.

2.2.1 Código MATLAB

h=0.1;
r=1:h:2;
t=0:h:2*pi+h;
[rr,tt]=meshgrid(r,t);
xx=rr.*cos(tt);
yy=rr.*sin(tt);
%Relación polares-cartesianas
u=sqrt(xx.^2+yy.^2);
v=atan(yy./xx); 
f=-log(0.1+u);
fx=-1./(u.*(0.1+u)).*xx;     % Derivada parcial en x
fy=-1./(u.*(0.1+u)).*yy;     % Derivada parcial en y
hold on
quiver(xx,yy,fx,fy)
axis([-3,3,-3,3]) 
contour(xx,yy,f,10)    %Curvas de nivel
axis([-3,3,-3,3])
hold off

3 Efecto de la vibración

A continuación vamos a aplicar sobre la placa una fuerza que origina una vibración de tal forma que los desplazamientos en un tiempo t0 de los puntos de la placa vienen dados por: [math] \vec u(\rho,\theta)=\frac{\sin(\pi \theta/2)}{30\rho}\vec g_{\rho}. [/math]

3.1 Campo vectorial de vibraciones

El campo vectorial tiene el siguiente aspecto sobre los puntos del mallado de la placa. Cada vector representa la fuerza aplicada en cada punto del sólido. De esta forma su dirección nos informa de la dirección que tendrán los desplazamientos de los puntos de la placa.

Campo vectorial u.

3.1.1 Código MATLAB

h=0.1;
r=1:h:2;
t=0:h:2*pi+h;
[rr,tt]=meshgrid(r,t);
xx=rr.*cos(tt);
yy=rr.*sin(tt);
u=sqrt(xx.^2+yy.^2);
v=atan(yy./xx); 
ux=sin(pi.*v./2)./(30.*u).*(xx./u);  % Componente x de u
uy=sin(pi.*v./2)./(30.*u).*(yy./u);   % Componente y de u
quiver(xx,yy,ux,uy)
axis([-3,3,-3,3])

3.2 Desplazamientos provocados

A continuación podemos obtener la posición final de los puntos de la placa al actuar sobre ella las vibraciones anteriores. Al comparar ambos estados, el inicial y el final, observamos que las partículas del sólido efectivamente se mueven en la dirección de los vectores del campo u, como antes habíamos supuesto. Estos desplazamientos se aprecian sobre todo en los puntos de θ=Π/2 y θ=3Π/2, donde existen una especie de discontinuidades en el contorno del sólido (ρ=1 y ρ=2).
Como estos desplazamientos son pequeños, en el segundo gráfico se muestra en detalle la discontinuidad producida en el contorno interior.

Comparativa entre la posición inicial de los puntos del anillo y la posición final después de haberle aplicado el campo vectorial u.

Comparativa en detalle.

3.2.1 Código MATLAB

h=0.1;                  
r=1:h:2;               
t=0:h:2*pi+h;            
[rr,tt]=meshgrid(r,t); 
xx=rr.*cos(tt);       
yy=rr.*sin(tt);
u=sqrt(xx.^2+yy.^2);
v=atan(yy./xx); 
%Desplazamientos originados por la vibración
ux=sin(pi.*v./2)./(30.*u).*(xx./u);
uy=sin(pi.*v./2)./(30.*u).*(yy./u);
%Posición final de los puntos de la placa
dx=xx+ux;
dy=yy+uy;
subplot(1,2,1)
mesh(xx,yy,0.*xx)       
axis([-3,3,-3,3])      
subplot(1,2,2)
mesh(dx,dy,0.*dx)
axis([-3,3,-3,3])

3.3 Divergencia

Para analizar el cambio de volumen local debido al desplazamiento anterior utilizamos el operador diferencial de la divergencia [math] \nabla \cdot \vec u [/math]. Su resultado nulo indica que este cambio de volumen local también es nulo y por tanto el sólido lo consideraremos incompresible. Debido a esto no se pueden destacar puntos notables.

3.4 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. Llegados a este punto, es obligado pasar a la visualización en 3 dimensiones por la propia definición de rotacional. Resulta el campo vectorial:[math] \nabla×\vec u=-\frac{\pi}{2}\frac{\cos(\pi \theta/2)}{30\rho^2}\vec g_{z}. [/math] Este campo vectorial está integrado por vectores en la dirección del eje z con sentido negativo. Pero donde de verdad se puede apreciar el giro que induce el campo u es en la visualización del módulo de los vectores del campo rotacional: en el segundo gráfico hay 2 figuras, una en 3 dimensiones y otra en la que se observa el plano XY en verdadera magnitud. De esta forma vemos claramente qué puntos del sólido experimentan una mayor tendencia al giro. Al observar este último con el tercer gráfico, que superpone el campo vectorial u al modulo del rotacional llegamos a la conclusión de que el sólido es menos proclive a girar en los puntos en los que se produce una inversión del sentido del campo vibratorio (θ=Π/2 y θ=3Π/2). Por el contrario, su tendencia al giro aumenta al acercarnos a θ=0 y θ=Π y, en menor medida, también aumenta al estar más próximos al origen.

Campo vectorial rotacional.

Módulo del rotacional en el sólido.

Superposición del campo vectorial u y el escalar del módulo del rotacional.

3.4.1 Código MATLAB

h=0.1;                  
r=1:h:2;               
t=0:h:2*pi+h;            
[rr,tt]=meshgrid(r,t); 
xx=rr.*cos(tt);        
yy=rr.*sin(tt);      
u=sqrt(xx.^2+yy.^2);
v=atan(yy./xx); 
ux=sin(pi.*v./2)./(30.*u).*(xx./u); 
uy=sin(pi.*v./2)./(30.*u).*(yy./u);  
%Función que resulta al multiplicar vectorialmente el operador nabla y el campo u
r=inline('(pi*(cos(pi/2*atan(yy./xx))))./(60*sqrt(xx.^2+yy.^2))','xx','yy');
Rot=r(xx,yy);
figure(1)
quiver3(xx,yy,0.*xx,0.*xx,0.*yy,-Rot)
figure(2)
subplot(1,2,1)
surf(xx,yy,Rot)
subplot(1,2,2)
mesh(xx,yy,Rot)
figure(3)
hold on
mesh(xx,yy,Rot)
quiver(xx,yy,ux,uy)
hold off

4 Deformaciones elásticas

Si ahora consideramos que estamos en un medio elástico lineal (la relación tensión-deformación es una relación lineal), isótropo (sus propiedades no dependen de la posición) y homogéneo (propiedades iguales en todos los puntos); podemos definir el tensor de tensiones a partir de los desplazamientos.

La fórmula que define [math]\sigma_{ij}[/math] es la siguiente: [math] \sigma_{ij}=\lambda \nabla \cdot \vec u \delta_{ij} + 2\mu \epsilon_{ij} [/math]

Donde [math] \lambda [/math] y [math] \mu [/math] son conocidos como los coeficientes de Lamé que dependen de las propiedades elásticas de cada material y [math] \epsilon_{ij} [/math] es la parte simétrica del tensor gradiente del campo vectorial [math] \vec u [/math].

Tomando [math] \lambda =1[/math] y [math] \mu =1[/math], el tensor de deformaciones en forma matricial resulta:

[math](\sigma_{ij})= \begin{bmatrix}
    \frac{-\sin(\pi \theta/2)}{15\rho^2} & \frac{\pi\cdot\cos(\pi \theta/2)}{60\rho} & 0 \\
    \frac{\pi\cdot\cos(\pi \theta/2)}{60\rho} & \frac{\sin(\pi \theta/2)}{15\rho^2} & 0 \\
    0 & 0 & 0
  \end{bmatrix}
[/math]

4.1 Tensiones normales

4.1.1 Tensiones normales en la dirección de [math]\vec g_\rho[/math]

Las tensiones normales en la dirección radial se obtienen mediante [math]\vec g_{\rho} \cdot \sigma \cdot \vec g_\rho[/math]. Resultan máximas en el punto [math](\rho,\theta)= (1,\pi)[/math] como se muestra en el gráfico:

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

Representación del tensor normal en la dirección de [math]\vec g_\rho[/math]

4.1.1.1 Código MATLAB

h=0.1;                         
r=1:h:2;               
t=0:h:2*pi+h;            
[rr,tt]=meshgrid(r,t);     
xx=rr.*cos(tt);     %Relación polares-cartesianas  
yy=rr.*sin(tt);
f=((-sin(pi*tt./2))./(15*(rr.^2))); % Campo escalar
surf(xx,yy,f)              
axis([-3,3,-3,3])          
axis equal

4.1.2 Tensiones normales en la dirección de [math]\vec g_\theta/\rho[/math]

Al calcular las tensiones en la dirección de [math]\vec g_\theta/\rho[/math], es decir, hacer la operación [math]\vec g_\theta/\rho \cdot \sigma \cdot \vec g_\theta/\rho[/math], vemos que resultan ----------- como se observa en el gráfico:

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

Representación del tensor normal en la dirección de [math]\vec g_\theta/\rho[/math]

4.1.2.1 Código MATLAB

h=0.1;                         
r=1:h:2;               
t=0:h:2*pi+h;            
[rr,tt]=meshgrid(r,t);     
figure(1)
xx=rr.*cos(tt);     
yy=rr.*sin(tt);
f=((sin(pi*tt./2))./(15*(rr.^2))); % Campo escalar
surf(xx,yy,f)              
axis([-3,3,-3,3])         
axis equal

4.2 Tensiones tangenciales

4.2.1 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·\cos(\pi \theta/2)}{60\rho}\vec g_\theta| = [/math]

h=0.1;                         
    r=1:h:2;               
    t=0:h:2*pi+h;            
    [rr,tt]=meshgrid(r,t);     
    xx=rr.*cos(tt);     %Relación polares-cartesianas  
    yy=rr.*sin(tt);

    f=(pi*(cos(pi*tt./2)))./(60.*rr); %Campo escalar

    surf(xx,yy,f)
    axis([-3,3,-3,3])          
    axis equal 
    colorbar
    view(2)

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

h=0.1;                         
    r=1:h:2;               
    t=0:h:2*pi+h;            
    [rr,tt]=meshgrid(r,t);     
    xx=rr.*cos(tt);     %Relación polares-cartesianas  
    yy=rr.*sin(tt);
    f=(pi*(cos(pi*vv./2)))./(60.*rr); %Campo escalar

    surf(xx,yy,f) 
    axis([-3,3,-3,3])          
    axis equal 
    colorbar
    view(2)