Diferencia entre revisiones de «Campos y deformaciones en 2D (Grupo 23B)»

Revisión actual del 13:30 12 dic 2023

El trabajo correspondiente a nuestro grupo es el número 5, que consiste en el estudio de campos y deformaciones en dos dimensiones. Para llevarlo a cabo, nos ayudaremos del lenguaje de programación M y del programa informático MATLAB, que nos permitirá entender los cálculos de una manera más visual.

Consideraremos una placa plana que ocupa un cuarto de un anillo circular centrado en el origen y comprendido entre los radios 1 y 2, en el plano y ≥ |x|, por lo que trabajaremos en coordenadas cilíndricas. Físicamente representa la sección transversal de un sólido para el cual se desprecian las variaciones en la dirección ortogonal a la sección considerada.

En ella, vamos a suponer que tenemos definidas dos cantidades físicas: la temperatura [math] T(x,y) [/math], que viene dada por [math] T(x,y) = x^2 + (y-1)^2 [/math], y los desplazamientos [math] \vec u(\rho,\theta) [/math], producidos por la acción de una fuerza [math] \vec u(\rho,\theta)=\frac{\rho}{20}\sin(6(\theta-\frac{\pi}{4}))\vec e_{\rho } [/math].

De esta forma, si definimos [math] \vec r_{0}(x,y)=x \vec i+y\vec j [/math] como el vector posición de los puntos de la placa antes de la deformación, la posición de cada punto [math] (x,y) [/math] de la placa después de la deformación vendrá dada por [math] \vec r_{d}(x,y)=\vec r_{0}(x,y)+\vec u(x,y) [/math].

1 Dibujo del sólido

Comenzaremos dibujando un mallado que represente los puntos interiores del sólido, tomando los ejes en el cuadrado [-3; 3] × [-1; 3] y como paso de muestreo h = 1/20 para las variables x e y.

REPRESENTACIÓN DEL MALLADO

%limpieza de programas anteriores
clear
clc
%mallado interior de la placa rectangular
u=1:0.05:2;
v=pi/4:0.05:3/4*pi;
[ro,teta]=meshgrid(u,v);
%cambio de coordenadas cilíndricas a cartesianas
Mx=ro.*cos(teta);
My=ro.*sin(teta);
%Mz tiene que ser una matriz nula del mismo tamaño que My
mesh(Mx,My,My*0);
axis([-3,3,-1,3]);
view(2)
title('Representación del mallado')
xlabel ('Eje x')
ylabel ('Eje y')

2 Dibujo de las curvas de nivel

Ahora dibujaremos también las curvas de nivel de la temperatura, que viene dada por el campo escalar [math] T(x,y) = x^2 + (y-1)^2 [/math]. Podemos observar la variación de los colores en función de la temperatura. En la zona más fría se emplean tonos azules, mientras que en la zona más cálida se utilizan tonos amarillentos. En la gráfica, encontramos los puntos aproximados de máxima temperatura: [-1.4,1.4] y [1.4,1.4].

CAMPO DE TEMPERATURAS Y CURVAS DE NIVEL

clear
clc
u=1:0.05:2;
v=pi/4:0.05:3/4*pi;
[ro,teta]=meshgrid(u,v);
Mx=ro.*cos(teta); 
My=ro.*sin(teta);
mesh(Mx,My,My*0);

%campo escalar de la temperatura
T=Mx.^2+(My-1).^2;

%este gráfico se coloca en la fila de arriba
subplot(2,1,1)
%aplicamos la función al mallado con un degradado de colores
surf(Mx,My,T);
axis([-3,3,-1,3]);
%para verlo desde arriba
view(2)
title('Campo de temperaturas')
%este gráfico se coloca en la fila de abajo
subplot(2,1,2)
%dibujamos las curvas de nivel
contour(Mx,My,T);
axis([-3,3,-1,3]);
title('Curvas de nivel')
%barra de indicación de colores
colorbar

3 Cálculo y dibujo del gradiente

A continuación, calcularemos el gradiente y lo representaremos. El gradiente se define como:

En nuestro caso, queda:

Podemos ver gráficamente cómo los vectores resultantes son ortogonales a las curvas de nivel ahora ya conocidas. La razón de esto es que el gradiente muestra la dirección máxima de variación en cada punto.

DIBUJO DEL GRADIENTE

clear
clc
u=1:0.1:2;
v=pi/4:0.1:3/4*pi;
[ro,teta]=meshgrid(u,v);
Mx=ro.*cos(teta); 
My=ro.*sin(teta);
mesh(Mx,My,My*0);
T=Mx.^2+(My-1).^2;

contour(Mx,My,T);
axis([-3,3,-1,3]);
colorbar

%para que se inserten los gráficos de las curvas de nivel y el gradiente en un mismo gráfico
hold on
%vectores del gradiente
Mxx=2*Mx;
Myy=2*(My-1);
%campo vectorial en 2D
quiver(Mx,My,Mxx,Myy);
%centramos los ejes
axis equal
axis([-3,3,-1,3]);
view(2)
title('Gradiente de T')
hold off

4 Dibujo del campo de vectores

Posteriormente, vamos a dibujar el campo de vectores en los puntos del mallado de la placa, tomando el campo dado en el enunciado. Se ha realizado un cambio de coordenadas cilíndricas a cartesianas, conociendo la expresión del vector [math] \vec e_{\rho} [/math]:

El campo en cilíndricas es el siguiente:

El campo en cartesianas nos queda:

DIBUJO DEL CAMPO

clear
clc
u=1:0.05:2;
v=pi/4:0.05:3/4*pi;
[ro,teta]=meshgrid(u,v);
Mx=ro.*cos(teta); 
My=ro.*sin(teta);
mesh(Mx,My,My*0);
axis([-3,3,-1,3]);
view(2)
title('Campo u')
xlabel ('Eje x')
ylabel ('Eje y')

hold on
mx=(ro.*sin(6.*(teta-pi/4)).*cos(teta))/20;
my=(ro.*sin(6.*(teta-pi/4)).*sin(teta))/20;
%campo vectorial en 2D
quiver(Mx,My,mx,my);
axis([-3,3,-1,3]);
hold off

5 Dibujo antes y después del desplazamiento

Seguidamente, representaremos el sólido antes y después del desplazamiento dado por el campo de vectores [math] \vec u [/math]. En vista de su expresión, el campo de desplazamientos sólo varía en la dirección del vector unitario [math] \vec e_{\rho} [/math].

SÓLIDO ANTES DEL DESPLAZAMIENTO. SÓLIDO DESPUÉS DEL DESPLAZAMIENTO. COMPARACIÓN DEL ANTES Y DEL DESPUÉS

clear
clc
u=1:0.1:2;
v=linspace(pi/4,3/4*pi,20);
[ro,teta]=meshgrid(u,v);
Mx=ro.*cos(teta); 
My=ro.*sin(teta);

%situación inicial
subplot(1,3,1)
surf(Mx,My,My*0);
axis([-3,3,-1,3]);
view(2)
title('Situación inicial')

%situación final
subplot(1,3,2)
mx=(ro.*sin(6.*(teta-pi/4)).*cos(teta))/20;
my=(ro.*sin(6.*(teta-pi/4)).*sin(teta))/20;
rx=Mx+mx;
ry=My+my;
surf(rx,ry,ry*0);
axis([-3,3,-1,3]);
view(2)
title('Situación final')

%comparación
subplot(1,3,3)
plot3(Mx,My,My*0);
hold on
plot3(rx,ry,ry*0);
axis([-3,3,-1,3]);
view(2)
title('Comparación')
hold off

6 Cálculo de la divergencia

A continuación, calcularemos la divergencia del campo [math] \vec u [/math], que es una medida del cambio de volumen local debido al desplazamiento. Su expresión es:

En nuestro caso, queda:

De este modo, deducimos que si un cuerpo se desplaza en un medio elástico y la divergencia del campo correspondiente es cero, esto quiere decir que no ha cambiado su volumen. En este caso, como obtenemos una divergencia distinta de cero, afirmamos que el volumen del sólido se ve afectado por el movimiento de sus moléculas.

En la siguiente gráfica, podemos ver que los puntos en los que la divergencia es mínima son aquellos que cumplen:

Los puntos en los que la divergencia es máxima vienen dados de dos formas:

Por otro lado, los puntos en los que la divergencia es nula quedan definidos por:

DIBUJO DE LA DIVERGENCIA

clear
clc
u=1:0.05:2;
v=pi/4:0.05:3/4*pi;
[ro,teta]=meshgrid(u,v);
Mx=ro.*cos(teta); 
My=ro.*sin(teta);

%representamos la divergencia
div=(sin(6.*(teta-pi/4)))/10;
surf(Mx,My,div);
axis([-3,3,-1,3]);
view(2)
title('Divergencia')
colorbar
%para mostrar un degradado de colores
shading interp

7 Cálculo y dibujo del rotacional

En este apartado, calcularemos el rotacional del campo [math] \vec u [/math], que viene dado por la siguiente expresión:

En nuestro caso, queda:

Mientras que el rotacional se representa como un campo vectorial, el módulo del rotacional se trata de un campo escalar:

El rotacional nos indica la capacidad que tienen los distintos puntos de nuestra placa para girar sobre otro punto determinado. En la siguiente gráfica, podemos observar que los puntos que sufren un mayor rotacional son aquellos representados en color amarillo. Estos puntos vienen dados por:

DIBUJO DEL MÓDULO DEL ROTACIONAL

clear
clc
u=1:0.05:2;
v=pi/4:0.05:3/4*pi;
[ro,teta]=meshgrid(u,v);
Mx=ro.*cos(teta); 
My=ro.*sin(teta);

%módulo del rotacional
rot=abs(cos(6.*(teta-pi/4)).*(3/10)); 
%dibujamos el rotacional
surf(Mx,My,rot); 
axis([-3,3,-1,3]);
view(2)
title('Módulo del rotacional')
colorbar
shading interp

8 Dibujo de las tensiones normales

A continuación, analizaremos las tensiones normales a las que se ve sometida la placa. En un medio elástico lineal, isótropo y homogéneo, los desplazamientos permiten escribir el tensor de tensiones [math]\sigma_{\rho \, \theta}[/math] a través de la fórmula:

donde:

• 1 es el tensor identidad en el conjunto de vectores libres del espacio [math] \mathbb{R}^{3}[/math]
• λ, µ son los conocidos como «coeficientes de Lamé», que dependen de las propiedades elásticas de cada material

Definimos [math]\epsilon (\vec u)=(\nabla\vec u+ \nabla \vec u^{t})/2 [/math] como la parte simétrica del tensor gradiente de [math] \vec u [/math] conocido como «tensor de deformaciones». Tomando λ=μ=1 debido a las propiedades de nuestro material, calculamos las tensiones normales en la dirección que marca el eje [math]\vec e_\rho [/math], en la dirección que marca el eje [math]\vec e_\theta [/math] y las correspondientes al eje [math]\vec e_z [/math], y las representamos en el dibujo.

8.1 CÁLCULO DEL TENSOR DE DEFORMACIONES Y DE TENSIONES

Para calcular el tensor deformaciones [math]\nabla{\vec u(ρ,θ,z)}[/math], necesitamos calcular:

• [math]\frac{\partial \vec u}{\partial \rho}=\frac{\partial}{\partial \rho}(u_{\rho})\vec e_{\rho}+u_{\rho}\frac{\partial}{\partial \rho}(\vec e_{\rho})=\frac{sin(6(\theta-\pi/4))}{20}\vec e_{\rho}+u_{\rho} \vec 0=\frac{sin(6(\theta-\pi/4))}{20}\vec e_{\rho}[/math]

• [math]\frac{\partial \vec u}{\partial \theta}=\frac{\partial}{\partial \theta}(u_{\rho})\vec e_{\rho}+u_{\rho}\frac{\partial}{\partial \theta}(\vec e_{\rho})=6 \rho\frac{cos(6(\theta-\pi/4))}{20}\vec e_{\rho}+\rho \frac{sin(6(\theta-\pi/4))}{20} \vec e_{\theta}[/math]

• [math]\frac{\partial \vec u}{\partial z}=\frac{\partial}{\partial z}(u_{\rho})\vec e_{\rho}+u_{\rho}\frac{\partial}{\partial z}(\vec e_{\rho})=0\vec e_{\rho}+\rho \frac{sin(6(\theta-\pi/4))}{20} \vec 0=\vec 0[/math]

Ahora sí, calculamos [math]\nabla{\vec u(ρ,θ,z)}[/math] y [math](\nabla{\vec u(ρ,θ,z)})^t[/math]:

• [math]\nabla{\vec u(ρ,θ,z)}= \left(\begin{matrix} \frac{ sin (6(\theta-\pi/4))}{20} & \frac{ 6 cos (6(\theta-\pi/4))}{20} & 0 \\ 0 & \frac{ sin (6(\theta-\pi/4))}{20} & 0 \\ 0 & 0 & 0 \end{matrix}\right)[/math]

• [math](\nabla{\vec u(ρ,θ,z)})^t = \left(\begin{matrix} \frac{ sin (6(\theta-\pi/4))}{20} & 0 & 0 \\ \frac{ 6 cos (6(\theta-\pi/4))}{20} & \frac{ sin (6(\theta-\pi/4))}{20}& 0 \\ 0 & 0 & 0 \end{matrix}\right)[/math]

Sabemos que [math] \lambda = \mu =1 [/math] y que [math] \nabla \cdot \vec u(\rho,\theta,z)=\frac{1}{10}\sin(6(\theta-\frac{\pi}{4})) [/math], por lo que el resultado del tensor de tensiones sería:

8.2 TENSIONES NORMALES EN LA DIRECCIÓN DEL EJE [math] \vec e_{\rho} [/math]

Las tensiones normales en la dirección que marca el eje [math]\vec e_{\rho}[/math] vienen dadas por:

Podemos ver que la solución se corresponde con la componente (1,1) de la matriz del tensor de tensiones.

TENSIONES NORMALES AL EJE [math]\vec e_{\rho}[/math] EN 2D

TENSIONES NORMALES AL EJE [math]\vec e_{\rho}[/math] EN 3D

%tensiones normales al eje e-ro
clear
clc
u=1:0.05:2;
v=pi/4:0.05:3/4*pi;
[ro,teta]=meshgrid(u,v);
Mx=ro.*cos(teta); 
My=ro.*sin(teta);

%elemento (1,1) de la matriz sigma
sigma1=(sin(6.*(teta-pi/4)).*(1/5));
surf(Mx,My,sigma1);
axis([-3,3,-1,3]);
view(2)
title('Tensiones normales en la dirección del eje e-ro')
colorbar
shading interp

8.3 TENSIONES NORMALES EN LA DIRECCIÓN DEL EJE [math] \vec e_{\theta} [/math]

Las tensiones normales en la dirección que marca el eje [math]\vec e_{\theta}[/math] vienen dadas por:

Podemos ver que la solución se corresponde con la componente (2,2) de la matriz del tensor de tensiones.

TENSIONES NORMALES AL EJE [math]\vec e_{\theta}[/math] EN 2D

TENSIONES NORMALES AL EJE [math]\vec e_{\theta}[/math] EN 3D

%tensiones normales al eje e-teta
clear
clc
u=1:0.05:2;
v=pi/4:0.05:3/4*pi;
[ro,teta]=meshgrid(u,v);
Mx=ro.*cos(teta); 
My=ro.*sin(teta);

%elemento (2,2) de la matriz sigma
sigma2=(sin(6.*(teta-pi/4)).*(1/5));
surf(Mx,My,sigma2);
axis([-3,3,-1,3]);
view(2)
title('Tensiones normales en la dirección del eje e-teta')
colorbar
shading interp

8.4 TENSIONES NORMALES EN LA DIRECCIÓN DEL EJE [math] \vec e_{z} [/math]

Las tensiones normales en la dirección que marca el eje [math]\vec e_{z}[/math] vienen dadas por:

Podemos ver que la solución se corresponde con la componente (3,3) de la matriz del tensor de tensiones.

TENSIONES NORMALES AL EJE [math]\vec e_z[/math] EN 2D

TENSIONES NORMALES AL EJE [math]\vec e_z[/math] EN 3D

%tensiones normales al eje e-z
clear
clc
u=1:0.05:2;
v=pi/4:0.05:3/4*pi;
[ro,teta]=meshgrid(u,v);
Mx=ro.*cos(teta); 
My=ro.*sin(teta);

%elemento (3,3) de la matriz sigma
sigma3=(sin(6.*(teta-pi/4)).*(1/10));
surf(Mx,My,sigma3);
axis([-3,3,-1,3]);
view(2)
title('Tensiones normales en la dirección del eje e-z')
colorbar
shading interp

A pesar de que los desplazamientos se realizan en el plano, las tensiones no tienen por qué ser planas y pueden existir en la dirección ortogonal al plano de la placa. En este caso, podemos ver que la tensión es más elevada en las direcciones [math] \vec e_{\rho} [/math] y [math] \vec e_{\theta} [/math] que en la dirección [math] \vec e_{z} [/math].

9 Cálculo de las tensiones tangenciales

En este apartado, calcularemos las tensiones tangenciales respecto al plano ortogonal a [math] \vec e_{\rho} [/math], es decir, [math]\left | \sigma\cdot \vec e_{\rho}-(\vec e_{\rho}\cdot \sigma \cdot\vec e_{\rho}) \vec e_{\rho} \right |[/math].

Podemos ver que la solución se corresponde con la componente (2,1) de la matriz del tensor de tensiones.

TENSIONES TANGENCIALES RESPECTO AL PLANO ORTOGONAL [math] \vec e_{\rho}\ [/math] EN 3D

%cálculo de las tensiones tangenciales 
clear
clc
u=1:0.05:2;
v=pi/4:0.05:3/4*pi;
[ro,teta]=meshgrid(u,v);
Mx=ro.*cos(teta);
My=ro.*sin(teta);

erho=abs(6*cos(6*(teta-pi/4))/20);
surf(Mx,My,erho);
axis equal
title('Tensiones tangenciales respecto al plano ortogonal a  e-rho')
colorbar
shading interp

10 Dibujo de la tensión de Von Mises

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

donde [math] \sigma_1 [/math], [math] \sigma_2 [/math] y [math] \sigma_3 [/math] son los autovalores de [math] \sigma [/math], también conocidos como «tensiones principales». Se trata de una magnitud escalar que se suele usar como indicador para saber cuándo un material inicia un comportamiento plástico. Para calcularlos, hemos empleado el comando eig en MATLAB.

En las siguientes gráficas, podemos observar que los puntos en que [math] \sigma_{VM} [/math] alcanza su mayor valor son aquellos representados en color amarillo. Este valor es 0.5614 (calculado por MATLAB) y se da en los puntos cuyas coordenadas cumplen una de las siguientes condiciones:

TENSIÓN DE VON MISES EN 3D

TENSIÓN DE VON MISES EN 2D

clear
clc
u=1:0.05:2;
v=pi/4:0.05:3/4*pi;
[ro,teta]=meshgrid(u,v);
Mx=ro.*cos(teta); 
My=ro.*sin(teta);

lu=length(u); 
lv=length(v);
%creamos las componentes de la matriz sigma
AA=inline('(1/10.*(sin(6.*(atan(My./Mx)-pi/4))))','Mx','My');
BB=inline('(1/5.*(sin(6.*(atan(My./Mx)-pi/4))))','Mx','My');
CC=inline('(1/10.*(sin(6.*(atan(My./Mx)-pi/4))))','Mx','My');
AB=inline('(6/20.*(cos(6.*(atan(My./Mx)-pi/4))))','Mx','My');
Msig=[];
%hallamos los autovalores
for i=1:lv
  for j=1:lu
    X=AA(Mx(i,j),My(i,j));
    Y=BB(Mx(i,j),My(i,j));
    Z=CC(Mx(i,j),My(i,j));
    O=AB(Mx(i,j),My(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
%dibujamos la tensión de Von Mises
surf(Mx,My,G);
title('Tensión de Von Mises')
colorbar
shading interp
%valor máximo de la tensión de Von Mises
max=max(max(G))

11 Cálculo y dibujo del campo de fuerzas

El campo de fuerzas [math] \vec F [/math], que actúa sobre la placa y que es el causante del desplazamiento observado, se calcula con la ecuación de Lamé:

donde [math] \Delta \vec u [/math] es el laplaciano del campo [math] \vec u [/math] y [math] \nabla(\nabla \cdot \vec u) [/math] es el gradiente de la divergencia del campo [math] \vec u [/math]. Por otro lado, sabemos que el campo [math] \vec u [/math] no depende del tiempo y que [math] \lambda=\mu=1 [/math] por las propiedades del material, por lo que la ecuación queda:

Ya conocemos el valor de la divergencia del campo [math] \vec u [/math], calculado en el apartado 6. Así, en coordenadas cartesianas, el gradiente de la divergencia nos quedaría:

Por otro lado, el laplaciano del campo vectorial [math] \vec u [/math] ha sido calculado, en coordenadas cartesianas, de la siguiente manera:

Tras introducir estas dos expresiones en la ecuación del campo de fuerzas [math] \vec F [/math], obtenemos el siguiente código M que nos permite representarlo:

CAMPO DE FUERZAS F QUE ACTÚA SOBRE LA PLACA

clear
clc
u=1:0.05:2;
v=pi/4:0.05:3/4*pi;
[ro,teta]=meshgrid(u,v);
Mx=ro.*cos(teta);
My=ro.*sin(teta);

fuerzai=(6*cos(6*(teta-pi/4)).*sin(teta))./(5*ro)-((sin(6*(teta-pi/4)).*cos(teta)/20)+((-36*sin(6*(teta-pi/4)).*cos(teta)-6*cos(6*(teta-pi/4)).*sin(teta)-6*cos(6*(teta-pi/4)).*sin(teta)-sin(6*(teta-pi/4)).*cos(teta))/20))./ro;
fuerzaj=(-6*cos(6*(teta-pi/4)).*cos(teta))./(5*ro)-((sin(6*(teta-pi/4)).*cos(teta)/20)+((-36*sin(6*(teta-pi/4)).*sin(teta)+6*cos(6*(teta-pi/4)).*cos(teta)+6*cos(6*(teta-pi/4)).*cos(teta)-sin(6*(teta-pi/4)).*sin(teta))/20))./ro;
quiver(Mx,My,fuerzai,fuerzaj);
axis([-3,3,-1,3]);
title('Campo de fuerzas F')
xlabel ('Eje x')
ylabel ('Eje y')