Visualización de campos escalares y vectoriales en elasticidad. (Grupo 5)

De MateWiki
Saltar a: navegación, buscar
Trabajo realizado por estudiantes
Título Visualización de campos escalares y vectoriales en elasticidad. (Grupo 5)
Asignatura Teoría de Campos
Curso 2023-24
Autores Laura León de Hoz
Julia García Pérez
Álvaro Sedano Hueso
Eduardo Serrano Icarán
Jaime Navalón Gómez
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

Consideramos 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|/2\). En esta placa tendremos definidas dos cantidades físicas:

  • La temperatura [math] T(x,y) [/math], que en nuestro caso depende tanto de la variable espacial y como de la variable x.
    [math]T(x,y)=cos((y-3)^2+x)[/math].
  • Los desplazamientos [math]\overrightarrow { u } (x,y)[/math] producidos por la acción de una fuerza determinada, tal que:

[math]\vec u(ρ,θ)=\frac{1}{2}e^{ρ-1} sin(2θ-\frac{π}{2})[/math]\(\vec e_θ\).
A continuación, estudiaremos el efecto que producen ambas cantidades físicas sobre la placa y las tensiones que éstas generan.

1 Mallado de los puntos interiores del sólido

Las condiciones que debe cumplir el anillo representado son:

  • El radio debe estar contenido entre 1 y 2
  • El anillo está comprendido en el plano \(y≥|x|/2\)
  • Los ejes en el cuadrado son: \((x, y) ∈ [-3; 3] × [−1; 3]\)
  • El paso de muestreo empleado es [math]h=\frac{2}{10}[/math] en las variables x e y que permitirá representar exactamente la placa en el espacio 2-D.

En coordenadas cilíndricas obtenemos: [math](ρ,θ)[/math] ∈ [1,2]×[arctg(1/2),[math]\pi[/math]-arctg(1/2)] 

%Apartado 1

%Mallado
h = 2/10;
rho = 1:h:2;
thetha = linspace(atan(1/2),pi-atan(1/2),pi/h);
[R,T]=meshgrid(rho,thetha);
x = R.*cos(T);
y = R.*sin(T);
mesh(x,y,0.*x);

%Visualicizacion del mallado
axis([-3,3,-1,3])
xlabel('Eje x')
ylabel('Eje y')
title ('Mallado')
axis equal
view(2)
Mallado de la placa

2 Curvas de nivel de la temperatura y gradiente de temperatura

2.1 Curvas de nivel de la temperatura

La distribución de temperaturas en cada punto del sólido viene dada por el campo escalar: [math]T(x,y)=cos((y-3)^2+x)[/math] 

%Apartado 2.1
h = 2/10;
rho = 1:h:2;
thetha = linspace(atan(1/2),pi-atan(1/2),pi/h);
[R,T] = meshgrid(rho,thetha);

%Temperatura
x = R.*cos(T);  %Cambio a coordenadas cilíndricas
y = R.*sin(T);
T = cos((y-3).^2+x);  %Función de temperaturas

%Representación el campo escalar de temperaturas
hold on
subplot(1,2,1)
surf(x,y,T) 
view(2)
axis ([-3,3,-1,3])
xlabel('Eje X')
ylabel('Eje Y')
colorbar
subplot(1,2,2)
contour(x,y,T,4/h)
title('Curvas de Nivel')
colorbar
axis ([-3,3,-1,3])
xlabel('Eje X')
ylabel('Eje Y')
hold off

Maximo = max(max(T))  %0.9999
Campo escalar de temperaturas y curvas de nivel

Punto objeto







Obtenemos una temperatura máxima de 0.999912 ºC

2.2 Gradiente de temperatura

Al tener las curvas de nivel de nuestro sólido y la función en la que varía la temperatura podemos obtener el gradiente de la misma, [math]∇T[/math].

La fórmula general del cálculo de [math]∇T[/math] en coordenadas cartesianas:
[math]\nabla T=\frac { \partial T(x,y) }{ \partial x } \overrightarrow { i } \quad +\frac { \partial T(x,y) }{ \partial y } \overrightarrow { j } \quad +\frac { \partial T(x,y) }{ \partial z } \overrightarrow { k } [/math].
Particularizándolo a nuestra temperatura: [math] ∇T=-sin((y-3)^2+x)\vec i-(2y-6)·sin((y-3)^2+x)\vec j [/math]

%Apartado 2.2

%Mallado
%Gradiente
h = 2/10;  %Paso de Muestreo
rho = 1:h:2;  %Asignación de las variables rho y theta
thetha = linspace(atan(1/2),pi-atan(1/2),pi/h);
[R,T] = meshgrid(rho,thetha);  %Creación del mallado
x= R.*cos(T);  %Cambio a coordenadas cilíndricas
y= R.*sin(T);

%Función de temperatura
T = cos((y-3).^2+x); 
Tx = -sin((y-3).^2 + x);
Ty = -sin((y-3).^2 + x).*(2*y-6)

%Representación del gradiente en 3D
figure(1)
subplot(1,2,1)  %Graficación en 3D
hold on
surf(x,y,T)
quiver3(x,y,T,Tx,Ty,0*T) 
view([75 25])
axis equal
xlabel('Eje X')
ylabel('Eje Y')
zlabel('Eje Z')
title('Gradiente 3D')
hold off

%Representacion del gradiente en 2D
subplot(1,2,2)
mesh(x,y,0*x);
contour(x,y,T);
view(2);
axis equal
axis([-3,3,-1,3]);
colorbar;
hold on
quiver(x,y,Tx,Ty);
axis equal
axis([-3,3,-1,3]);
view(2);
title('Gradiente en 2D');
hold off
Gradiente en 2D y 3D

Como podemos observar el [math]∇T[/math] es ortogonal a las curvas de nivel del anillo.

3 Ley de Fourier y energía calorífica

La Ley de Fourier dice que la energía calorífica [math]\vec Q~[/math] viaja de acuerdo con la fórmula:

[math]\vec Q=-k∇T[/math]

Tal que k es la constante de conductividad térmica de la placa que supondremos k=1. Y quedaría:

[math]\vec Q=-∇T[/math]
%Apartado 3

%Mallado
%Gradiente
h = 2/10;  %Paso de Muestreo
rho = 1:h:2;  %Asignación de las variables rho y theta
thetha = linspace(atan(1/2),pi-atan(1/2),pi/h);
[R,T] = meshgrid(rho,thetha);  %Creación del mallado
x= R.*cos(T);  %Cambio a coordenadas cilíndricas
y= R.*sin(T);

%Función de temperatura
T = cos((y-3).^2+x); 
Tx = sin((y-3).^2 + x);
Ty = sin((y-3).^2 + x).*(2*y-6)

%Representación del gradiente en 3D
figure(1)
subplot(1,2,1)  %Graficación en 3D
hold on
surf(x,y,T)
quiver3(x,y,T,Tx,Ty,0*T) 
view([75 25])
axis equal
xlabel('Eje X')
ylabel('Eje Y')
zlabel('Eje Z')
title('Energía calorífica 3D')
hold off

%Representacion del gradiente en 2D
subplot(1,2,2)
mesh(x,y,0*x);
contour(x,y,T);
view(2);
axis equal
axis([-3,3,-1,3]);
colorbar;
hold on
quiver(x,y,Tx,Ty);
axis equal
axis([-3,3,-1,3]);
view(2);
title('Energía calorífica en 2D');
hold off
Gradiente en 2D y 3D

4 Campo de vectores en el mallado del sólido

Este apartado no es objeto del trabajo al ser dato el campo [math]\vec u[/math].
Observamos como es la gráfica de este campo vectorial [math]\vec u(ρ,θ)=\frac{1}{2}e^{ρ-1} sin(2θ-\frac{π}{2})[/math]\(\vec e_θ\) en matlab: 

%Apartado 4

h = 2/10;
rho = 1:h:2;
thetha = linspace(atan(1/2),pi-atan(1/2),pi/h);
[R,T] = meshgrid(rho,thetha);
x = R.*cos(T); %Cambio a coordenadas cilíndricas 
y = R.*sin(T);
hold on
mesh(x,y,0.*x);
A = (1/2) .* exp(R-1) .* sin(2.*T-(pi/2)) .* -sin(T); %Campo vectorial
B = (1/2) .* exp(R-1) .* sin(2.*T-(pi/2)) .* cos(T);
quiver(x,y,A,B)
axis([-3,3,-1,3]);
title('Campo de vectores')
hold off
Campo vectorial

5 Desplazamiento del sólido debido al campo de vectores

Analizamos el desplazamiento producido por el campo [math]\vec u(ρ,θ)=\frac{1}{2}e^{ρ-1} sin(2θ-\frac{π}{2})[/math]\(\vec e_θ\) en el cuarto de anillo, y comparamos ambas graficas para poder observarlo bien. 

%Apartado 5

%Mallado
h = 2/10;
rho = 1:h:2;
thetha = linspace(atan(1/2),pi-atan(1/2),pi/h);
[R,T] = meshgrid(rho,thetha);
x = R.*cos(T);
y = R.*sin(T);

% Situación inicial
subplot(2,2,1)
surf(x,y,y*0);
axis([-3,3,-1,3]);
view(2)
title('NO DESPLAZADA')

% Situación final
subplot(2,2,2)
A = (1/2) .* exp(R-1) .* sin(2.*T-(pi/2)) .* -sin(T); %Cambio a cartesianas
B = (1/2) .* exp(R-1) .* sin(2.*T-(pi/2)) .* cos(T);
X = x+A;
Y = y+B;
surf(X,Y,Y*0);
axis([-3,3,-1,4]);
view(2)
title('DESPLAZADA')

% Comparación de las 2 situaciones
subplot(2,2,3)
surf(X,Y,Y*0);
hold on
surf(x,y,y*0);
axis([-3,3,-1,4]);
view(2)
title('COMPARACIÓN')
hold off
Desplazamiento debido al campo de vectores

6 Divergencia

Mediante la divergencia estudiamos el cambio de volumen local provocado por el campo de desplazamientos. Tomaremos el campo [math]\vec u(ρ,θ)=\frac{1}{2}e^{ρ-1} sin(2θ-\frac{π}{2})[/math]\(\vec e_θ\).
La divergencia en cilíndricas se expresa tal que:
[math]\nabla\cdot\vec u=\frac{1}{\rho}[\frac{\partial}{\partial{\rho}}({\rho·u_ρ})+\frac{\partial}{\partial{\theta}}(u_θ)] [/math]
Calculamos la divergencia de nuestro campo: [math] \nabla\cdot\vec u= \frac{1}{ρ}e^{ρ-1}·cos(2θ- \frac{π}{2}) [/math] y a continuación la representamos.

%Apartado 6

%Mallado
h = 2/10; 
rho = 1:h:2; 
thetha = linspace(atan(1/2),pi-atan(1/2),pi/h);
[R,T] = meshgrid(rho,thetha); 
x = R.*cos(T); 
y = R.*sin(T);

%Divergencia
div = (exp(R -1) .* cos(2.*T-pi/2)) ./ R;
subplot(1,2,1) 

 %Gráfica 2D
surf(x,y,div)
shading interp
view(2)
axis([-3,3,-1,3])
colorbar
title('Divergencia 2D')
xlabel('Eje X')
ylabel('Eje Y')

 %Gráfica 3D
subplot(1,2,2) 
view(2)
surf(x,y,div)
shading interp
axis([-3,3,-1,3])
colorbar
axis vis3d
title('Divergencia 3D')
xlabel('Eje X')
ylabel('Eje Y')
zlabel('Eje Z')

% Encontrar el máximo, mínimo y nulo de la divergencia
Maximo_div = max(max(div)) %1.3591
Minimo_div = min(min(div)) %-1.3591
Nulo_div = find(div == 0)
Divergencia en 2D y 3D

La divergencia (∇·[math]\overrightarrow {u})[/math] nos permite conocer en que puntos la sección se comprime (tonos fríos) y en cuales se expande (tonos cálidos).
La divergencia mínima es de -1.35906 en el punto (-1.42186,1.40652)
La divergencia máxima es de 1.35906 en el punto (1.42186,1.40652)
La divergencia nula se da en el punto (0,1)

7 Rotacional

El rotacional representa el giro de la superficie cuando aplicamos el campo de desplazamientos [math]\overrightarrow { u }[/math]
En cilíndircas el rotacional se expresa como:
[math]\nabla×\vec u(ρ,θ) = \frac{1}{ρ}\left|\begin{matrix} \vec e_ρ & ρ\vec e_θ & \vec e_z \\ \frac{\partial}{\partial ρ} & \frac{\partial}{\partial θ} & \frac{\partial}{\partial z} \\ \vec v_ρ & \vec ρv_θ & \vec v_z \end{matrix}\right|[/math].
Tomaremos el campo [math]\vec u(ρ,θ)=\frac{1}{2}e^{ρ-1} sin(2θ-\frac{π}{2})[/math]\(\vec e_θ\).
Calculamos el rotacional de este campo, tal que: [math] \nabla×\vec u = [(\frac{1}{2}e^{ρ-1}·sen(2θ-\frac{π}{2}))·(1+ \frac{1}{ρ})] \vec e_z [/math]
y a continuación lo representamos:

%Apartado 7

%Mallado
h = 2/10;
rho = 1:h:2;
thetha = linspace(atan(1/2),pi-atan(1/2),pi/h);
[R,T]=meshgrid(rho,thetha);
x = R.*cos(T);
y = R.*sin(T);

%Rotacional
rot = (0.5 * exp(R -1) .* sin(2.*T - pi/2)) .* (1./R + 1);

%Representación en 2D
subplot(1,2,1)
surf(x,y,rot)
shading interp
axis([-3,3,-1,3])
view(2)
title('Rotacional en 2D') 
xlabel('Eje X')
ylabel('Eje Y')

%Representación en 3D
subplot(1,2,2)
surf(x,y,rot)
shading interp
colorbar
title('Rotacional en 3D')
xlabel('Eje X')
ylabel('Eje Y')
zlabel('Eje Z')
Rotacional en 2D y 3D

El máximo rotacional es de 2.03871 en el punto (0,2)
El mínimo rotacional es de -1.22323 en el punto (1.78885,0.894427)

8 Tensor de deformaciones

Tomamos [math]ε(\vec u)=\frac {\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. En un medio elástico lineal, isótropo y homogéneo los desplazamientos permiten escribir el tensor de tensiones [math]{\sigma }_{ij}[/math] a través de la fórmula [math]σ=λ\nabla·\vec u 1+2με[/math].
Obtenemos un tensor de dreformaciones tal que:
[math] σ = \begin{pmatrix} \frac{1}{ρ}e^{ρ-1}·cos(2θ-\frac{π}{2}) & \frac{-e^{ρ-1}·sen(2θ-\frac{π}{2})}{2ρ} & \ 0 \\ \frac{-e^{ρ-1}·sen(2θ-\frac{π}{2})}{2ρ} & \frac{1}{ρ}e^{ρ-1}·cos(2θ-\frac{π}{2})+\frac{2·e^{ρ-1}·cos(2θ-\frac{π}{2})}{ρ} & \ 0 \\ \ 0 &\ 0 & \frac{1}{ρ}e^{ρ-1}·cos(2θ-\frac{π}{2})) \end{pmatrix}[/math]

%Apartado 8

%Mallado
h = 2/10;
rho = 1:h:2;
thetha = linspace(atan(1/2),pi-atan(1/2),pi/h);
[R,T] = meshgrid(rho,thetha);
x = R.*cos(T);
y = R.*sin(T);
mesh(x,y,0.*x);

%Definimos la matriz
a = (exp(R - 1) .* cos(2.*T - pi./2)) ./R;  %Elemento (1,1) y (3,3) de la matriz
b = (-exp(R - 1) .* sin(2.*T - pi./2)) ./ 2.*R; %Elemento (1,2) y (2,1) de la matriz
c = (2 .*(exp(R - 1) .* cos(2.*T - pi./2))) ./ R;
d = a+c  %Elemento (2,2) de la matriz

%Representación 2D
%Tensión normal en e_rho
subplot(1,3,1)
surf(x,y,a)
axis equal
shading interp
view(2)
colorbar
title('Tensión normal eje e sub rho')

%Tensión normal en e_theta
subplot(1,3,2)
surf(x,y,d)
axis equal
shading interp
view(2)
colorbar
title('Tensión normal eje e sub theta')

%Tensión normal en e_z
subplot(1,3,3)
surf(x,y,a)
axis equal
shading interp
view(2)
colorbar
title('Tensión normal eje e_z')

%Representación 3D
%Tensión normal en e_rho
figure
subplot(1,3,1)
surf(x,y,a)
axis equal
shading interp
colorbar
title('Tensión normal eje e sub rho')

%Tensión normal en e_theta
subplot(1,3,2)
surf(x,y,d)
axis equal
shading interp
colorbar
title('Tensión normal eje e sub theta')

%Tensión normal en e_z
subplot(1,3,3)
surf(x,y,a)
axis equal
shading interp
colorbar
title('Tensión normal eje e_z')
Tensión normal 2D
Tensión normal 3D

9 Tensiones tangenciales

La tensión tangencial respecto al plano ortogonal a [math]\vec{i}[/math] (en t=0), viene dada por la siguiente expresión: [math]\ \left| \sigma ·\vec i -(\vec i·\sigma·\vec i)\vec i \right| [/math]
Así, el tensor de tensiones tangenciales particularizado es: [math]\frac{-e^{ρ-1}·sin(2θ-\frac{π}{2})}{2ρ}[/math] 

%Apartado 9

%Mallado
h = 2/10;
rho = 1:h:2;
thetha = linspace(atan(1/2),pi-atan(1/2),pi/h);
[R,T] = meshgrid(rho,thetha);
x = R.*cos(T);
y = R.*sin(T);
m = mesh(x,y,0.*x);

%Cambio a coordenadas cilíndricas
x = R.*cos(T);
y = R.*sin(T);

%Función de la tensión tangencial
tang = (-exp(R-1).*sin(2.*T-pi/2))./(2.*R);
tantan = 0;

%Eliminamos los ceros
for i=1:15
    for j=1:6
        ii= tang(i,j);
        if ii~=0
            tantan(i,j)=[tang(i,j)];
         end 
    end
end

%Representación en 2D
hold on    
subplot(2,1,1)
surf(x,y,tantan)
shading interp
axis equal 
xlabel('x')
ylabel('y')
zlabel('z')
view(2)
colorbar
title('Tensión tangencial a la dirección del eje e sub rho')

%Representación en 3D
subplot(2,1,2)
surf(x,y,tantan)
shading interp
axis equal
xlabel('x')
ylabel('y')
zlabel('z')
view(3)
colorbar
title('Tensión tangencial a la dirección del eje e sub rho')
hold off
Tensión tangencial

10 Tensión de Von Mises

La tensión de Von Mises se define como: [math]{ \sigma }_{ VM }=\sqrt { \frac { { \left( { \sigma }_{ 1 }-{ \sigma }_{ 2 } \right) }^{ 2 }+{ \left( { \sigma }_{ 2 }-{ \sigma }_{ 3 } \right) }^{ 2 }+{ \left( { \sigma }_{ 3 }-{ \sigma }_{ 1 } \right) }^{ 2 } }{ 2 } } [/math]. Sabemos que [math]{ \sigma }_{ 1 } ,{ \sigma }_{ 2 }, { \sigma }_{ 3 }[/math] son las tensiones principales, conocidas como los autovalores de [math]\sigma[/math].
La tensión de von Mises es una magnitud escalar que nos indica cuando un material inicia un comportamiento plástico. 

%Aprtado 10

%Mallado
h = 2/10;
rho = 1:h:2;
thetha = linspace(atan(1/2),pi-atan(1/2),pi/h);
[R,T] = meshgrid(rho,thetha);
LR = length(rho);
LT = length(thetha);
X = R.*cos(T);  % Cambio a coordenadas cilíndricas
Y = R.*sin(T);

%Componentes de nuestra matriz
a = @(R,T) (1./R).*exp(R-1).*cos(2.*T-(pi/2)); %Elemento (1,1) y (3,3) de la matriz
b = @(R,T) (1/2).*exp(R-1).*sin(2.*T-(pi/2)); %Elemento (1,2) y (2,1) de la matriz
c = @(R,T) 3.*(1./R).*exp(R-1).*cos(2.*T-(pi/2)); %Elemento (2,2) de la matriz
sig=[];
 
%Autovalores
for i=1:LT
    for j=1:LR
        A = a(R(i,j),T(i,j));
        B = b(R(i,j),T(i,j));
        C = c(R(i,j),T(i,j));
        sig = [A,B,0;B,C,0;0,0,A];
        autovalores = eig(sig);
        V = sqrt(((autovalores(1)-autovalores(2))^2+((autovalores(2)-autovalores(3))^2+((autovalores(3)-autovalores(1))^2))*1/2));
        W(i,j) = V;
    end
end

%Tensión de Von Mises
subplot(1,2,1)
surf(X,Y,W);
axis equal
xlabel('x')
ylabel('y')
view(2)
axis([-3,3,-1,3])
title('Tensión de Von Mises')
colorbar
shading interp
subplot(1,2,2)
surf(X,Y,W);
axis equal
xlabel('x')
ylabel('y')
zlabel('z')
view(3)
title('Tensión de Von Mises')
colorbar
shading interp
%Máximo tensión de Von Mises
Maximo = max(max(W))

%El maximo es 3.3291
Tensión de Von Mises
Punto objeto










11 Campo de fuerzas

12 Módulo del desplazamiento transversal

Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Visualización_de_campos_escalares_y_vectoriales_en_elasticidad._(Grupo_5)&oldid=66158»