Revisión actual del 22:43 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) o en coordenadas cilindricas [math] \rho =2 y \theta = \pi / 2[/math]

Esto ocurre para los puntos definidos por [math]{\frac{sin \theta}{5 \rho}} \cdot (2 \rho - 1) [/math] ya que así el coseno es máximo y el polinomio dependiente de ρ también. Así se puede apreciar en la imagen que cuanto más nos acercamos a esos valores (más amarillo) más rotados han sido los puntos.

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] y multiplicandolo con la matriz Identidad será : [math]\begin{pmatrix}\frac{(\rho -1)cos(\theta)}{\rho 5} & 0 & 0 \\ 0 & \frac{(\rho -1)cos(\theta)}{\rho 5} & 0 \\ 0 & 0 & \frac{(\rho -1)cos(\theta)}{\rho 5} \end{pmatrix} [/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:

Queremos hallar las componentes de la matriz de la derivada covariante:

Haciendo su traspuesta tenemos : [math]\nabla u + \nabla u^t = \begin{pmatrix} 0 & -\frac{\rho -1}{5 \rho} \cdot sin(\theta) & 0 \\\frac{sin(\theta)}{5} & \frac{\rho -1}{5 \rho} \cdot cos(\theta) & 0 \\ 0 & 0 & 0 \end{pmatrix} + \begin{pmatrix}0 & \frac{sin(\theta)}{5} & 0 \\ -\frac{\rho -1}{5 \rho} \cdot sin(\theta) & \frac{\rho -1}{5} \cdot cos(\theta) & 0 \\ 0 & 0 & 0 \end{pmatrix}[/math]

y finalmente obtenemos la tensión: [math]\sigma=\lambda\nabla\cdot u \cdot I +2\mu e = \begin{pmatrix}\frac{\rho -1}{5 \rho} \cdot cos(\theta) & \frac{sin(\theta)}{5 \rho} & 0 \\ \frac{sin(\theta)}{5 \rho} & 3\frac{\rho -1}{5 \rho} \cdot cos(\theta)& 0 \\ 0 & 0 & \frac{\rho -1}{5 \rho} \cdot cos(\theta) \end{pmatrix} [/math]

%A B son componentes de la matriz del tensor de tensiones
A=((uu-1)./(5.*uu)).*cos(vv);
B=(sin(vv)./(5.*uu);

Tgrho=(A);      %el vector (tensor aplicado a e rho) proyectado a la direccion wrho
Tgtht=(A+2*A);  %el vector (tensor aplicado a e tht) proyectado a la direccion wtht
Tgz=(A);        %el vector (tensor aplicado a e z) proyectado a la dirección wz

t1=Tgrho.*cos(vv);
t2=Tgrho.*sin(vv);

8.1 Tensiones normales en el eje [math]\vec e_\rho [/math]

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] Como vemos eso será el componente (1,1) de la matriz solución.

%dibujamos Tgrho
figure (8)
subplot(1,3,1)
hold on
surf(xx,yy,zz);       %me crea la figura
mesh(xx,yy,Tgrho);
quiver(xx,yy,t1,t2)
title('proyección de la tension direccion erho ')
colorbar
hold off

8.2 Tensiones normales en el eje [math]\vec e_\theta [/math] =

[math] \vec e_\theta \cdot \sigma \cdot \vec e_\theta [/math] Siendo este el componente (2,2) de nuetra matriz solución Con [math] \vec e_ \theta = \begin{pmatrix} 0\\ 1 \\ 0 \end{pmatrix} [/math]

%dibujamos Tgtht
subplot(1,3,2)
hold on
surf(xx,yy,zz);%me crea la figura
mesh(xx,yy,Tgtht);
title('proyección de la tension direccion etht')
colorbar
hold off

8.3 Tensiones normales en el eje [math]\vec e_z [/math] =

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

subplot(1,3,3)
hold on
surf(xx,yy,zz);%me crea la figura
mesh(xx,yy,Tgz);
title('proyección de la tension direccion ez')
colorbar
hold off

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] y también la solución a [math] e_\rho \cdot \sigma \cdot e_\rho = \frac{(\rho -1)cos(\theta)}{\rho 5}[/math]

[math]|\sigma\cdot \vec e_\rho -(\vec e_\rho \cdot\sigma \cdot\vec e_\rho)\vec e_\rho|=\begin{pmatrix}\frac{\rho -1}{5 \rho} \cdot cos(\theta) & \frac{sin(\theta)}{5 \rho} & 0 \\ \frac{sin(\theta)}{5 \rho} & 3\frac{\rho -1}{5 \rho} \cdot cos(\theta)& 0 \\ 0 & 0 & \frac{\rho -1}{5 \rho} \cdot cos(\theta) \end{pmatrix} \cdot \begin{pmatrix} 1\\0\\0 \end{pmatrix} - \frac{(\rho -1)cos(\theta)}{\rho 5} \cdot \begin{pmatrix} 1\\0\\0 \end{pmatrix} = \frac{sin(\theta)}{5 \rho} [/math]

Erho=abs (sin(vv)./(5.*uu));
subplot(1,2,1)
hold on
mesh(xx,yy,Erho)
hold off

subplot(1,2,2)
hold on
mesh(xx,yy,Ehro);
hold off

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.

Esta tension nos dirá que parte de la placa tiene un estrés mayor, y por ende por qué punto se producirá la rotura. En nuestro caso podemos ver que esta se producirá en el punto (0,2) es decir justo por la mitad y la parte exterior del semi anillo.

h=0.1; 
u=1:h:2; 
v=0:h:pi;  
 figure(10)
[uu,vv]=meshgrid(u,v);  
 
val=zeros(32,11);
 
for i=1:32*11
   uum=uu(i);
   vvm=vv(i);
   a=((uum-1).*cos(vvm)./(5.*uum));
   b=-(sin(vvm))./5;
   c=(a./(uum)).*(1+1./uum);
   mat=[a,b,0;b,c,0;0,0,a];
   [v,d]=eig(mat); 
   val(i)=sqrt(((d(1,1)-d(2,2))^2+(d(2,2)-d(3,3))^2+(d(3,3)-d(1,1))^2)/2);  
 end
 
xx=uu.*cos(vv);
yy=uu.*sin(vv);
surf(xx,yy,val)
colorbar
view(2)

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

Queremos sacar la gráfica del campo de fuerzas [math]\vec F = - \nabla \cdot \sigma [/math]

[math]\vec F = - \nabla \cdot \begin{pmatrix}\frac{\rho -1}{5 \rho} \cdot cos(\theta) & \frac{sin(\theta)}{5 \rho} & 0 \\ \frac{sin(\theta)}{5 \rho} & 3\frac{\rho -1}{5 \rho} \cdot cos(\theta)& 0 \\ 0 & 0 & \frac{\rho -1}{5 \rho} \cdot cos(\theta)\end{pmatrix} [/math]

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

Resultando en :

[math]F1 = -1/\rho\cdot [\frac{\partial }{\partial \rho}\cdot(\rho\frac{\rho -1}{5 \rho} \cdot cos(\theta))+\frac{\partial }{\partial\theta}\cdot( \frac{sin(\theta)}{5 \rho})+\frac{\partial }{\partial z}\cdot( 0)] = -\frac{2 cos(\theta)}{5\rho} [/math]

[math]F2= -1/\rho\cdot[\frac{\partial }{\partial \rho}\cdot(\rho\frac{sin(\theta)}{5 \rho})+\frac{\partial }{\partial\theta}\cdot(3\frac{\rho -1}{5 \rho} \cdot cos(\theta))+\frac{\partial }{\partial z}\cdot(0)] = 3\frac{\rho -1}{5 \rho^2} \cdot sin(\theta) [/math]

[math]F3= 0[/math]

Como podemos ver la fuerza dependerá tanto de ρ como de θ , en nuestro caso las fuerzas tendrán un sentido interior por los costados mientras que parecen salir por la parte superior.

h=0.1; 
u=1:h:2; 
v=0:h:pi;  
[uu,vv]=meshgrid(u,v)

wx=(-2.*cos(vv))./(5.*uu)
wy=(3.*(uu-1).*sin(vv))./(5.*uu.^2)
figure(12)
hold on     
axis([-3,3,-1,3])
mesh(xx,yy,0*xx) %mallado del solido
quiver(xx,yy,wx,wy,1.5,'k'); %campo vectorial F(uu,vv)
hold off