Visualización de campos escalares y vectoriales sobre semianillo (Grupo 21-A)

Trabajo realizado por estudiantes
Título	Visualización de campos escalares y vectoriales sobre semianillo (Grupo21-A)
Asignatura	Teoría de Campos
Curso	2022-23
Autores	José Serna García, Galo Pastor Brezmes, Javier López Gonzalez y Miguel Millaruelo Frontela
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Representación del semianillo

Para empezar este artículo, vamos a representar el anillo con el que trabajaremos. Para ello emplearemos coordenadas cilíndricas y lo representaremos entre los ejes [-3,3]x[-1,3].

%Primero parametrizaremos la superficie en función de rho y theta.
r=1:0.2:2;
tt=pi/4:pi/40:(3*pi)/4;
%El termino tt empieza en el valor de pi/4 y no en 0 para conseguir la estructura como en la imagen puesta en el enunciado.

%Después vamos a definir el mallado correspondiente. de forma básica para estructurar los datos y cálculos posteriores.
[R,TT]=meshgrid(r,tt);

%Ahora vamos a definir las variables X e Y en coordenadas cilíndricas (en función de R y TT). Ya que las bases de los ejes donde se nos obliga presentar estan 
%en coordenadas cartesianas.
X=R.*cos(TT);
Y=R.*sin(TT);

%Por último, representaremos la superficie deseada entre los ejes [-3,3]x[-1,3].
figure(1);
subplot(1,2,1);
mesh(X,Y,0*X);
view(2);
axis([-3,3,-1,3]);
subplot(1,2,2);
mesh(X,Y,0*X);
view(3);
axis([-3,3,-1,3]);

2 Curvas de nivel del campo de temperaturas

A continuación veremos cómo actúa el campo escalar de temperaturas "T" sobre nuestro anillo y veremos representadas las curvas de nivel del campo sobre él. Como podemos observar en la imagen, a partir de las curvas de nivel se deduce que el punto del anillo en el cual la temperatura es máxima es el punto (0,2), representado con los colores amarillos dentro de la pieza.

%Primero vamos a parametrizar la superficie igual que en el apartado anterior y a crear el mallado.
r=1:0.2:2;
tt=pi/4:pi/40:(3*pi)/4;
[R,TT]=meshgrid(r,tt);

%Definimos las variables X e Y en función de R y TT.
X=R.*cos(TT);
Y=R.*sin(TT);

%Definimos el campo de temperaturas "T" en función de X e Y, definidas anteriormente a partir de las coordenadas parametrizadas de la placa.
T=((X-3).^2)+((10.*(Y-0.5)).^2);

%Representamos las curvas de nivel del campo T sobre la superficie designada.
figure(1);
mesh(X,Y,0*X);
subplot(1,2,1);
surf(X,Y,T);
axis([-3,3,-1,3]);
view(2);
title('Campo de temperaturas');
subplot(1,2,2);
contour(X,Y,T);
axis([-3,3,-1,3]);
view(2);
title('Curvas de nivel');
colorbar;

3 Gradiente del campo de temperaturas

El gradiente de este campo de temperaturas será: ∇T=(2x-6)i+(200y-100)j. Gracias a esto podemos demostrar que la dirección de aumento máximo de la temperatura es perpendicular a las curvas de nivel del campo. Recordamos que el gradiente nos dará un vector, entonces lo que importa es la dirección que lleva y lo que le da el significado fisico.

%Parametrizamos la superficie y creamos el mallado.
r=1:0.2:2;
tt=pi/4:pi/40:(3*pi)/4;
[R,TT]=meshgrid(r,tt);

%Definimos las variables X e Y en función de R y TT.
X=R.*cos(TT);
Y=R.*sin(TT);

%Definimos el campo de temperaturas "T" en función de X e Y. También las matrices Tx y Ty de su gradiente.
T=((X-3).^2)+((10.*(Y-0.5)).^2);
Tx=(2*X)-6;
Ty=(200*Y)-100;

%Representamos el gradiente del campo T sobre las curvas de nivel del campo mismo.
figure(1)
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 del campo T');
hold off

4 Campo de vectores en los puntos del mallado del sólido

El campo de vectores dado será: u(ρ,θ)=(1/2)log(ρ)sin(2θ − π/2) en la dirección de eρ

%intervalos variables
u=1:0.1:2;
v=pi/4:pi/50:3.*pi/4;
%mallado
[U,V]=meshgrid(u,v);
%parametrización
X=U.*cos(V);
Y=U.*sin(V);
FX=log(U).*sin(2.*V-pi/2).*cos(V);
FY=log(U).*sin(2.*V-pi/2).*sin(V);
%visualizacion del dibujo
quiver(X,Y,FX,FY);
axis([-3,3,-1,3]);
view(2)
xlabel('Eje x')
ylabel('Eje y')
title('Campo de vectores en los puntos del mallado del sólido')

5 Sólido previo y posterior a los desplazamientos

Ahora representaremos las deformaciones producidas en nuestra placa debido a la actuación del campo vectorial u(ρ,θ), mostrando la placa antes y después de las deformaciones, así como una tercera representación comparando ambos casos.

centro

%creamos los intervalos de las variables
u=1:0.1:2,
v=pi/4:pi/50:3.*pi/4;  
%mallado
 [U,V] = meshgrid(u,v);
 X = U.*cos(V);
 Y = U.*sin(V);
%función
 FX=log(U).*sin(2.*V-pi/2).*cos(V);
 FY=log(U).*sin(2.*V-pi/2).*sin(V);
%representación gráfica
%antes de la deformación
 subplot(1,3,1)
 surf(X,Y,0*X);
 view(2)
 axis equal
 axis([-3,3,-1,3])
 title('Antes de la deformación');
%después de la deformación
 subplot(1,3,2)
 UX=X+FX;
 UY=Y+FY;
 surf(UX,UY,0*UX);
 view(2)
 axis equal
 axis([-3,3,-1,3])
 title('Después de la deformación');
%comparación
 subplot(1,3,3)
 hold on
 plot(X,Y,'r');
 plot(UX,UY,'b');
 view(2)
 axis equal
 axis([-3,3,-1,3])
 title('Comparación');
 hold off

6 Divergencia máxima, mínima y nula del campo u sobre la placa

La divergencia del campo u sobre la placa nos muestra la diferencia de areas debidas al desplazamientode la misma. La divergencia del vector u dado sería: 1/(2ρ)sin(2θ -pi/2)(log(ρ)+1)

centro

%creamos los intervalos de las variables
 w=30;
 u=linspace(1,2,w);
 v=linspace(pi/4,3.*pi/4,w); 
 %mallado
  [U,V] = meshgrid(u,v);
  X = U.*cos(V);
  Y = U.*sin(V);
 %divergencia
  DIVu=(1/2.*U).*sin(2.*V -pi./2).*(log(U)+1);
 %gráfica
  surf(X,Y,DIVu)
  view(2)
  axis equal
  axis([-3,3,-1,3])
  colorbar
 %divergencias máximas y mínimas
  DIVuMAX=max(DIVu);
  DIVuMIN=min(DIVu);
  title('Divergencia del campo vectorial')

7 Rotacional

A continuación representaremos el rotacional en los diferentes puntos de nuestra placa.

El rotacional del campo será: [math]\triangledown \times \vec{u} \left ( \rho ,\theta \right )=\left ( \frac{d\vec{e_{\theta }}\rho u_{\rho }}{dz} -\frac{d u_{\rho}}{d\theta } \vec{e_{z}} \right ) [/math] y dado que: [math]\frac{d\vec{e_{\theta }}\rho u_{\rho }}{dz} = 0[/math], podemos concluir que: [math]\left | \triangledown \times \vec{u} \left ( \rho ,\theta \right ) \right | = \frac{\log(\rho )cos(2 \theta -\frac{\pi }{2})}{\rho }[/math].

num = 0.2; 
u = 1:num:2; 
v = linspace(pi/4,3*pi/4,6); 

[rho,theta] = meshgrid(u,v);
x = rho.*cos(theta);
y = rho.*sin(theta);

rotacional = abs((cos(2*theta-pi/2)./(rho)).*(log(rho)));

subplot(1,2,1);
surf(x,y,rotacional);
axis([-3,3,-1,3]);view(2);
colorbar;
title('Rotacional');

centro

8 Vectores normales direccionales

El tensor de tensiones [math]\sigma _{ij}[/math] viene dado por: [math]\sigma = \lambda \triangledown \cdot \vec{u} \mathbf{1} + 2\mu \epsilon [/math], siendo [math]\lambda = \mu = 1[/math]

Siendo [math]\epsilon [/math] definido como: [math]\epsilon \left ( \vec{u} \right ) = \frac{\left ( \triangledown \vec{u} + \triangledown \vec{u}^{t}\right )}{2}[/math] y cuyas componentes de la diagonal principal son: [math]\epsilon_{11} = \frac{sen(2\theta -\frac{\pi }{2})}{2\rho }[/math], [math]\epsilon_{22} = \frac{sen(2\theta -\frac{\pi }{2})\cdot log(\rho )}{2}[/math] y [math]\epsilon_{22} = 0 [/math].

La divergencia del campo será: [math]\triangledown \cdot \vec{u} \left ( \rho ,\theta \right ) = \frac{1}{\rho }\cdot \frac{d\rho u_{\rho }}{d\rho } = \frac{sen(2\theta -\frac{\pi }{2})}{2\rho } \cdot \left ( log(\rho ) +1\right )[/math]

La tensión en dirección [math]\vec{\rho} [/math] será: [math]\sigma _{11} = \triangledown \cdot \vec{u} + 2\epsilon _{11} = \frac{sen(2\theta -\frac{\pi }{2})}{2\rho }\cdot\left ( log(\rho ) + 3\right )[/math]

La tensión en dirección [math]\vec{\theta} [/math] será: [math]\sigma _{22} = \triangledown \cdot \vec{u} + 2\epsilon _{22} = \frac{sen(2\theta -\frac{\pi }{2})}{2\rho }\cdot\left ( log(\rho ) + 1\right ) + log(\rho )\cdot sen(2\theta -\frac{\pi }{2})[/math]

La tensión en dirección [math]\vec{k} [/math] será: [math]\sigma _{33} = \triangledown \cdot \vec{u} + 2\epsilon _{33} = \frac{sen(2\theta -\frac{\pi }{2})}{2\rho }\cdot\left ( log(\rho ) + 1\right )[/math]

num = 0.2; 
u = 1:num:2; 
v = linspace(pi/4,3*pi/4,6); 

[rho,theta] = meshgrid(u,v); 
X = rho.*cos(theta);        
Y = rho.*sin(theta);
Z = zeros(size(X)); 
clf

a = (sin(2*theta-pi()/2)./2*rho).*(log(rho)+3);
b = (sin(2*theta-pi()/2)./2*rho).*(log(rho)+1)+log(rho).*sin(2*theta-pi()/2); 
c = (sin(2*theta-pi()/2)./2*rho).*(log(rho)+1);

subplot(1,3,1)
surf(X,Y,a)
axis([-3,3,-1,3]);
view(2)
colorbar
xlabel('x')
ylabel('y')
zlabel('z')
title('Tensión normal direccional Rho')

subplot(1,3,2)
surf(X,Y,b)
axis([-3,3,-1,3]);
view(2)
colorbar
xlabel('x')
ylabel('y')
zlabel('z')
title('Tensión normal direccional Theta')

subplot(1,3,3)
surf(X,Y,c)
axis([-3,3,-1,3]);
view(2)
colorbar
xlabel('x')
ylabel('y')
zlabel('z')
title('Tensión normal direccional k')

centro

9 Tensiones tangenciales respecto a [math]\vec{\rho} [/math]

La tensión tangencial respecto a [math]\vec{\rho} [/math] viene dada por la siguiente expresión:

[math]\left | \sigma \cdot \vec{e_{\rho }} -\left ( \vec{e_{\rho }} \cdot \sigma \cdot \vec{e_{\rho }} \right ) \cdot \vec{e_{\rho }} \right | = \left | \begin{pmatrix} \sigma _{11} \\ \sigma _{21} \\ \sigma _{31} \end{pmatrix} - \sigma _{11} \cdot \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix} \right | = \left | \begin{pmatrix} \sigma _{11} \\ \sigma _{21} \\ \sigma _{31} \end{pmatrix} - \begin{pmatrix} \sigma _{11} \\ 0 \\ 0 \end{pmatrix} \right | = \left | \begin{pmatrix} 0 \\ \sigma _{21} \\ \sigma _{31} \end{pmatrix}\right |[/math]

Y dado que:

[math]\sigma _{21} = \frac{\frac{du_{\theta }}{d\rho } + \frac{du_{\rho }}{d\theta } - u_{\theta }}{2} = 0[/math]

[math]\sigma _{31} = \frac{\frac{du_{z}}{d\rho } + \frac{du_{\rho }}{dz}}{2} = 0[/math]

Se puede concluir que todas las tensiones tangenciales respecto a [math]\vec{\rho} [/math] son nulas.

10 Tension de Von Mises

Tenemos en cuena que la tension de Von Mises se calcula mediante la formula: [math]\sigma _{VM}= \sqrt{\frac{(\sigma_{1}-\sigma_{2})^2 +(\sigma_{2}-\sigma_{3})^2 +(\sigma_{3}-\sigma_{1})^2}{2}}[/math] donde [math]\sigma_{1},\sigma_{2} y \sigma_{3}[/math] representan las tensiones principales como autovalores de la matriz de tensiones. En este caso estas tensiones ya han sido calculadas como:

Dirección [math]\vec{\rho} [/math] es: [math]\sigma _{1} = \triangledown \cdot \vec{u} + 2\epsilon _{11} = \frac{sen(2\theta -\frac{\pi }{2})}{2\rho }\cdot\left ( log(\rho ) + 3\right )[/math]

Dirección [math]\vec{\theta} [/math] es: [math]\sigma _{2} = \triangledown \cdot \vec{u} + 2\epsilon _{22} = \frac{sen(2\theta -\frac{\pi }{2})}{2\rho }\cdot\left ( log(\rho ) + 1\right ) + log(\rho )\cdot sen(2\theta -\frac{\pi }{2})[/math]

Dirección [math]\vec{k} [/math] es: [math]\sigma _{3} = \triangledown \cdot \vec{u} + 2\epsilon _{33} = \frac{sen(2\theta -\frac{\pi }{2})}{2\rho }\cdot\left ( log(\rho ) + 1\right )[/math]

Por lo que ya son los autovalores de la matriz total de tensiones que debería de calcularse de manera anterior. Por lo que se debe hacer es usar la formula de la tension de Von Mises según este siguiente codigo.

paso = 0.2; 
u = 1:paso:2; 
v = linspace(pi/4,3*pi/4,6); 

%Creamos la matriz de paso mediante la parametrizacion de r y tt
[r,tt] = meshgrid(u,v); 
X = r.*cos(tt);        
Y = r.*sin(tt);
Z = zeros(size(X)); 
clf

sig1 = (sin(2*tt-pi()/2)./2*r).*(log(r)+3);
sig2= (sin(2*tt-pi()/2)./2*r).*(log(r)+1)+log(r).*sin(2*tt-pi()/2); 
sig3= (sin(2*tt-pi()/2)./2*r).*(log(r)+1);

%Calculamos la tension de Von Mises usando la formula, teniendo en cuenta que los parametros son vectores.
sigVM=sqrt((sig1-sig2).^2 +(sig2-sig3).^2 +(sig3-sig1).^2/2)

%El termino se usa en la matriz como valor en cada punto conocido de la superficie.
surf(X,Y,sigVM)
axis([-3,3,-1,3]);
view(2)
title("Tension de Von Mises")

Este grafico nos muestra claramente la acción de las distintas tensiones dentro de la superficie dada para una deformación general del sistema de campos.

11 Fuerzas que actúan sobre la placa

Según los términos del problema debemos aproximar las fuerzas del sistema a la formula dada la Ecuación de Lamé, que es:

[math]\vec{F}[/math]=[math]\frac{\sigma^2\vec{u}}{\sigma t^2}-\mu\triangle\vec{u}-(\lambda+\mu)\triangledown(\triangledown\cdot\vec{u})[/math]

Donde encontramos los términos [math]\lambda[/math] y [math]\mu[/math] representan los coeficientes de la ecuación de Lame que presentan los módulos de corte y de compresibilidad, pero se pueden calcular fácilmente usando variaciones infinitesimales dentro del sistema de coordenadas cilíndricas

En el caso nuestro el coeficiente de compresibilidad [math]\lambda[/math] será representado en el eje de [math]\rho[/math] y el coeficiente de corte [math]\mu[/math] en dirección del eje de [math]\phi[/math] y creando el limite para que los diferenciales sean 0

[math]\lambda[/math]=[math]\frac{MM{r}}{\triangle \rho}[/math]=[math]\frac{\triangle \rho}{\triangle \rho}[/math]=1

[math]\mu[/math]=[math]\frac{MM{\phi}}{\triangle \phi}[/math]=[math]\frac{\triangle \phi*\rho}{\triangle \rho}[/math]=[math]\rho[/math]

Recordamos que ambos coeficientes están ya creados en coordenadas cilíndricas.

En este ejercicio podemos no contar en los calculos el termino [math]\frac{\sigma^2\vec{u}}{\sigma t^2}[/math] debido a que la formula del campo [math]\vec{u}[/math] no viene determinada en ningun momento por el tiempo, entonces el valor de este será de 0

Entonces nos queda por último el termino del Laplaciano representado por: [math]\triangle\vec{u}[/math] que aunque se nos recomienda hacer cambios de variable en la ecuacion del campo [math]\vec{u}[/math] podemos saltarnos ello haciendo el Laplaciano directamente desde las coordenadas cilindricas

[math]\triangle\vec{u}[/math]=[math]\frac{1}{\rho}\frac{\sigma\vec{u}}{\sigma\rho}+\frac{\sigma^2\vec{u}}{\sigma\rho^2}+\frac{1}{\rho^2}\frac{\sigma^2\vec{u}}{\sigma\phi^2}[/math]

Por lo que calculando estos terminos conseguiremos el valor de las Fuerzas que actuan sobre la placa.

syms rho phi

mu=rho;
lambda=1;

%funcion del campo
u=(cos(2*phi-pi/2)./(rho)).*(log(rho));

%derivadas parciales del sistema
Urho=diff(u,rho);
Uphi=diff(u,phi);

%Derivadas pariales segundas del sistema
Urho2=diff(Urho,rho);
Uphi2=diff(Uphi,phi);

%Unimos los valores de las diferentes funciones para calcular el laplaciano
LPL=(1/rho)*Urho +Urho2 +(1/rho^2)*Uphi2;

%Con el Laplaciano calculado Calculamos la ecuancion de Lame.

LAME= -mu*LPL-(lambda+mu)*gradient(Urho+Uphi)

%Calcularemos ahora el resto y mostraremos.
num=0.2;
RHO= 1:num:2;
PHI= linspace(pi/4,3*pi/4,6);

[r,P] = meshgrid(RHO,PHI);
x = r.*cos(P);
y = r.*sin(P);

FUERZA= -RHO.*LPL(RHO,PHI)-(1.+RHO).*gradient(Urho(RHO,PHI)+Uphi(RHO,PHI))

Tambien se puede determinar esto calculando previamente las derivadas para depues calcular todo de manera numerica

%Calculamos la parametrizacion de las variables
num=0.2;
RHO= 1:num:2;
PHI= linspace(pi/4,3*pi/4,6);

[r,P] = meshgrid(RHO,PHI);
x = r.*cos(P);
y = r.*sin(P);

mu=RHO;
lambda=1;

%Funciones de las derivadas

r1=sin(2*PHI-pi/2).*0.5.*1./RHO.*log(10);
r2=sin(2*PHI-pi/2).*(-1./RHO.^2).*(1/2.*RHO.*log(10));

P1=(log10(RHO)./2).*2.*cos(2.*PHI-pi/2);
P2=(log10(RHO)./2).*4.*sin(2.*PHI-pi/2);

%Laplaciano del sistema

LPL=r2+(1./RHO).*r1 + (1./RHO.^2).*P2;

LAME=-mu.*LPL-(lambda+mu).*gradient(r1+P1);