Revisión del 17:37 10 dic 2021

En este articulo vamos a considerar una placa plana que ocupa medio anillo circular centrado en el origen, esta en el plano [math]y\le 0 [/math] y esta comprendido entre los radios 1 y 2 . También disponemos de una función temperatura tal que : [math]T(x,y)=log((x-3)^2 +2)[/math]

Y un campo dado : [math]\vec u(\rho,\theta)=\frac{\rho-1}{5}\cdot sin(\theta)e_\theta^\rightarrow [/math]

1 Mallado

Realizaremos en Matlab un mallado de los puntos de la placa sobre la que vamos a trabajar. Dibujaremos la curva en los ejes [math](x,y) ∈ [-3,3]\times [-1,3][/math] con un paso de muestro [math]h=1/10[/math]. dado que nuestro solido es medio anillo circular debemos proceder usando coordenadas cilindricas.

Mallado del anillo

h=0.1;                 %se define h para dividir los intervalos
u=1:h:2;               %intervalo rho y teta
v=0:h:pi;  
 
[uu,vv]=meshgrid(u,v);  %malla de u x v
 
xx=uu.*cos(vv);         %parametrización con los puntos de la malla
yy=uu.*sin(vv);
zz=0.*uu;

figure(1)
mesh(xx,yy,0*xx);        %dibujamos la grafica
axis([-3,3,-1,3]);       %región de la grafica

2 Lineas de nivel de la temperatura

La temperatura nos viene dada por el campo [math]T(x,y)=log((x-3)^2 +2)[/math] que viene en coordenadas cartesianas, como nosostros estamos trabajando en coordenadas cilindricas debemos hacer un cambio. O cambiamos nuestras coordenadas a cartesianas o cambiamos el campo a cilindricas resultando en [math]T(\rho,\theta)=log((\rho\cdot cos\theta-3)^2+ 2)[/math]

Nosotros hemos optado por la primera opción. Como se puede apreciar en la gráfica, el punto en el que la temperatura es máxima será (-2,0). Es importante notar que dado que nuetra temperatura solo depende de la variable [math] x [/math] nuestras curvas resultarán lineas rectas verticales.

T=log(((xx-3).^2)+2);  %Campo escalar temperatura
 figure(2)
hold on
contour(xx,yy,T,125)    %contour para las superficies de nivel;
colorbar
hold off

3 Cálculo del gradiente de la Temperatura

REVISAR Calculamos el gradiente del campo (escalar) de la temperatura que tenemos en el apartado anterior, siendo este: [math]T(x,y)=log((x-3)^2 +2)= log(x^2 -6x+11)[/math]

[math]\nabla T=\frac{\partial T}{\partial x}\cdot\vec i+\frac{\partial T}{\partial y}\cdot\vec j = \frac{2x-6}{x^2 +11 -6x} \vec i + 0 \vec j[/math]

T = log(((xx-3).^2)+2);               %Campo escalar temperatura 
GradienteX=(2*xx-6)./(xx.^2-6*xx+11); %derivada de la temperatura respecto a x 
GradienteY= GradienteX.*0;            %derivada de la temperatura respecto a y
figure(3)   
subplot(1,2,1)
hold on
quiver3(xx,yy,T,GradienteX,GradienteY,0*T,0.5); %representación del gradiente en 3 dimensiones 
surf(xx,yy,T);                                  %surf para superficies
contour(xx,yy,T,20);
axis([-3,3,-1,3]) 
hold off

4 Dibujo del campo de Vectores

Dibujamos nuestro campo de vectores en el mallado del solido. Nuestro campo en cilindricas será:

[math]\vec u=\frac{\rho-1}{5}\cdot sin(\theta)e_\theta^\rightarrow [/math]

Campo de vectores

% Campo vectorial u(uu,vv)
ux=-(((uu-1)./5).*sin(vv)).*sin(vv);  
uy=(((uu-1)./5).*sin(vv)).*cos(vv);
figure(4)
hold on     
axis([-3,3,-1,3])
mesh(xx,yy,0*xx) %mallado del solido
quiver(xx,yy,ux,uy,1.5,'b'); %campo vectorial u(uu,vv)
hold off

5 Solido antes y despues del desplazamiento

Aplicando a nuestro solido el campo [math]\vec u [/math] tendrá el desplazamiento ilustrado más abajo. Como podemos observar tiene una leve deformación en el centro, pero esta pasa casi desapercibida.

figure(5)
subplot(1,2,1)
mesh(xx,yy,0*xx)  %mallado del solido
axis([-3,3,-1,3])
view(2)
axis equal
title('Antes del desplazamiento')
xlabel('Eje x')
ylabel('Eje y')
x1=xx+ux;         %escribimos los desplazamientos
x2=yy+uy;
subplot(1,2,2)
mesh(x1,x2,0*x1); %mallado del sólido desplazado
axis([-3,3,-1,3])
view(2)
axis equal
xlabel('Eje x')
ylabel('Eje y')
title('Después del desplazamiento')

6 Representación de la divergencia

Deberemos calcular la divergencia de [math] \nabla\cdot u [/math] y determinar dónde esta es máxima,mínima y nula. Esta divergencia es un cambio de volumen local dado al desplazamiento pero dado que nuestra figura es una placa lo que observaremos será un cambio en su área.

[math]\nabla\cdot u = 1/\rho\cdot[{\frac{\partial }{\partial \rho}\cdot(\rho\vec u_\rho)+\frac{\partial }{\partial\theta}\cdot(\vec u_\theta)+\frac{\partial }{\partial z}\cdot(\rho\vec u_z)]=0+\frac{(\rho -1)cos(\theta)}{\rho 5}+0}[/math]

DIV=((uu-1)./(5.*uu)).*cos(vv);
figure(6)
surf(xx,yy,DIV)  %Representación de la divergencia
axis([-3,3,-1,3])
axis equal
view(2)
colorbar
title('Divergencia')
maxdivergencia=max(max(DIV));
mindivergencia=min(min(DIV));

7 Cálculo de determinante del rotacional

Debemos calcular [math]|\nabla \times u|[/math] siendo de esta forma un campo escalar y ver qué puntos sufren un mayor rotacional

[math]\nabla×\vec u(\rho,\theta) = \frac{1}{\rho} \begin{pmatrix} \vec e_ρ & \rho\vec e_θ & \vec e_z \\ \frac{\partial}{\partial \rho} & \frac{\partial}{\partial \theta} & \frac{\partial}{\partial z} \\ u_ρ & \rho · u_θ & u_z \end{pmatrix} = {\frac{sin \theta}{5 \rho}} \cdot (2 \rho - 1) \vec e_z [/math]

Los puntos que sufren mayor rotacional como podemos ver en la imagen son los situados alrededor de (0,2)

Rotacional

rot=sin(vv).*(2.*uu-1)./(uu.*5);
figure(7)
surf(xx,yy,rot);
axis([-3,3,-1,3])
axis equal
colorbar
xlabel('Eje x')
ylabel('Eje y')
title('Rotacional')

8 Cálculo de las tensiones

En este caso tenemos la parte simétrica del tensor u tal que [math]e(\vec u)=(\nabla u +\nabla u^t)/2[/math] , esta e la usaremos para calcular los desplazamientos , estos viene dados por la siguiente formula: [math]\sigma=\lambda\nabla\cdot u \cdot I +2\mu e(\vec u)[/math] En la que [math] \lambda , \mu[/math] son coeficientes de Lamé iguales a 1.

Como hemos calculado anteriormente [math]\nabla\cdot u=\frac{(\rho -1)cos(\theta)}{\rho 5} [/math]

Dado que estamos en coordenadas cilindricas y queremos calcular el gradiente en un campo vectorial tendremos que usar los simbolos de Christoffel tal que:

Haciendo su traspuesta tenemos :

y finalmente obtenemos la tensión: [math]\sigma=\lambda\nabla\cdot u \cdot I +2\mu e = [/math]

x

8.1 Tensiones normales en el eje \vec e_\rho

Tal que [math] \vec e_\rho \cdot \sigma \cdot \vec e_\rho [/math]

Siendo [math] \vec e_\rho = \begin{pmatrix} 1\\ 0 \\ 0 \end{pmatrix} [/math]

8.2 Tensiones normales en el eje \vec e_\theta =

[math] \vec e_\theta \cdot \sigma \cdot \vec e_\theta [/math] Con [math] \vec e_ \theta = \begin{pmatrix} 0\\ 1 \\ 0 \end{pmatrix} [/math]

8.3 Tensiones normales en el eje \vec e_z =

[math] \vec e_z \cdot \sigma \cdot \vec e_z [/math] Con [math] \vec e_z = \begin{pmatrix} 0\\ 0 \\ 1 \end{pmatrix} [/math]

9 Cálculo de la tensión tangencial a [math]\vec e_\theta[/math]

Tendremos que calcular [math]|\sigma\cdot \vec e_\rho -(\vec e_\rho \cdot\sigma \cdot\vec e_\rho)\vec e_\rho|[/math]

Como en el apartado anterior tomaremos [math] \vec e_\rho = \begin{pmatrix} 1\\ 0 \\ 0 \end{pmatrix} [/math]

10 Representación y cálculo de la tensión de Von Miles

La tensión de Von Mises viene dada por la formula [math]\sigma_VM = \sqrt {\frac{(\sigma_1 - \sigma_2 )^2 + (\sigma_2 - \sigma_3 )^2 + (\sigma_3 - \sigma_1 )^2}{2}}[/math]

Siendo [math] \sigma_1 , \sigma_2 y \sigma3[/math] autovalores de la tensión previamente calculada.

11 Calculo, dibujo e interpretación del campo de fuerzas [math]\vec F[/math]

Queremos sacra la grafica del campo de fuerzas [math]\vec F = - \nabla \cdot \sigma [/math]

De esta forma tomaremos las filas de [math] \sigma[/math] para poder aplicar la divergencia y así obtener el campo vectorial buscado.