Deformación de un semianillo plano Grupo 4-A

De MateWiki
Saltar a: navegación, buscar


Trabajo realizado por estudiantes
Título Deformación de un semianillo plano Grupo 4-A
Asignatura
Curso 2021/22
Autores Jaime Guerrero Suárez

Sergio Míguez González Pedro Michelini

Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura.

1 Enunciado

Consideramos una placa plana que ocupa la mitad de un anillo circular centrado en el origen y comprendido entre los radios 1 y 2, en el plano y ≥ 0. En ella vamos a suponer que tenemos definidas dos cantidades físicas: la temperatura \(T(x,y)\), que viene dada por [math]T(x, y) = log((x − 3)^2 + 2)[/math], y los desplazamientos \(\vec u(x,y)\) producidos por la acción de una fuerza determinada. De esta forma, si definimos \(\vec r\)0(x,y) el vector de posición de los puntos de la placa en reposo, la posición de cada punto \((x,y)\) de la placa después de la deformación viene dada por \(\vec r\)(x, y) = ~\(\vec r\)0(x, y) + ~\(\vec u\)(x, y). Vamos a suponer que la fuerza aplicada sobre la placa ha provocado un desplazamiento de los puntos de la misma dado por el vector de desplazamientos [math]\vec u(ρ, θ)=\frac{ρ − 1}{5} sin(θ)~[/math]\(\vec e\)θ.

2 Mallado

Comenzamos la representación del sólido definiendo un mallado siguiendo las directrices del enunciado: "Tomar los ejes (comando axis) en el rectángulo (x, y) ∈ [-3; 3] × [−1; 3] y como paso de muestreo h = 1/10 para las variables x e y. Al tratarse de una figura plana en 2D, la componente z es nula.

Representación del mallado
h= 0.1;                           %Paso de muestreo
r= 1:h:2;
tt= 0:h:pi;
[RR,TT]= meshgrid(r,tt);       %Mallado
x=RR.*cos(TT);
y=RR.*sin(TT);                                                        
mesh(x,y,0.*x)                  %Visualización de la placa
view(2)                           
axis ([-3,3,-1,3])
xlabel('Eje X')
ylabel('Eje Y')


3 Representación de la temperatura

Una vez tenemos definido el mallado (que representa al propio sólido), pasamos a la representación de la temperatura a partir de la función [math]T(x, y) = log((x − 3)^2 + 2)[/math]. Para llevarla a cabo, representamos el campo escalar. Además, hemos definido una serie de curvas de nivel para poder apreciar mejor el cambio de temperatura a lo largo de la placa.

Temperaturas
h= 0.1;                           %Paso de muestreo
r= 1:h:2;
tt= 0:h:pi;
[RR,TT]= meshgrid(r,tt);       %Mallado
x=RR.*cos(TT);
y=RR.*sin(TT);  
T=log10((x-3).^2 + 2);
hold on
subplot(1,2,1)                %Dividimos la pantalla en dos
surf(x,y,T)                   %Representamos el campo escalar de temperaturas
view(2)
axis ([-3,3,-1,3])
xlabel('Eje X')
ylabel('Eje Y')
colorbar                      %Mostramos la escala
subplot(1,2,2)                %Escribimos en la segunda pantalla
contour(x,y,T,30)
title('Curvas de Nivel','Fontsize',16)      %Líneas de nivel
colorbar                      %Mostramos las escala
axis ([-3,3,-1,3])
xlabel('Eje X')
ylabel('Eje Y')
hold off

Del gráfico se deduce que la temperatura aumenta en la dirección negativa de x. Por lo tanto, el valor máximo se alcanza cuando x=-2.

4 Gradiente de Temperatura

El siguiente paso después de observar el campo de la temperatura es calcular el gradiente de esta misma función y dibujarlo como campo vectorial. El cálculo del gradiente consiste en:

Fórmula del gradiente de una función

En nuestro caso, el gradiente calculado manualmente nos da [math]\nabla[/math]T=[math]\frac{2(x-3)}{(x-3)^2+2}[/math]\(\vec i\). Como sabemos, por definición el gradiente es siempre perpendicular a la función original. Una vez que lo tenemos ya calculado, podemos observarlo y confirmar que, efectivamente, se trata de un campo vectorial perpendicular a las curvas de nivel formadas por la temperatura.

Gradiente de la temperatura comparado con sus curvas de nivel
h= 0.1;                           %Paso de muestreo
r= 1:h:2;
tt= 0:h:pi;
[RR,TT]= meshgrid(r,tt);       %Mallado
x=RR.*cos(TT);
y=RR.*sin(TT);                                                                                               
T=log10((x-3).^2 + 2); 
hold on       %Campo escalar de temperatura            
contour(x,y,T,60)             %Líneas de nivel
colorbar                      %Mostramos las escala
axis ([-3,3,-1,3])
ti=inline('(2.*(x-3))./(log(10).*((x-3).^2)+2)','x','y'); %Derivada parcial con respecto a x
tj=inline('0.*x','x','y');                                %Derivada parcial con respecto a y   
Tx=ti(x,y);
Ty=tj(x,y);
quiver(x,y,Tx,Ty);
axis ([-3,3,-1,3])
title ('Gradiente de la Temperatura','Fontsize',18);
set(h,'maxheadsize',0.33)
hold off


5 Representación del campo vectorial de deformaciones

Ahora nos centraremos en el campo \(\vec u\) de deformaciones, definido en el enunciado en coordenadas cilíndricas (factor a tener en cuenta). No necesitamos calcularlo, pero sí convertir la dirección de eθ en las direcciones \(\vec i\) y \(\vec j\) para poder representarlo sin problema en la misma cuadrícula. Aplicamos la equivalencia eθ=[math]\frac{-x_2}{ρ}[/math]\(\vec i\)+[math]\frac{x_1}{ρ}[/math]\(\vec j\)=-sin(θ)\(\vec i\)+cos(θ)\(\vec j\). El nuevo campo \(\vec u\) nos queda de la siguiente forma: \(\vec u\)=-[math]\frac{ρ-1}{5}[/math]sin2(θ)\(\vec i\)+[math]\frac{ρ-1}{5}[/math]cos(θ)sin(θ)\(\vec j\)

Campo de deformaciones sobre los diversos puntos del mallado.
h= 0.1;                     %Paso de muestreo                              
rr= 1:h:2;                  %Usamos coordenadas polares.
tt= 0:h:pi;
[RR,TT]= meshgrid(rr,tt);   %Mallado
x=RR.*cos(TT);
y=RR.*sin(TT);  
hold on  
mesh(x,y,0.*x)                  %Visualización de la placa
view(2)                 
axis ([-3,3,-1,3])
a=-((RR-1)./5).*((sin(TT)).^2);
b=((RR-1)./5).*(cos(TT).*sin(TT));
w=quiver(x,y,a,b);
axis ([-3,3,-1,3])
title('Campo de vectores','Fontsize',18);
set(w,'maxheadsize',0.33)


6 Desplazamiento del mallado producido por el campo vectorial

Ya tenemos el campo vectorial del desplazamiento, por lo que vamos a poder observar la deformación en la placa antes y después de la acción del campo \(\vec u\):

La placa antes y después de la acción del campo de desplazamientos
Comparación antes (verde) y después del desplazamiento (rojo)
h= 0.1;                          %Paso de muestreo                                 
rr= 1:h:2;                       %Usamos coordenadas polares
tt= 0:h:pi;
[RR,TT]= meshgrid(rr,tt);        %Mallado
x=RR.*cos(TT);
y=RR.*sin(TT);
                                 %Sólido antes de los desplazamientos
subplot(1,2,1)                  
i=mesh(x,y,0*x);
view(2)
set(i,'EdgeColor','g');
axis ([-3,3,-1,3])
xlabel('Eje X')
ylabel('Eje Y')
title('Placa no desplazada','Fontsize',16); 
                                 %Sólido después de los desplazamientos
subplot(1,2,2)                     
A=-((RR-1)./5).*((sin(TT)).^2);
B=((RR-1)./5).*(cos(TT).*sin(TT)); 
X=x+A;
Y=y+B;
j=mesh(X,Y,0*X);
view(2)
set(j,'EdgeColor','r');
axis ([-3,3,-1,3])
xlabel('Eje X')
ylabel('Eje Y')
title('Placa desplazada','Fontsize',16);
%Comparación de ambas representaciones. 
hold on
figure                                 
j=mesh(X,Y,0*X);
axis ([-3,3,-1,3])
view(2)
set(j,'EdgeColor','r');
title('Desplazamiento de la placa','Fontsize',16);
hold on
i = mesh(x,y,0*x);
set(i,'EdgeColor','g');
hold off


Como podemos apreciar al comparar estas dos gráficas, el desplazamiento producido por el campo \(\vec u\) ha creado una ligera torsión radial en la placa y una elongación del extremo derecho.

7 Divergencia del campo vectorial

A continuación, calculamos la divergencia del campo vectorial. La divergencia es una medida del cambio de volumen local debido al desplazamiento, es decir, que una parte del sólido aumentará de volumen y otra lo perderá (dilatación o contracción). Al calcularla nos queda: DivU.png

Divergencia
h= 0.1;                        %Paso de muestreo
                               %Usamos coordenadas polares
rr= 1:h:2;
tt= 0:h:pi;
[RR,TT]= meshgrid(rr,tt);      %Mallado
x=RR.*cos(TT);
y=RR.*sin(TT);
DIVu=((1)./RR).*(cos(TT).*((RR-1)./5)); %Divergencia de u
surf(x,y,DIVu);
view(2);
axis ([-3,3,-1,3])
colorbar;
title('Divergencia','Fontsize',16);
xlabel('Eje X')
ylabel('Eje Y')


Del gráfico se puede deducir que la divergencia toma valores máximos en el punto (ρ, θ) = (2,0) y valores mínimos en el punto (ρ, θ) = (2,π). Por otro lado, vemos que la divergencia es nula cuando ρ=1 (cualquiera sea el θ) y cuando θ=π/2 (cualquiera sea el ρ). Como vemos en la gráfica, la divergencia es positiva cuando θ se encuentra entre 0 y π/2 (se dilata); es negativa con θ entra π/2 y 2π (se contrae) y al ser nula entorno a π/2, esa región del sólido se mantiene constante en volumen.


8 Rotacional del campo vectorial

Continuamos calculando el rotacional, con resultado: Formularottt.png en \(\vec e\)z. Calculamos su módulo y graficamos:

Rotacional
h=0.1;                                       %Paso de muestreo                                
rr=1:h:2;                                    %Usamos coordenadas polares
tt=0:h:pi;
[RR,TT]= meshgrid(rr,tt);                    %Mallado
x=RR.*cos(TT);
y=RR.*sin(TT);
ROTu=abs((sin(TT)./(5.*RR)).*(2-((1)./RR))); %Modulo del rotacional de u
                                             %Rotacional en 2-D
subplot(1,2,1);
surf (x,y,ROTu);
axis equal
view(2);
colorbar;
title('Rotacional en 2-D', 'Fontsize',16);
xlabel('Eje X')
ylabel('Eje Y')
                                            %Rotacional en 3-D
subplot(1,2,2);
surf(x,y,ROTu);               
colorbar
title('Rotacional en 3-D', 'Fontsize',16);
xlabel('Eje X')
ylabel('Eje Y')
zlabel('Eje Z')


Observamos en el gráfico que los puntos que sufren un mayor rotacional son los que se aproximan al punto (ρ, θ) = (1,π/2). Recientemente hemos aprendido a interpretarlo mediante el teorema de Stokes, y viendo la gráfica en 3D, sabemos que la velocidad en θ=0 y en θ=π es nula en cuanto a rotación, y es máxima en el punto mencionado antes, lo que genera una tensión plástica que lo deforma.

9 Tensiones normales a los ejes coordenados

Una vez definida la parte simétrica del tensor gradiente de u conocido como tensor de deformaciones. siendo este ε(u). Entonces conoceremos el tensor de tensiones definido como σ. Sabemos que σ se puede obtener a partir de la fórmula:

centro

Tomamos los coeficientes de Lamé (λ y μ) ambos como iguales a 1, y el tensor de deformaciones como la suma del gradiente de \(\vec u\) con su matriz traspuesta divididos por 2. EL resultado de sigma lo podemos ver incluido en el código del programa.

Después podremos representar las tensiones normales en la dirección que marca el eje i, es decir [math]\vec i [/math] = \(\vec g_{i}\) * σ * \(\vec g_{i}\), las tensiones normales en la dirección que marca el eje j, es decir [math]\vec j [/math] = \(\vec g_{j}\) * σ * \(\vec g_{j}\) y las correspondientes al eje k, es decir [math]\vec k[/math] = \(\vec g_{k}\) * σ * \(\vec g_{k}\). Resultan las siguientes gráficas:

Tensiones normales a los ejes
Vista plano XY
r=1:h:2;
t=0:h:pi;
[rr,tt]=meshgrid(r,t); 
clf
subplot(1,3,1);
xx=rr.*cos(tt);
yy=rr.*sin(tt);
a=(rr-1).*cos(tt)./(5.*rr); %Elementos (1,1) y (3,3) de la matriz sigma 
b=(2-rr).*sin(tt)./5;    %elementos (2,1) y (1,2)
c= 2.* (rr-1).*cos(tt)./5;   %Elemento (2,2) de la matriz Épsilon
d=a+c; %Elemento (2,2) de la matriz Sigma
surf(xx,yy,a)
axis equal
view(2)
colorbar
subplot(1,3,2);   %Elemento (2,2) de la matriz sigma
surf(xx,yy,d)
axis equal
view(2)
colorbar
subplot(1,3,3)
surf(xx,yy,a)
axis equal
view(2)
colorbar
figure
subplot(3,1,1);
surf(xx,yy,a)
colorbar
subplot(3,1,2);   %Elemento (2,2) de la matriz sigma
surf(xx,yy,d)
colorbar
subplot(3,1,3);   %Elemento (3,3) de la matriz sigma
surf(xx,yy,a)
colorbar


10 Tensiones tangenciales respecto al plano ortogonal a \(\vec i\)

La siguiente gráfica representa las tensiones superficiales a la placa. Como se ve en el resultado, estas incrementan cuando θ tiende a π/2.

Tensiones tangenciales.
h=0.1;
r=1:h:2;
t=0:h:pi;
[rr,tt]=meshgrid(r,t);
xx=rr.*cos(tt);
yy=rr.*sin(tt);
a=(2-rr).*sin(tt)./5;
b=(2-rr).*sin(tt)./5;
surf(xx,yy,a)
axis equal
view(2)
colorbar
subplot(1,3,2);   
surf(xx,yy,b)
axis equal
view(2)
colorbar


11 Tensión de Von Mises

Debemos calcular la Tensión de Von Mises. Se trata de una magnitud escalar que se suele usar como indicador para saber cuándo un material inicia un comportamiento plástico. Emplearemos la siguiente fórmula, que toma como valores de σ(1), σ(2) y σ(3) los autovalores de la matriz σ calculada en el apartado 9. En el programa incluimos una calculadora de autovalores. TVONMISES.png

Tensión de Von Mises
h=0.1;
rr=[1:h:2];
tt=[0:h:pi];
[RR,TT]=meshgrid(rr,tt); 
sigma=zeros(3,3);%Creación matriz sigma y de Von Mises              
VM=zeros(32,11); %Tomando cada punto del mallado (dimensiones)
for i=1:32*11
     r=RR(i)';
     t=TT(i)';
     sigma(1,1)=(r-1)./(5.*r).*cos(t);
     sigma(1,2)=(2-r)./5.*sin(t);
     sigma(2,1)=(2-r)./5.*sin(t);
     sigma(2,2)=(2.*r.^2-r-1)./(5.*r).*cos(t);
     sigma(3,3)=(r-1)./(5.*r).*cos(t);
     [v,d]=eig(sigma);                           %obtenemos autovalores
     VM(i)=sqrt(((d(1,1)-d(2,2))^2+(d(2,2)-d(3,3))^2+(d(3,3)-d(1,1))^2)/2);  %Fórmula de Tensión de Von Mises
end
                %Pasamos a cartesianas  
xx=RR.*cos(TT);        
yy=RR.*sin(TT);
surf(xx,yy,VM)
title('Tension de Von Mises')
axis([-3,3,-1,3])
xlabel('Eje X')
ylabel('Eje Y')
colorbar
view(2)
end


Como se observa, el mayor valor se alcanza cuando (ρ, θ) = (2,0) y en (2,π).

12 Campo de fuerzas \(\vec F\) que actúa sobre la placa

El campo de fuerzas que actúa sobre la placa sigue la siguiente fórmula F = −∇ · σ. En la gráfica podemos apreciar claramente cómo actúa sobre el sólido este campo y podríamos intuir qué desplazamiento y deformación sufriría.

Campo de fuerzas que recibe la placa
h= 0.1;                     %Paso de muestreo                              
rr= 1:h:2;                  %Usamos coordenadas polares.
tt= 0:h:pi;
[RR,TT]= meshgrid(rr,tt);   %Mallado
x=RR.*cos(TT);
y=RR.*sin(TT);  
z=0.*x;  
hold on
mesh(x,y,0.*x)                  %Visualización de la placa
view(2)                 
axis ([-3,3,-1,3])             
axis ([-3,3,-1,3])
a=(2.*RR.^2-RR.^3-1)./(5.*RR.^2).*cos(TT);
b=-(2.*RR.^2-2.*RR-1)./(5.*RR).*sin(TT);
c=-1./RR.*cos(TT);
w=quiver3(x,y,z,a,b,c);
axis ([-3,3,-1,3,-1,1])
title('Campo de fuerzas','Fontsize',18);
set(w,'maxheadsize',0.33)
view(2)
hold off

Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Deformación_de_un_semianillo_plano_Grupo_4-A&oldid=48330»