Visualización de campos en elasticidad en una placa plana semicircular

De MateWiki
Saltar a: navegación, buscar
Trabajo realizado por estudiantes
Título Visualización de campos en elasticidad en una placa plana semicircular. Grupo 15, Trabajo 5
Asignatura Teoría de Campos
Curso 2020-21
Autores

Alcolea Herreros, Javier (Grupo C)
España Nata, Camila (Grupo A)
Del Fresno Oca, Mariano (Grupo C)
Vitoria Verdejo, Juan (Grupo C)

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

En este artículo se analizan los efectos que causan distintos campos escalares y vectoriales sobre un sólido. Consideramos la placa plana que ocupa la mitad de un anillo circular centrado en el origen y comprendido entre los radios [math]1[/math] y [math]2[/math], en el plano [math]y≥0[/math].
En la placa hay dos expresiones físicas que dependen de las variables [math]ρ[/math] y [math]θ[/math]:

  • La temperatura \(T(ρ,θ)\).
  • El campo de desplazamientos [math]\vec u(ρ,θ)[/math] producido por la acción de una fuerza aplicada sobre la placa.


1 Mallado de la placa

Para la representación del mallado de la placa hemos parametrizado la superficie usando coordenadas cilíndricas. El dibujo se ha representado usando Matlab. El paso de muestreo elegido ha sido de [math]h=0.1[/math], para las variables [math]x[/math] e [math]y[/math]. La placa ha sido representada en el cuadrado [math][-3,3]×[-1,3] [/math].

Figura 1
% Región de la placa
u=1:0.1:2;
h=round(pi/0.1);
v=linspace(0,pi,h);
 
% Matriz de u y v
[u,v]=meshgrid(u,v);
x=u.*cos(v);
y=u.*sin(v);
z=u.*0;
 
%Mallado
i=mesh(x,y,z);
 
%Región de dibujo
axis([-3,3,-1,3])
view(2)
title('Placa plana');

2 Temperatura de la placa

El siguiente paso será representar a partir del campo escalar de temperaturas : [math]T(x,y)=x^2+(y-1)^2[/math] (como puede verse, en coordenadas cartesianas), la temperatura en cada punto de la placa usando de un sistema de colores, de temperaturas más frías (colores azules) a temperaturas más altas (colores amarillos). Se han hecho dos representaciones, la de la izquierda muestra el campo escalar de temperaturas y la de la derecha representa las curvas de nivel del campo [math]T(x,y)[/math] en la placa.

Figura 2
function temperatura 
%Parametrización u,v
h=10; %paso de muestro
u = linspace(1,2,h); 
v = linspace(0,pi,h); 
[U,V] = meshgrid(u,v); 
x = U.*cos(V);
y = U.*sin(V);
 
%función temperatura y curvas de nivel
T =(x.^2)+(y-1).^2;
subplot(1,2,1)
surf(x,y,T)
title ('Temperatura')
view (2)
axis equal
axis([-3,3,-1,3])
subplot(1,2,2)
contour (x,y,T,60)
title ('Curvas de nivel')
axis equal
axis([-3,3,-1,3])
colorbar %barra de colores del gradiente de Ta
Tmax=max(max(T)) %Temperatura máxima de la placa

3 Gradiente de la placa

En este apartado, una vez calculada la temperatura, calculamos el gradiente, que representa la variación de la temperatura en un movimiento en una dirección determinada. Derivamos respecto de cada variable la función de temperatura [math]T(x,y)=x^2+(y-1)^2[/math], resultando [math]\nabla(T)=2x\vec i + 2(y-1)\vec j[/math]. Además calculamos el campo vectorial que hemos obtenido sobre el campo escalar gráficamente. Los vectores son perpendiculares a las curvas de nivel, como se puede ver en la imagen.

Figura 3
% Región de la placa 
u=1:0.1:2; 
h=round(pi/0.1); 
v=linspace(0,pi,h); 
% Matriz de u y v 
[u,v]=meshgrid(u,v); 
x=u.*cos(v); 
y=u.*sin(v); 
z=u.*0;   
%gradiente de la temperatura
T =x.^2+(y-1).^2;
Tx=(2.*x);
Ty=2.*(y-1);
%dibujo
contour(x,y,T,h);
view (2)
axis equal
axis([-3 3 -1 3])
hold on
quiver(x,y,Tx,Ty);
axis equal
axis([-3 3 -1 3])
view (2)
hold off

4 Deformación de la placa

Hemos considerado el campo de desplazamientos [math]\vec u[/math]. Lo hemos calculado sabiendo que:

  • Los puntos situados en el eje [math]ρ=1[/math] no sufren desplazamiento.
  • [math]\nabla·\vec u(ρ,θ)=(2-\frac{1}{ρ})/5[/math].

Dada la segunda condición nombrada, y sabiendo cual es la definición de divergencia obtenemos la siguiente igualdad:

Figura 4a

Al integrar a ambos lados:

Figura 4b

Y dada la primera condición que dice que los puntos en ρ=1 no sufren desplazamiento:

Figura 4c

Por lo que el campo de desplazamientos se expresará como:

Figura 4d

5 Campo de vectores

A continuación, representamos el campo de vectores (desplazamientos) del apartado anterior. Usamos una base ortonormal.

Figura 5

El campo de vectores es mayor donde son mayores las flechas de la imagen, en nuestro caso en la parte exterior de la figura (en [math]ρ=2[/math]). 

%Discretización de los parámetros de la superficie 
u=1:0.1:2;  
h=round(pi/0.1);  
v=linspace(0,pi,h); 
%Creación del mallado 
[U,V]=meshgrid(u,v); 
%Componentes en las direcciones de i y de j del campo de desplazamientos 
x=U.*cos(V); 
y=U.*sin(V); 
z=U.*0; 
a=(1/5).*(U-1).*cos(V); 
b=(1/5).*(U-1).*sin(V); 
c=b.*0; 
figure 
%Dibujo del mallado 
i=mesh(x,y,z); 
set(i,'EdgeColor','g'); 
hold on 
%Campo de desplazamientos 
quiver3(x,y,z,a,b,c); 
axis equal 
axis([-3,3,-1,3]); 
view(0,90) 
title('Placa y campo de desplazamientos') 
hold off

6 Sólido antes y después del desplazamiento

El desplazamiento del sólido dado por el campo de vectores [math]\vec u[/math]. Lo hemos representado antes y después de que se produzca, así como una comparación entre ambos. A pesar de que se indica que las imágenes sean ambas en una misma figura, lo hemos representado en tres para una mejor visualización.

Figura 6a
Figura 6b
Figura 6c

Aunque en la última imagen no se aprecia bien, el desplazamiento sigue la dirección del campo de vectores (desplazamientos) del apartado anterior. 

%DESPLAZAMIENTO ANTES Y DESPUÉS:
%parametrización u,v
h=50; %paso de muestro
u = linspace(1,2,h); %valor de u en el intervalor [1,2]según el muestro
v = linspace(0,pi,h); %valor de v en el intervalor [0,pi]según el muestro
[U,V] = meshgrid(u,v); %matrices de las componentes polares de la placa
x = U.*cos(V);
y = U.*sin(V);
%Definición de las funciones que describen los desplazamientos que tienen
%lugar en la placa tomando como incógnitas x e y
Ux=(1/5).*(U-1).*cos(V)+x;
Uy=(1/5).*(U-1).*sin(V)+y;
%Representación gráfica:
%ANTES
subplot(1,3,1)
surf(x,y,0*x);
view(2)
axis equal
axis([-3 3 -1 3])
title('Antes del desplazamiento');
%DESPUÉS
subplot(1,3,2)
surf(Ux,Uy,0*Ux);
view(2)
axis equal
axis([-3 3 -1 3])
title('Después del desplazamiento')
%COMPARACIÓN
subplot(1,3,3)
hold on
plot3(x,y,x*0);
plot3(Ux,Uy,0*Ux);
view(2)
axis equal
axis([-3 3 -1 3])
title ('comparación')
hold off

7 Divergencia del campo

Primero calculamos analíticamente la divergencia de [math]\vec u[/math]: centro

La divergencia es mínima en [math]ρ=1[/math] y máxima en [math]ρ=2[/math]. Esto se puede calcular intuitivamente viendo la ecuación de divergencia en coordenadas cilíndricas. En nuestra figura, la divergencia no se anula. Pasamos la divergencia a coordenadas cartesianas, para introducirlo en Matlab: [math]div(\vec u)= \frac {1} {5} *(2-\frac {1} {\sqrt {x^2+y^2}})[/math]

Figura 7
%Divergencia de u: donde sufre más deformación el cuerpo: 2
%parametrización u,v
h=10; %paso de muestro
u = linspace(1,2,h); 
v = linspace(0,pi,h);
[U,V] = meshgrid(u,v);  
x = U.*cos(V);
y = U.*sin(V);
 
%Definición de la divergencia del campo de desplazamientos u tomando como incógnitas x e y
DIVu=(1/5).*(2-(1./(sqrt(x.^2+y.^2))));
 
%Representación gráfica de la divergencia
surf(x,y,DIVu)
view (2)
axis equal
axis([-3,3,-1,3])
colorbar
 
%obtención de la divergencia máxima
DIVuMAX=max(max(DIVu))
%obtención de la divergencia mínima
DIVuMIN=min(min(DIVu))

Los resultados que da Matlab son los siguentes: 

DIVuMAX=0.3000
DIVuMIN=0.2000

La divergencia es una medida del cambio de volumen local debido al desplazamiento, lo apreciamos en la siguiente gráfica: centro

8 Rotacional de [math]\vec u[/math]

Cálculo del valor absoluto del rotacional del vector [math]\vec u[/math]: [math]|\nabla×\vec u|[/math]

centro
Figura 8
u=1:0.1:2;  
h=100;  
v=linspace(0,pi,h);  
[U,V ]=meshgrid(u,v); % matriz de las componentes 
x = U.*cos(V); 
y = U.*sin(V); 
Rot= 0.*x+0.*y; 
surf (x,y,Rot) 
view(2) 
axis equal 
axis ([-3,3,-1,3]) 
Colorbar

¿Qué puntos sufren mayor rotacional? En nuestro caso, partiendo del valor nulo obtenido del rotacional, consideramos que ningún punto sufre la tendencia a ser girado por el rotacional.

9 Estudio del tensor de tensiones

La parte simétrica del tensor gradiente del vector de u, o también conocido como tensor de deformaciones se definirá a partir de la siguiente fórmula:

centro

En un medio elástico lineal, isótopo y homogéneo, como el de nuestra placa, los desplazamientos permiten escribir el tensor de tensiones a través de la siguiente ecuación:

centro

Para nuestro caso supondremos las coeficientes de Lamé ([math]μ[/math] y [math]λ[/math]) de valor [math]1[/math].

Figura 9
%parametrización u,v 
u=1:0.1:2;  
h=100;  
v=linspace(0,pi,h);  
[U,V] = meshgrid(u,v); %matrices de las componentes polares de la placa 
x = U.*cos(V); 
y = U.*sin(V); 
%Representación de la función normal 
f=1/5.*(4-1./(sqrt(x.^2+y.^2))); 
subplot(1,3,1) 
surf(x,y,f) 
axis equal 
%Representación de la función normal 
f=1/5.*(2-1./(sqrt(x.^2+y.^2))); 
subplot(1,3,2) 
surf(x,y,f) 
axis equal 
%Representación de la función normal 
f=1/5.*(2-1./(sqrt(x.^2+y.^2))); 
subplot(1,3,3) 
surf(x,y,f) 
axis equal

10 Tensiones tangenciales

El cálculo de las tensiones tangenciales de nuestra placa para la dirección [math]\vec g_p[/math] es:

Figura 10

La tensiones tangenciales son nulas, esto se debe a que nuestro campo de desplazamientos depende solo de la distancia al origen como se puede apreciar en la figura 5. Si el campo no depende del ángulo la tensión tangencial será siempre nula, debido a que va a ser igual en todas las direcciones.

11 Tensión de Von Mises

La tensión de Von Mises viene dada por la fórmula siguiente:

centro

Las tensiones principales, [math]σ[/math] autovalores, de [math]σ´[/math]se representan como [math]σ´_1[/math], [math]σ´_2[/math] y [math]σ´_3[/math]. Esta formula nos da una magnitud escalar que se suele usar como indicador para saber cuándo un material inicia un comportamiento plástico.

Figura 11
Figura 12
h=0.1; % paso de muestreo para u y v 
u=1:h:2; %rhó 
v=0:h:pi; %theta 
[U,V]=meshgrid(u,v); % las matrices U y V 
x = U.*cos(V); 
y = U.*sin(V); 
lu=length(u); 
lv=length(v); 
%se crean las componentes de la matriz sigma. 
AA= inline('1/5.*(4-1./(sqrt(x.^2+y.^2)))','x','y'); 
BB= inline('1/5.*(4-3./(sqrt(x.^2+y.^2)))','x','y'); 
CC= inline('1/5.*(2-1./(sqrt(x.^2+y.^2)))','x','y'); 
AB= inline('(0.*x+0.*y)','x','y'); 
Msig=[]; 
% hallar los autovalores. 
for i=1:lv  
    for j=1:lu 
        X=AA(x(i,j),y(i,j)); 
        Y=BB(x(i,j),y(i,j)); 
        Z=CC(x(i,j),y(i,j)); 
        O=AB(x(i,j),y(i,j)); 
        Msig=[X,O,0;O,Y,0;0,0,Z]; 
        autovalores=eig(Msig);  
        VM=sqrt(((autovalores(1)-autovalores(2))^2+((autovalores(2)-autovalores(3))^2+((autovalores(3)-autovalores(1))^2))*1/2));  
        G(i,j)=VM; 
end 
end 
%dibujo de la tensión de Von Mises. 
surf(x,y,G) 
title('Tensión de Von Mises') 
colorbar

En nuestro dibujo se alcanza el máximo en todos los puntos que están situados en 0.4

12 Masa total de la planta

Dada la densidad de la placa en coordenadas cartesianas: [math] d(x,y) = 1 + x·y·log(1 + x + y^2) [/math]
Calcularemos la masa total de la placa aproximando la integral numéricamente. Además, la densidad se hace negativa, por lo que tomaremos la función en valor absoluto.
La masa viene dada por la siguiente integral: [math]m=\int_{S} d·dS=\int_{0}^{π}\int_{1}^{2}d(u,v)|\vec n(u,v)|dudv=\int_{0}^{π}\int_{1}^{2}d(u,v)|\vec r_u(u,v)×\vec r_v(u,v)|dudv[/math] 

h=10; %pasos de muestreo
u = linspace(1,2,h); %valor de u 
v = linspace(0,pi,h); %valor de v 
[U,V ]=meshgrid(u,v); %matriz de la PLACA EN COORDENADAS POLARES
x = U.*cos(V); %coordenadas cartesianas
y = U.*sin(V); %coordenadas cartesianas
d=(1+x.*y.*abs(log(1+x+y.^2))).*sqrt(x.^2+y.^2); %funcion de densidad d(x,y)*area placa
a=h^2*d; %nos da la masa en cada punto
mas=sum(sum(a)) %este comando nos suma los elementos de la matriz a la masa total
Masa= 1.5163e+04

La masa total es [math]1.5163·10^4[/math] unidades de masa.

Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Visualización_de_campos_en_elasticidad_en_una_placa_plana_semicircular&oldid=45591»