Visualización de campos escalares y vectoriales en elasticidad (Grupo 4B)

De MateWiki
Saltar a: navegación, buscar
Trabajo realizado por estudiantes
Título Deformaciones de un sólido de sección semicircular. Grupo 4-B
Asignatura Teoría de Campos
Curso 2020-21
Autores Miguel Ángel Díaz Hernández
John Cuenca Uyaguari
Jesús Navarro Amador
Antoni Capó Villalonga
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura


El trabajo a realizar consta del estudio de la transformación que un sólido sufre.Se trata de un desplazamiento.

Se tienen definidas dos cantidades 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 determinada. Para dicho campo observaremos su gradiente, divergencia y rotacional.


Más adelante veremos las tensiones a las que el cuerpo ha sido sometido

1 Introducción del sólido


Suponemos el siguiente cuerpo. Para poder observar las fuerzas de desplazamiento que ocurren sobre este tomaremos una sección conveniente. En este caso es preferible cortar por un plano perpendicular al eje del cilindro.
El resultado de la sección es el anillo circular centrado en el origen y comprendido entre los radios [math]1[/math] y [math]2[/math], es decir [math]1≤x^2+y^2≤4[/math], y en el plano [math]y≥0[/math], que en coordenadas cilíndricas equivale a la región [math](ρ,θ)∈[1,2]×[0,π][/math].

left Placamallado4b.png 


























%Inicializamos las variables
rho = [1:0.1:2];
h = round(pi/0.1);
theta= linspace(0,pi,h);

%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta); 

X=RHO.*cos(THETA);        
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)  
axis([-3,3,-1,3]); 

view(2);
title('Mallado y placa')



2 Temperatura y gradiente


Pasamos a analizar la primera de las cantidades físicas del proyecto, la temperatura \(T(ρ,θ)\)
Ésta viene dada por el campo escalar:

[math]T(x,y)=ln(y^2+2)[/math].


Así mismo, el gradiente de la temperatura \(T(x,y)\) es el campo vectorial:

[math]\nabla T(x,y) =\frac{2y}{y^2+2} \vec j[/math].


Temperatura:
Temperatura4b.png
La imagen de la izquierda muestra las curvas de nivel, en las cuales la temperatura es constante.

En la imagen de la derecha, podemos ver como la temperatura crece sobre la placa y va aumentando a medida que subimos por esta, hasta alcanzar su máximo, con una temperatura de [math] 1,8 [/math]°C ,

Gradiente:
izquierda derecha


























Para visualizar el gradiente de la temperatura \(T(x,y)\) en MatLab hemos creado tres subventanas. En las vistas 2D podemos apreciar la dirección del campo según [math]\vec j[/math] y una vista en planta del mismo
El gradiente nos indica la dirección en la que la temperatura varía más rápidamente.
En nuestro caso, podemos observar que a más altura, mayor es el gradiente, y por tanto, más varía la temperatura sobre la placa. Se puede apreciar también la ortogonalidad del gradiente respecto a las curvas de nivel.

Código matlab temperatura 

%Inicializamos las variables
rho = [1:0.1:2];
h = round(pi/0.1);
theta= linspace(0,pi,h);

%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta); 

X=RHO.*cos(THETA);        
Y=RHO.*sin(THETA);

%CAMPO ESCALAR DE TEMPERATURAS
Temperatura=log((Y.^2)+2);
Maximo=max(max(Temperatura));

%Representación de la malla en 2D en la subventana primera
subplot(1,2,1)
hold on

view(2)
axis([-3,3,-1,3]);
axis equal;
title('Lineas de nivel de temperatura en 2D')
contour(X,Y,Temperatura,10)
hold off
subplot(1,2,1)

%Representación de la malla en 3D en la subventana segunda
subplot(1,2,2);

surf(X,Y,Temperatura);
colorbar;
axis([-3,3,-1,3]); 
axis vis3d
title('Temperatura en 3D')


Código matlab gradiente 

%Inicializamos las variables
rho = [1:0.3:2];
h = round(pi/0.2);
theta= linspace(0,pi,h);

%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta); 

X=RHO.*cos(THETA);        
Y=RHO.*sin(THETA);

%GRADIENTE
Grad=(2*Y)./((Y.^2)+2);

%CAMPO ESCALAR DE TEMPERATURAS
Temperatura=log((Y.^2)+2);

%Creación de la subventana primera: gradiente en el campo de temperaturas

subplot(2,2,1)
hold on
surf(X,Y,Temperatura);
axis([-3,3,-1,3]);
%Dibujo del gradiente
quiver3(X,Y,Temperatura,0*X,Grad,0*Temperatura);
axis([-3,3,-1,3]);
axis vis3d
view([-60,10])
title('Gradiente 3D')
hold off

subplot(2,2,3)
hold on
surf(X,Y,Temperatura);
axis([-3,3,-1,3]);
%Dibujo del gradiente
quiver3(X,Y,Temperatura,0*X,Grad,0*Temperatura);
axis([-3,3,-1,3]);
axis vis3d
view([-90,0])
title('Gradiente 2D (vista lateral)')
hold off
%Creación de la segunda subventana: campo vectorial del gradiente

subplot(2,2,2)
quiver(X,Y,0*X,Grad); 
colorbar
axis([-3,3,-1,3]);
title('Gradiente 2D')
hold on

view(2)
axis([-3,3,-1,3]);
axis equal;
title('Lineas de nivel de temperatura en 2D')
contour(X,Y,Temperatura,10)
hold off


3 Calculo de campo de desplazamientos


Ahora tendremos que evaluar la otra magnitud. Para hallar el campo de desplazamientos [math] \vec u(ρ,θ) [/math] deberemos tener en cuenta la siguiente ecuación:

  • [math]\vec u(ρ,θ) = senθ f(ρ) \vec g_θ [/math]


Tendremos que hallar f(ρ).
Que los puntos en ρ=1 no sufren desplazamiento significa que

  • [math]\vec u(ρ,θ) = senθ f(1) \vec g_θ = 0 ; f(1) = 0 [/math] ;


  • [math]\nabla×\vec u(ρ,θ) = \frac{1}{ρ}\left|\begin{matrix} \vec g_ρ & \vec g_θ & \vec g_z \\ \frac{\partial}{\partial ρ} & \frac{\partial}{\partial θ} & \frac{\partial}{\partial z} \\ 0 & ρ^2f(ρ)senθ & 0 \end{matrix}\right| = \frac{1}{ρ}(ρ^2f(ρ))´senθ · \vec g_z[/math].


Vamos a suponer que la función [math]ρ^2f(ρ) [/math] es creciente, [math] (ρ^2f(ρ))´≥0[/math].

  • [math]|\nabla×\vec u(ρ,θ)| =\frac{1}{ρ}(ρ^2f(ρ))´senθ = \frac{2ρ - 1}{10}senθ[/math].


Procedemos a integrar para resolver la EDO:
[math] (ρ^2f(ρ))´= \frac{2ρ - 1}{10} ↔ ρ^2f(ρ) = \frac{2ρ^3}{30} - \frac{ρ^3}{20} +c ↔ f(ρ) = \frac{2ρ}{30} - \frac{1}{20} + \frac{c}{ρ^2} [/math].
Imponemos la condicón [math] f(1) = 0 [/math]
[math] f(1) = \frac{2}{30} - \frac{1}{20} + c = 0 ↔ c = -\frac{1}{60} [/math]
[math] f(ρ) = \frac{ρ}{15} - \frac{1}{20} + \frac{1}{60ρ^2} [/math]
Finalmente queda que el campo de desplazamientos es:

  • [math] \vec u(ρ,θ) = ( \frac{ρ}{15} - \frac{1}{20} + \frac{1}{60ρ^2} ) senθ \vec g_θ [/math]


4 Campo de desplazamientos


Para visualizar el campo de desplazamientos [math]\vec u(ρ,θ)[/math] en MatLab hemos sustituido el valor de [math] \vec g_θ [/math] :

[math]\vec u(ρ,θ) =(\frac{ρ}{15}-\frac{1}{20}-\frac{1}{60ρ^2}) \vec g_θ = (\frac{ρ}{15}-\frac{1}{20}-\frac{1}{60ρ^2})[-ρsenθ\vec i + ρcosθ \vec j] [/math]


Así mismo observamos las dos componentes del campo vectorial:

[math]Fx = (\frac{ρ}{15}-\frac{1}{20}-\frac{1}{60ρ^2})[-ρsenθ\vec i] [/math]
[math]Fy = (\frac{ρ}{15}-\frac{1}{20}-\frac{1}{60ρ^2})[ρcosθ \vec j][/math]



En la siguiente imagen vemos el desplazamiento que sufrirá la placa.

Campodesp.png 
%Inicializamos las variables
rho = [1:0.1:2];
h = round(pi/0.1);
theta= linspace(0,pi,h);

%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta); 

X=RHO.*cos(THETA);        
Y=RHO.*sin(THETA);

%CAMPO VECTORIAL DE DESPLAZAMIENTOS
%Componente x(X,Y) del campo vectorial
FX=(RHO/15-1/20-1./(60*(RHO.^2))).*(sin(THETA)).*(-Y); 

%Componente y(X,Y) del campo vectorial
FY=(RHO/15-1/20-1./(60*(RHO.^2))).*sin(THETA).*X; 

%Representamos el campo vectorial de desplazamientos
quiver(X,Y,FX,FY)

axis([-3,3,-1,3]);
view(2)

xlabel('Eje x')
ylabel('Eje y')
title('Campo de desplazamientos sobre la placa')


5 Placa antes y después del desplazamiento


centro
En esta imagen podemos comprobar que la zona central de la placa es la que más se desplaza, por el contrario, en la parte inferior de la placa el desplazamiento es mínimo. 

u = [1:0.1:2];
h = round(pi/0.1);
v = linspace(0,pi,h);

[RHO,THETA]=meshgrid(u,v);

x = RHO.*cos(THETA);
y = RHO.*sin(THETA);
z = RHO.*0;


b =(RHO/15-1/20-1./(60*(RHO.^2))).*(sin(THETA)).*(-y);  
c =(RHO/15-1/20-1./(60*(RHO.^2))).*sin(THETA).*x; 

xd = x+b;
yd = y+c;

%PLACA CON DESPLAZAMIENTO
subplot(2,1,2)
j = mesh(xd,yd,z);
set(j,'EdgeColor','r');
hold on


axis([-3 3 -1 3]);
view(0,90)
title('Placa Desplazada');
hold off

%PLACA CON DESPLAZAMIENTO Y CAMPO
subplot(2,1,1)


hold on


m = mesh(x,y,z);
set(m,'EdgeColor','b');
axis([-3 3 -1 3]);
view(0,90)
title('Desplazamiento de la placa');

%CAMPO VECTORIAL DE DESPLAZAMIENTOS
%Componente x(X,Y) del campo vectorial
FX=(RHO/15-1/20-1./(60*(RHO.^2))).*(sin(THETA)).*(-y); 

%Componente y(X,Y) del campo vectorial
FY=(RHO/15-1/20-1./(60*(RHO.^2))).*sin(THETA).*x; 
 

%Representamos el campo vectorial de desplazamientos
quiver(x,y,FX,FY,'k')

hold off


6 Divergencia del campo desplazamientos


6.1 Cálculo de la divergencia


Con la divergencia mediremos la diferencia entre el flujo saliente y el entrante del campo vectorial de desplazamiento. Más adelante en la imagen podremos ver la tendencia que tiene la superficie a desviarse de su posición de partida
La divergencia del campo de desplazamientos [math]\vec u(ρ,θ)[/math] viene dada por el campo escalar:

[math]\nabla · \vec u(ρ,θ) = \frac{1}{ρ} \frac{\partial}{\partial θ}(ρ(\frac{ρ}{15}-\frac{1}{20}-\frac{1}{60ρ^2})senθ)=(\frac{ρ}{15}-\frac{1}{20}-\frac{1}{60ρ^2})cosθ[/math].


6.2 Visualización de la divergencia


Divergencia4b.png Div2.png 

%Inicializamos las variables
%Inicializamos las variables
rho = [1:0.1:2];
h = round(pi/0.1);
theta= linspace(0,pi,h);

%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta); 

%Divergencia
subplot(1,2,1)
X=RHO.*cos(THETA);        
Y=RHO.*sin(THETA);
%Función de divergencia
f=(RHO/15-1/20-1./(60*(RHO.^2))).*cos(THETA);
%Dibujamos la función de divergencia

surf(X,Y,f)
axis([-3,3,-1,3]);
colorbar;
%Vista y título de la gráfica
axis vis3d
title('Divergencia del campo de desplazamiento')

subplot(1,2,2)


surf(X,Y,f)
axis([-3,3,-1,3]);
colorbar;
%Vista y título de la gráfica
view(2)
title('Divergencia del campo de desplazamiento')


6.2.1 Comparación desplazamiento con divergencia


centro

%Inicializamos las variables
rho = [1:0.1:2];
h = round(pi/0.1);
theta= linspace(0,pi,h);

%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta); 

X=RHO.*cos(THETA);        
Y=RHO.*sin(THETA);
z = RHO.*0;
subplot(1,3,1)
%Dibujamos el mallado
mesh(X,Y,0*X)  
axis([-3,3,-1,3]); 
view(2);
title('placa')

b =(RHO/15-1/20-1./(60*(RHO.^2))).*(sin(THETA)).*(-Y);  
c =(RHO/15-1/20-1./(60*(RHO.^2))).*sin(THETA).*X; 

xd = X+b;
yd = Y+c;

%PLACA CON DESPLAZAMIENTO
subplot(1,3,2)
j = mesh(xd,yd,z);
set(j,'EdgeColor','r');
title('placa desplazada')
view(2)
axis([-3,3,-1,3]); 
f=(RHO/15-1/20-1./(60*(RHO.^2))).*cos(THETA);
%Dibujamos la función de divergencia
subplot(1,3,3)
surf(X,Y,f)
axis([-3,3,-1,3]);
title('DIVERGENCIA')


7 Rotacional del campo desplazamientos


7.1 Cálculo del rotacional


El rotacional del campo de desplazamientos [math]\vec u(ρ,θ)[/math] viene dado por el campo vectorial:

[math]\nabla×\vec u(ρ,θ) = \frac{1}{ρ}\left|\begin{matrix} \vec g_ρ & \vec g_θ & \vec g_z \\ \frac{\partial}{\partial ρ} & \frac{\partial}{\partial θ} & \frac{\partial}{\partial z} \\ 0 & ρ^2f(ρ)senθ & 0 \end{matrix}\right| = \frac{1}{ρ}(ρ^2f(ρ))´senθ · \vec g_z[/math].


Por lo tanto su módulo será igual a:

[math]|\nabla×\vec u(ρ,θ)| =\frac{1}{ρ}(ρ^2f(ρ))´senθ = \frac{2ρ - 1}{10}senθ[/math].


7.2 Visualización del módulo del rotacional


La imagen siguiente muestra el valor absoluto del rotacional.
El rotacional nos permite ver la tendencia que tiene un campo vectorial a inducir rotación alrededor de un punto.

Rot1.pngRot2.png






















En nuestro ejercicio, la placa tendrá mayor rotacional en la parte superior central; mientras que la zona inferior no tendrá prácticamente rotacional. 

%Inicializamos las variables
rho = [1:0.1:2];
h = round(pi/0.1);
theta= linspace(0,pi,h);

%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta);
X=RHO.*cos(THETA);  
Y=RHO.*sin(THETA);
%Función del rotacional
rot=abs((2*RHO-1)/10.*sin(THETA));
%Dibujamos la función del  MÓDULO  del rotacional
subplot(1,2,1)
surf(X,Y,rot)
axis([-3,3,-1,3]);
colorbar;
%Vista y título de la gráfica
axis vis3d
title('Rotacional del campo de desplazamiento')

subplot(1,2,2)
surf(X,Y,rot)
axis([-3,3,-1,3]);
colorbar;
%Vista y título de la gráfica
view(2)
title('Rotacional del campo de desplazamiento')


7.2.1 Comparación desplazamiento con rotacional


Comprot.png 


%Inicializamos las variables
rho = [1:0.1:2];
h = round(pi/0.1);
theta= linspace(0,pi,h);

%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta); 

X=RHO.*cos(THETA);        
Y=RHO.*sin(THETA);
z = RHO.*0;
subplot(1,3,1)
%Dibujamos el mallado
mesh(X,Y,0*X)  
axis([-3,3,-1,3]); 
view(2);
title('placa')

b =(RHO/15-1/20-1./(60*(RHO.^2))).*(sin(THETA)).*(-Y);  
c =(RHO/15-1/20-1./(60*(RHO.^2))).*sin(THETA).*X; 

xd = X+b;
yd = Y+c;

%PLACA CON DESPLAZAMIENTO
subplot(1,3,2)
j = mesh(xd,yd,z);
set(j,'EdgeColor','r');
title('placa desplazada')
view(2)
axis([-3,3,-1,3]); 
rot=(2*RHO-1)/10.*sin(THETA);
%Dibujamos la función de divergencia
subplot(1,3,3)
surf(X,Y,rot)
axis([-3,3,-1,3]);
title('rotacional')


8 Tensiones


Para el cálculo de las tensiones que el sólido sufre [math]σ[/math] vamos a tener en cuenta que disponemos del siguiente croquis de un sólido cualquiera
centro
La tensión [math]σ[/math] se debe expresar a través de la sección que el plano produce, nuestra placa. Para dicho plano se toman en cuenta la tensión expresada en dos direcciones. La tensión normal y la tensión tangencial.
Analizaremos las tensiones en tres direcciones, es decir, las expresaremos en la base en la base {[math]\vec g_ρ,\frac{\vec g_θ}{ρ},\vec g_z[/math]}

8.1 Cálculo del tensor de deformaciones y de tensiones


  • [math]\frac{\partial \vec u}{\partial ρ}=(\frac{1}{15}+\frac{1}{30ρ^3})senθ \vec g_θ + (\frac{ρ}{15} - \frac{1}{20} - \frac{1}{60ρ^2})senθ \frac{\partial \vec g_θ }{\partial ρ} = (\frac{1}{15}+\frac{1}{30ρ^3})senθ \vec g_θ + (\frac{ρ}{15} - \frac{1}{20} - \frac{1}{60ρ^2}) )senθ \frac{\partial \vec g_θ }{ρ} =[/math]


[math]=(\frac{ρ}{15}+\frac{1}{30ρ^2})senθ \vec w_θ + (\frac{ρ}{15} - \frac{1}{20} - \frac{1}{60ρ^2}))senθ \vec w_θ [/math].



  • [math]\frac{\partial \vec u}{\partial θ}=(\frac{ρ}{15} - \frac{1}{20} - \frac{1}{60ρ^2})cosθ \vec g_θ + (\frac{ρ}{15} - \frac{1}{20} - \frac{1}{60ρ^2})senθ \frac{\partial \vec g_θ }{\partial θ} = (\frac{ρ}{15} - \frac{1}{20} - \frac{1}{60ρ^2})cosθ \vec g_θ + (\frac{ρ}{15} - \frac{1}{20} - \frac{1}{60ρ^2})senθ \vec g_ρ =[/math]


[math]=(\frac{ρ}{20}+\frac{1}{60ρ} - \frac{ρ^2}{15})senθ \vec w_ρ + (\frac{ρ^2}{15} - \frac{ρ}{20} - \frac{1}{60ρ}))cosθ \vec w_θ [/math].



Se define el tensor de deformaciones [math]e(\vec u(ρ,θ))[/math] como la parte simétrica del tensor [math] \nabla{\vec u(ρ,θ)} [/math], llamado también gradiente del campo de desplazamientos [math]\vec u(ρ,θ) [/math]:

[math]e(\vec u(ρ,θ)) =\frac{\nabla{\vec u(ρ,θ)}+(\nabla{\vec u(ρ,θ)})^t}{2}[/math].


Para ello se calculan [math]\nabla{\vec u(ρ,θ)}[/math] y [math](\nabla{\vec u(ρ,θ)})^t[/math].

  • [math]\nabla{\vec u(ρ,θ)}= \left(\begin{matrix} 0 & (\frac{ρ}{20}+\frac{1}{60ρ}-\frac{ρ^2}{15})senθ & 0 \\ (\frac{2ρ}{15}-\frac{1}{20}+\frac{1}{60ρ^2})senθ & (\frac{ρ^2}{15}-\frac{ρ}{20}-\frac{1}{60ρ})cosθ & 0 \\ 0 & 0 & 0 \end{matrix}\right)[/math].



  • [math](\nabla{\vec u(ρ,θ)})^t = \left(\begin{matrix} 0 & (\frac{2ρ}{15}-\frac{1}{20}+\frac{1}{60ρ^2})senθ & 0 \\ (\frac{ρ}{20}+\frac{1}{60ρ}-\frac{ρ^2}{15})senθ & (\frac{ρ^2}{15}-\frac{ρ}{20}-\frac{1}{60ρ})cosθ & 0 \\ 0 & 0 & 0 \end{matrix}\right)[/math].


Por lo tanto el tensor de deformaciones será igual a:

  • [math]e(\vec u(ρ,θ))= \left(\begin{matrix} 0 & (-\frac{ρ^2}{30}+\frac{11ρ}{120}-\frac{1}{40}+\frac{1}{120ρ}+\frac{1}{120ρ^2})senθ & 0 \\ (-\frac{ρ^2}{30}+\frac{11ρ}{120}-\frac{1}{40}+\frac{1}{120ρ}+\frac{1}{120ρ^2})senθ & (\frac{ρ^2}{15}-\frac{ρ}{20}-\frac{1}{60ρ})cosθ & 0 \\ 0 & 0 & 0 \end{matrix}\right)[/math].



En un medio elástico lineal, isótropo y homogéneo los desplazamientos permiten escribir el tensor de tensiones [math]\sigma(ρ,θ)[/math] a través de la fórmula:

[math]\sigma(ρ,θ) = λ(\nabla · \vec u(ρ,θ))I+2μe(\vec u(ρ,θ)) [/math].


Donde [math]λ[/math] y [math]μ[/math] son conocidos como coeficientes de Lamé que dependen de las propiedades elásticas de cada material. Concretamente:

[math] λ=\frac{Eν}{(1+ν)(1-2ν)}, 2μ=\frac{E}{(1+ν)} [/math],

siendo [math]E[/math] el módulo de Young y [math]ν[/math] el parámetro de Poisson.
El tensor de tensiones quedaría como:

  • [math]\sigma_{ij}(ρ,θ)=\begin{pmatrix}( \frac{ρ}{15}- \frac{1}{20} - \frac{1}{60ρ^2})cosθ& (-\frac{ρ^2}{15} +\frac{11ρ}{60}-\frac{1}{20}+\frac{1}{60ρ}+\frac{1}{60ρ^2})sinθ& 0 \\ (-\frac{ρ^2}{15} +\frac{11ρ}{60}-\frac{1}{20}+\frac{1}{60ρ}+\frac{1}{60ρ^2})sinθ & (\frac{2ρ^2}{15} -\frac{ρ}{30}-\frac{1}{20}-\frac{1}{30ρ}-\frac{1}{60ρ^2})cosθ & 0 \\ 0 & 0 & ( \frac{ρ}{15}- \frac{1}{20} - \frac{1}{60ρ^2})cosθ \end{pmatrix}[/math].


8.2 Cálculo de las tensiones normales


Vamos a comenzar por las tensiones en la dirección normal al plano. De acuerdo a nuestro croquis podemos apreciar que se trata de la proyección de la tensión [math] \vec σ[/math] al eje normal del plano, que tiene dirección:

[math]\vec n[/math].

La proyección no es más que el tensor:

[math] T = \vec n\otimes\vec n · \vec σ[/math]


derecha

8.2.1 Tensiones normales en la dirección que marca el vector [math]\vec g_ρ [/math] y [math] \vec g_z [/math]


Las tensiones normales en la dirección que marca el vector [math]\vec g_ρ [/math] y [math] \vec g_z [/math] que son de la misma magnitud, se obtienen a partir de la expresión:

|[math](σ·\vec g_ρ))\vec g_ρ| = [/math]|[math]( \frac{ρ}{15}- \frac{1}{10} - \frac{1}{60ρ^2})cosθ\vec g_ρ + (-\frac{ρ^2}{15} +\frac{11ρ}{60}-\frac{1}{20}+\frac{1}{60ρ}+\frac{1}{60ρ^2})sinθ\vec g_θ = ( \frac{ρ}{15}- \frac{1}{20} - \frac{1}{60ρ^2}) cosθ[/math]|



8.2.2 Tensiones normales en la dirección que marca el vector [math]\frac{\vec g_θ}{ρ} [/math]


Las tensiones normales en la dirección que marca el vector [math]\frac{\vec g_θ}{ρ} [/math] se obtienen a partir de la expresión:

|[math](σ·\vec g_θ))\vec g_θ| = [/math]|[math] (-\frac{ρ^2}{15} +\frac{11ρ}{60}-\frac{1}{20}+\frac{1}{60ρ}+\frac{1}{60ρ^2})sinθ\vec g_θ + (\frac{2ρ^2}{15} -\frac{ρ}{30}-\frac{1}{20}-\frac{1}{30ρ}-\frac{1}{60ρ^2})cosθ \vec g_θ = (\frac{2ρ^2}{15} -\frac{ρ}{30}-\frac{1}{20}-\frac{1}{30ρ}-\frac{1}{60ρ^2})cosθ[/math]

Para el cálculo de las tensiones normales hemos utilizado el tensor de tensiones calculado anteriormente [math]\sigma(ρ,θ)[/math] respecto de la base {[math]\vec g_ρ,\frac{\vec g_θ}{ρ},\vec g_z[/math]}.

8.3 Visualización de las tensiones normales


La primera imagen muestra el módulo de las tensiones que marca el vector [math]\vec g_ρ [/math] y [math]\vec g_z [/math].
La segunda imagen muestra el módulo de las tensiones que marca el vector [math]\frac{\vec g_θ}{ρ} [/math].



Las tres tensiones son debidas a la deformación de la placa.

centro


centro


Como podemos ver, las tensiones normales en la dirección [math] \frac{\vec g_θ}{ρ} [/math] son mayores que en la dirección de [math]\vec g_ρ [/math] y [math]\vec g_z [/math]. 

%Inicializamos las variables
rho = [1:0.1:2];
h = round(pi/0.1);
theta= linspace(0,pi,h);

%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta); 

subplot(2,2,1)

X=RHO.*cos(THETA);  
Y=RHO.*sin(THETA);
t=abs((RHO./15-1/20-1./(60*RHO.^2)).*cos(THETA));
surf(X,Y,t)
colorbar;
axis([-3,3,-1,3]); 
view(2);
title('Tensiones normales en la dirección g sub rho y g sub z')
subplot(2,2,2)
mesh(X,Y,t)
colorbar;
axis([-3,3,-1,3]); 
view([0,0]);
title('Tensiones normales en la dirección g sub rho y g sub z')


subplot(2,2,3)
tt=abs((2*RHO.^2/15-RHO./30-1/20-1/30.*RHO-1/60.*RHO.^2).*cos(THETA));
surf(X,Y,tt)
colorbar;
axis([-3,3,-1,3]); 
view(2);
title('Tensiones normales en la dirección g sub theta/rho')

subplot(2,2,4)
mesh(X,Y,tt)
colorbar;
axis([-3,3,-1,3]); 
view([0,0]);

title('Tensiones normales en la dirección g sub theta/rho')


9 Tensiones Tangenciales

derecha


Las tensiones en la dirección tangencial, de acuerdo a nuestro croquis, es la proyección de la tensión [math]\vec σ [/math] al plano, que tiene dirección:

[math]\vec t[/math].

La proyección no es más que el tensor:

[math] T = (\vec σ - \vec n\otimes\vec n [/math]·σ)[math]\vec t[/math]


9.1 Cálculo del módulo de las tensiones tangenciales


9.1.1 Módulo de la tensiones tangenciales con respecto a los planos ortogonales a {[math]\vec g_ρ,\frac{\vec g_θ}{ρ}[/math]}


  • |[math]σ·\vec g_ρ-(\vec g_ρ·(σ·\vec g_ρ))\vec g_ρ|=( \frac{ρ}{15}- \frac{1}{10} - \frac{1}{60ρ^2})cosθ\vec g_ρ + (-\frac{ρ^2}{15} +\frac{11ρ}{60}-\frac{1}{20}+\frac{1}{60ρ}+\frac{1}{60ρ^2})sinθ\vec g_θ - ( \frac{ρ}{15}- \frac{1}{10} - \frac{1}{60ρ^2})cosθ\vec g_ρ = (-\frac{ρ^2}{15}+\frac{11ρ}{60}-\frac{1}{20}+\frac{1}{60ρ}+\frac{1}{60ρ^2})sinθ [/math]




  • |[math]σ·\vec g_ρ-(\vec g_ρ·(σ·\vec g_ρ))\vec g_ρ| = (-\frac{ρ^2}{15}+\frac{11ρ}{60}-\frac{1}{20}+\frac{1}{60ρ}+\frac{1}{60ρ^2})sinθ [/math]

9.1.2 Módulo de la tensiones tangenciales con respecto a los planos ortogonales a {[math]\vec g_z[/math]}


  • |[math]σ·\vec g_z-(\vec g_z·(σ·\vec g_z))\vec g_z[/math]| = 0

9.2 Visualización de las tensiones tangenciales


centro 

%Inicializamos las variables
rho = [1:0.1:2];
h = round(pi/0.1);
theta= linspace(0,pi,h);


%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta); 
X=RHO.*cos(THETA);  
Y=RHO.*sin(THETA);
subplot(2,1,1)
t=abs(-RHO.^2/15+11*RHO./60-1/20+1./(60*RHO)+1./(60*RHO.^2)).*sin(THETA);
surf(X,Y,t)
colorbar;
axis([-3,3,-1,3]); 
view(2);
title('Tensiones tangenciales de g sub rho respecto al plano g sub rho y g sub theta')

subplot(2,1,2)
surf(X,Y,t)
colorbar;
axis([-3,3,-1,3]); 
axis vis3d
title('Tensiones tangenciales de g sub rho respecto al plano g sub rho y g sub theta')


10 Tensión de Von Mises


10.1 Cálculo de la tensión de Von Mises


La tensión de Von Misses se define por la fórmula:

[math]\ σ_{VM}= \sqrt {\frac{(σ_1-σ_2)^2+(σ_2-σ_3)^2+(σ_3-σ_1)^2}{2}}[/math],

donde [math]\ σ_1,σ_2,σ_3, [/math] son los autovalores del tensor de tensiones [math]σ(ρ,θ)[/math], cuya matriz de componentes respecto de la base {[math]\vec g_ρ,\frac{\vec g_θ}{ρ},\vec g_z[/math]} es:

[math]\sigma_{ij}(ρ,θ)=\begin{pmatrix}( \frac{ρ}{15}- \frac{1}{20} - \frac{1}{60ρ^2})cosθ& (-\frac{ρ^2}{15} +\frac{11ρ}{60}-\frac{1}{20}+\frac{1}{60ρ}+\frac{1}{60ρ^2})sinθ& 0 \\ (-\frac{ρ^2}{15} +\frac{11ρ}{60}-\frac{1}{20}+\frac{1}{60ρ}+\frac{1}{60ρ^2})sinθ & (\frac{2ρ^2}{15} -\frac{ρ}{30}-\frac{1}{20}-\frac{1}{30ρ}-\frac{1}{60ρ^2})cosθ & 0 \\ 0 & 0 & ( \frac{ρ}{15}- \frac{1}{20} - \frac{1}{60ρ^2})cosθ \end{pmatrix}[/math].


Se trata de una magnitud escalar que se suele usar como indicador para saber cuando un material inicia un comportamiento plástico y no elástico puro.

10.2 Visualización de la tensión de Von Mises


La tensión de Von Mises representa la máxima tensión que puede soportar la placa antes de iniciar un comportamiento plástico.
En el ejercicio esta tensión máxima tiene un valor de 0,31667, y se alcanza en la zona externa inferior de la placa, concretamente en [math]x=-2[/math] y [math]x=2[/math].

centro 



%tension de von misses
rho = 1:0.1:2;
h = round(pi/0.1);
theta = linspace(0,pi,h);
[RR,TT]=meshgrid(rho,theta);
xx=RR.*cos(TT);
yy=RR.*sin(TT);
M11 = @(R,T) (R/15-1/20-1/(60*R^2))*cos(T);
M12 = @(R,T) (-R^2/15+11*R/60-1/20+1/(60*R)+1/(60*R^2))*sin(T);
M22 = @(R,T) (2*R^2/15-R/30-1/20-1/(30*R)-1/(60*R^2))*cos(T);
%matriz de tensiones
 sig = [];
 vm = zeros(length(theta), length(rho));
 
for i = 1:length(theta)
    for j = 1:length(rho)
      sig(1,1) = M11(RR(i,j),TT(i,j));
      sig(1,2) = M12(RR(i,j),TT(i,j));
      sig(1,3) = 0;
      sig(2,1) = M12(RR(i,j),TT(i,j));
      sig(2,2) = M22(RR(i,j),TT(i,j));
      sig(2,3) = 0;
      sig(3,1) = 0;
      sig(3,2) = 0;
      sig(3,3) = M11(RR(i,j),TT(i,j));
       
        [w,d] = eig(sig);
        
         vm(i,j) = sqrt(((d(1,1)-d(2,2))^2+(d(2,2)-d(3,3))^2+(d(3,3)-d(1,1))^2)/2);
        
    end
end

subplot(1,3,2);
surf(xx,yy,vm);
axis([-3,3,-1,3])
axis vis3d
colorbar;
xlabel('Eje x')
ylabel('Eje y')
zlabel('Eje z')
title('Tensión de Von Mises');
 
subplot(1,3,3)
hold on
 
m = vm(1:size(vm,1), 11);
e = xx(1:length(xx),11);
 
surf(xx,yy,vm);
view(0,0)
axis([-3,3,-1,3])
colorbar;
title('XOZ');

maxvm = max(m);

for k = 1:length(m)
   if m(k) == maxvm
       plot3(e(k),0,maxvm,'xr','markersize',10)
   end
 end
 
 xlabel('Eje x')
 zlabel('Eje z')
 txt= ['valor max = ' num2str(maxvm)]; 
 text(-1.85,maxvm+0.03,0.33,txt)
 hold off
 
 subplot(1,3,1)

 
 
surf(xx,yy,vm);
view(2)
axis([-3,3,-1,3])
colorbar;
title('Tensión de Von Mises');



11 Círculo de Mohr


El círculo de Mohr de tensiones es un círculo dibujado en el plano σ-τ en el que cada punto de su circunferencia representa las tensiones normales y cortantes en un plano AB con una inclinación cualquiera.

El círculo de Mohr se utiliza como recurso gráfico para el análisis de las tensiones en estados tensionales biaxiales.

Gracias a esto, podemos definir cuales son las direcciones de las tensiones principales.
En este caso, calculamos las tensiones en el punto (U(ρ,θ)[2,π/4]).
Mohr2d.png derecha



































11.1 Cálculo del círculo a partir de las tensiones

%Mohr en el punto (U(ρ,θ)[2,π/4])
clc

       %DATOS(tensiones)%pasando de cartesianas a cilindricas%metiendo el punto dentro de las formulas obtenidas anteriormente%%

RHO=2;
THETA=pi/4;
tensionx=abs((RHO./15-1/20-1./(60*RHO.^2)).*cos(THETA));
tensiony=abs((2*RHO.^2/15-RHO./30-1/20-1/30.*RHO-1/60.*RHO.^2).*cos(THETA));
ttangen=abs(-RHO.^2/15+11*RHO./60-1/20+1./(60*RHO)+1./(60*RHO.^2)).*sin(THETA);

%centro
sigma = (tensionx+tensiony)/2;

%radio[R]
R = (((tensionx-tensiony)/2)^2+ttangen^2)^0.5;
tau1 = +R;
tau2 = -R;
%tension principal
sigma1 = sigma + R;
sigma2 = sigma - R;

%angulo phi[phi_p]
phi_p = 0.5*atan(2*ttangen/((tensionx-tensiony)))*180/(pi);
%angulo complementario a phi_p
phi_s = 0.5*atan(-((tensionx-tensiony)/(2*ttangen)))*180/(pi);


%datos circulo
x0 = sigma;         % x centro circulo
y0 = 0;             % y centro circulo
N = 500;            % numero de puntos

% 
Theta = linspace(0,2*pi,N);
r = R * ones(1,N);

%pasamos a cartesianas
[x, y] = pol2cart(Theta, r);

%centrado
x = x + x0;         
y = y + y0;

%%%%%%%%%%%%%%%%%%%%%%Linea phi%%%%%%%%%%%%%%%%%%%%%%%%%%%
x1 = sigma2:0.25:sigma1+0.25;
y1 =(ttangen/(R-sigma))*x1 - (ttangen/(R-sigma))*sigma;

x2= sigma2:0.15:sigma1;
y2=0;
x3=sigma;
y3= -R:0.15:R;

%dibujo
plot(x,y,'b-',x1,y1,'r-', x2,y2,'k.', x3 ,y3, 'k.');   
title('CÍRCULO DE MOHR')
xlabel('\sigma \rightarrow ')        
ylabel('\tau^{cw} \rightarrow ')
txt4 = ['\leftarrow RADIO = ' num2str(R)];
text(tau1+x0-0.05,0.5,txt4)
grid on
axis equal        
%
txtC='\sigma';
text(sigma,0,txtC)

%Sigma 1

txt1 = '\leftarrow \sigma_{1}';
text(sigma1,0,txt1)

%Sigma 2

txt2 = '\leftarrow \sigma_{2}';
text(sigma2,0,txt2)

%%%%%punto de interseccion circulo%%%%%
x3 = R;
y3 = ttangen;
txt3 = 'R, \tau_{xy}';
text(x3,y3,txt3)


11.2 Interpretación del resultado

12 Cálculo de la masa total de la placa


Si la placa tiene como densidad [math]d(x,y)=1+(xy)ln(1+x+y^2)[/math], para calcular la masa de dicha placa se procede de la siguiente manera:

  • Se parametriza la superficie en coordenadas cartesianas:

[math]\ x(u,v)=ucosv[/math], [math]\ y(u,v)=usenv[/math], [math]\ z(u,v)=0[/math], [math]\ u∈[1,2] [/math] y [math]\ v∈[0,π][/math].

  • Se calculan las coordenadas cilíndricas de la placa:

[math]\ ρ(u,v)= \sqrt {(x(u,v))^2+(y(u,v))^2}=u [/math], [math]\ θ(u,v)=arctg(\frac{y(u,v)}{x(u,v)})=v[/math], [math]\ z(u,v)=z(u,v)=0[/math], [math]\ u∈[2,3] [/math] y [math]\ v∈[0,π][/math]. De manera que la parametrización de la placa en coordenadas cilíndricas será igual a: [math]\vec r(u,v)=u\vec g_ρ+v\vec g_θ\ [/math].

  • Se calculan los vectores tangentes a la placa:

[math]\vec r_u(u,v)=\frac{\partial \vec r(u,v)}{\partial u}=\vec g_ρ[/math], [math]\vec r_v(u,v)=\frac{\partial \vec r(u,v)}{\partial v}=\vec g_θ[/math].

  • Se calcula el vector normal a la placa:

[math]\vec n(u,v)=\vec r_u(u,v)×\vec r_v(u,v)=ρ(u,v)\vec g_z=u\vec g_z[/math], cuyo módulo es: [math]|\vec n(u,v)|=|\vec r_u(u,v)×\vec r_v(u,v)|=u[/math].


La masa de la placa vendrá dada por la integral:[math]\int_{S} ddS=\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] [math]=\int_{0}^{π}\int_{1}^{2}u(1+(x_1(u,v)x_2(u,v))ln(1+x_1(u,v)+(x_2)^2)(u,v))dudv=9,3273 [/math] u.


Como en este caso, la primitiva es difícil de calcular, se ha utilizado, como método numérico, el siguiente programa para aproximar la integral. 

%Calculo de la masa con un programa de matlab. Queremos hacernos una idea
%del resultado que obtendremos. Aqui he utilizado que sin
%(2T)=2sin(T)*cos(T)
f=@(R,T) R.*(1+R.^2.*abs(sin(2*T))/2.*log(1+R.*abs(cos(T))+(R.*sin(T)).^2));
res=integral2(f,1,2,0,pi)


%Ahora la vamos a calcular utilizando la aproximación numérica mas
%sencilla. Hago una malla en [1,2]x[0,pi]
%obtendré hxk rectángulo en esa malla. Creo una matriz A de ceros y en cada
%elemento de esa matriz guardo el valor del volumen del paralelepípedo
%formado por el rectángulo y el valor del integrando en uno de los
%extremos del rectángulo. Por fin sumo todos los elementos de la matriz y
%obtengo un valor aproximado de la integral
%Calculo de la masa
h=1/100; k=pi/300;
r=1:h:2;
t=0:k:pi;
A=zeros(length(r)-1,length(t)-1);
for i=2:length(r)
    for j=2:length(t)
        A(i,j)=h*k*r(i)*(1+r(i)^2*abs(sin(2*t(j)))/2*log(1+r(i)*abs(cos(t(j)))+(r(i)*sin(t(i)))^2));
    end
end
masa=sum(sum(A))



Con nuestra aproximación la masa sería M=9,2185 u.

Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Visualización_de_campos_escalares_y_vectoriales_en_elasticidad_(Grupo_4B)&oldid=45969»