Visualización de campos escalares y vectoriales en elasticidad. Grupo A-15

Consideramos una placa plana que ocupa el anillo circular centrado en el origen y comprendido entre las circunferencias de radios 2 y 1. En ella vamos a suponer que tenemos definidas dos cantidades físicas: la temperatura \(T(x,y)\), que depende de las dos variables espaciales \((x,y)\), y los desplazamientos \(\vec u\) producidos por la acción de una fuerza determinada. A partir de esto, estudiaremos el efecto de la temperatura dada por \(T(x,y)\) para después plantear la acción de las vibraciones \(\vec u(\rho,\theta)\) y las tensiones que éstas producen.

1 Visualización de la placa

Anillo inicial

Consideramos una placa plana con forma de anillo, de radio interno uno y de radio externo dos.

h=0.1;                         % Paso de muestreo
rho=1:h:2;                     % Definición del intervalo de rho [1,2]
theta=0:h:2*pi+h;              % Definición del intervalo de theta [0,2*pi]
[rr,tt]=meshgrid(rho,theta);   % Definición de la placa en coordenadas polares
xx=rr.*cos(tt);                % Parametrización según los ejes a xx y yy
yy=rr.*sin(tt);
mesh(xx,yy,0*xx)               % Visualización de la placa
axis([-3,3,-3,3])              % Elección de los ejes

2 Distribución de temperaturas del sólido

La distribución de temperaturas en cada punto del sólido viene dada por el campo escalar:

                                            \(T(x,y)=e^{-y}\)

Esta variación de temperaturas vendrá representada por curvas de nivel formadas por puntos que se encuentran a una misma temperatura. Considerando como foco de calor el origen de coordenadas y siendo nuestra distribución de temperaturas independiente de \(x\); observamos gráficamente que, cuanto mayor sea \(y\), menor será la temperatura.

Distribución de temperaturas en la placa

h=0.1;                         % Paso de muestreo
rho=1:h:2;                     % Definición del intervalo de rho [1,2]
theta=0:h:2*pi+h;              % Definición del intervalo de theta [0,2*pi]
[rr,tt]=meshgrid(rho,theta);   % Definición de la placa en coordenadas polares
xx=rr.*cos(tt);                % Parametrización
yy=rr.*sin(tt);
T=exp(-yy);                    % Expresión del campo escalar de temperaturas
surf(xx,yy,T)                  % Visualización de la malla en función de la temperatura en 2D
view(2)
axis([-3,3,-3,3])              % Región de visualización

Por tanto, de la visualización bidimensional del campo de temperaturas, resulta que el punto del sólido donde \(y=-2\) será el punto donde la temperatura es máxima.

3 Estudio del gradiente de temperaturas \(\nabla T\)

Al calcular el gradiente del campo escalar \(T(x,y)=e^{-y}\) , obtenemos un campo vectorial \(\nabla T\) expresado como:

                                                 \(\nabla T\)=\(- e^{-y}\vec j\)

Este campo vectorial gradiente nos indica la dirección en la cual la temperatura en cada punto del sólido aumenta más rápido y su módulo nos indica la velocidad con la que aumenta la temperatura en dicha dirección. Según esto, podemos ver gráficamente que a medida que nos acercamos a temperaturas más elevadas, el módulo de \(\nabla T\) en cada punto va siendo mayor.

Campo vectorial gradiente \(\nabla T\)

h=0.1;                         % Paso de muestreo
rho=1:h:2;                     % Definición del intervalo de rho [1,2]
theta=0:h:2*pi+h;              % Definición del intervalo de theta  [0,2*pi]
[rr,tt]=meshgrid(rho,theta);   % Definición de la placa en coordenadas polares
xx=rr.*cos(tt);                % Parametrización
yy=rr.*sin(tt);
T=exp(-yy);                    % Expresión del campo escalar de temperaturas
gradx=0*xx; 	               % Derivada parcial en x
grady=-exp(-yy);               % Derivada parcial en y
hold on
surf(xx,yy,T)                  % Visualización de la malla en función de la temperatura
axis([-3,3,-3,3])
quiver(xx,yy,gradx,grady)      % Visualización del campo vectorial gradiente
axis([-3,3,-3,3]) 
hold off

Las direcciones en las que el campo de temperaturas varía más rápidamente coincidirán con el vector que pase de forma más rápida de una curva de nivel a la siguiente, o lo que es lo mismo, el gradiente será perpendicular a las curvas de nivel.

Campo vectorial gradiente ortogonal a las curvas de nivel

h=0.1;                         % Paso de muestreo
rho=1:h:2;                     % Definición del intervalo de rho [1,2]
theta=0:h:2*pi+h;              % Definición del intervalo de theta  [0,2*pi]
[rr,tt]=meshgrid(rho,theta);   % Definición de la placa en coordenadas polares
xx=rr.*cos(tt);                % Parametrización
yy=rr.*sin(tt);
T=exp(-yy);                    % Expresión del campo escalar de temperaturas
gradx=0*xx; 	               % Derivada parcial en x
grady=-exp(-yy);               % Derivada parcial en y
hold on
quiver(xx,yy,gradx,grady)      % Visualización del campo vectorial gradiente
axis([-3,3,-3,3]) 
contour(xx,yy,T,30)            % Visualización de las curvas de nivel
axis([-3,3,-3,3])
hold off

4 Efectos del campo vectorial \(\vec u\) sobre la placa

A continuación, vamos a estudiar el efecto de aplicar sobre la placa un campo vectorial:

                                         \(\vec u(\rho,\theta)=(1-\rho)^2 \vec g_{\theta}\)

interpretado como una fuerza que producirá la acción de una serie de vibraciones sobre el sólido y, consecuentemente, unos desplazamientos.

4.1 Campo vectorial de desplazamientos

La deformación que sufre la placa viene definida por la función \(\vec u(\rho,\theta)\), que nos indica el desplazamiento que sufren los puntos del solido debido a la vibración producida por la fuerza aplicada:

                                         \(\vec u(\rho,\theta)= (1-\rho)^2 \vec g_\theta\)

El campo de desplazamientos de los puntos de la placa viene definido en la siguiente imagen:

Visualización del campo vectorial sobre el año

figure(4)
ua=1-uu;             
ux=ua.^2.*(-yy);        % Definicion de los desplazamientos horizontales
uy=ua.^2.*xx;           % Definicion de los desplazamientos verticales
quiver(xx,yy,ux,uy)     
axis([-3,3,-3,3,])      % Eleccion de los ejes

Para apreciar mejor el desplazamiento de los puntos de la placa, representaremos el sólido antes y después del desplazamiento de sus puntos:

figure(5)
xd=xx+ux;                          % Posiciones de los puntos de la malla despues del desplazamiento horizontal
yd=yy+uy;                          % Posiciones de los puntos de la malla despues del desplazamiento verticales
subplot(1,2,1),mesh(xx,yy,0*xx)
axis([-4,4,-4,4])
view(2)
subplot(1,2,2), mesh(xd,yd,0*xd)
axis([-4,4,-4,4])
view(2)

4.2 Desplazamientos provocados en la placa

Una vez interpretado el campo vectorial, podemos obtener la posición de los puntos de la placa al actuar sobre ella; siendo la posición final de los puntos tras el desplazamiento, la suma de la posición inicial de los puntos del mallado más el desplazamiento (componentes de \(\vec u\))

h=0.1;                         % Paso de muestreo
rho=1:h:2;                     % Definición del intervalo de rho [1,2]
theta=0:h:2*pi+h;              % Definición del intervalo de theta [0,2*pi]
[rr,tt]=meshgrid(rho,theta);   % Definición de la placa en coordenadas polares

subplot(1,3,1)                 % Vamos a poner varios gráficos en una sola imagen
xx=rr.*cos(tt);                % Parametrización de las variables xx y yy en función de rr y tt
yy=rr.*sin(tt);
mesh(xx,yy,0*xx)               % Dibujamos la malla
title('Sólido inicial')
axis([-3,3,-3,3])              % Seleccionamos los ejes
view(2)                        % Vemos la imagen desde arriba

subplot(1,3,2)                 % Comenzamos el segundo gráfico en la misma imagen
ux=(((1-rr).^2).*(-yy));       % Componente x del campo vectorial u
uy=(((1-rr).^2).*(xx));        % Componente y del campo vectorial u
dx=xx+ux;                      % Suma del anillo inicial más la deformación producida por el campo en la dirección x
dy=yy+uy;                      % Suma del anillo inicial más la deformación producida por el campo en la dirección y
mesh(dx,dy,0*dx)               % Dibujamos la malla del sólido final
title('Sólido con el desplazamiento')
axis([-3,3,-3,3])              % Seleccionamos los ejes
view(2)                        % Vemos la imagen desde arriba

subplot(1,3,3)                 % Comenzamos el tercer gráfico en la misma imagen
mesh(dx,dy,0*xx)               % Dibujamos la malla del sólido final
axis([-3,3,-3,3])              % Seleccionamos los ejes
view(2)                        % Vemos la imagen desde arriba
hold on                        % Dibujamos dos gráficas en el mismo gráfico
quiver(xx,yy,ux,uy)            % Dibujamos el campo vectorial
title('Sólido con el desplazamiento y el campo del desplazamiento')
axis([-3,3,-3,3])              % Seleccionamos los ejes
view(2)                        % Vemos la imagen desde arriba

Al comparar ambos estados, el inicial y el final tras los desplazamientos, observamos que las partículas efectivamente se mueven en la dirección de los vectores del campo \(\vec u\) como habíamos supuesto.

Comparativa del anillo antes y después de haberse desplazado por la acción del campo vectorial u

4.3 Estudio de la divergencia [math]\nabla \cdot \vec u[/math]

Para analizar el cambio de volumen local debido a los desplazamientos anteriores, empleamos el operador diferencial de la divergencia:

La divergencia sirve para conocer el aumento de volumen de un cuerpo por efecto del desplazamiento. En nuestro caso, hemos calculado: [math]\bigtriangledown{}\cdot{}\overrightarrow{u} =\frac{1}{\rho}[(\frac{\partial}{\partial\theta})(\rho\cdot{}(1-\rho)^2)] =0 [/math]

Luego la placa no aumenta su volumen por efecto del desplazamiento que sufre. Aunque, como se ve en otros apartados, aumenta realmente su superficie, esto se debe a que es una placa, luego no aumenta su volumen.

Veamos como queda la divergencia en la placa a través del siguiente código:

h=0.1;
rho=1:h:2;                      % mallado del intervalo [1,2] con tamaño h= 0.1
theta=0:h:2*pi+h;               % mallado del intervalo [0,2*pi] con tamaño h=0.1
[rr,tt]=meshgrid(rho,theta);    % matrices de coordenadas u y v
xx=rr.*cos(tt);                 % parametrizacion
yy=rr.*sin(tt);

%Gráfica en 2D
surf(xx,yy,0*xx)

4.4 Estudio del rotacional [math]|\nabla \times \vec u|[/math]

Para analizar la tendencia del campo vectorial a inducir rotación en un punto, empleamos el operador diferencial del rotacional:

                       [math]\nabla \times \vec u=\frac{1}{√g}\begin{bmatrix}  \vec g_{\rho} & \vec g_{\theta} & \vec g_{z} \\ \frac{∂}{∂ \rho} & \frac{∂}{∂ \theta} & \frac{∂}{∂ z}  \\  u_{\rho} & u_{\theta} & u_{z} \end{bmatrix} [/math]

Al calcularlo, observamos que los puntos con mayor rotacional son aquellos en los que la malla ha sufrido el desplazamiento.

h=0.1;                         % Paso de muestreo
rho=1:h:2;                     % Definición del intervalo de rho [1,2]
theta=0:h:2*pi+h;              % Definición del intervalo de theta [0,2*pi]
[rr,tt]=meshgrid(rho,theta);   % Definición de la placa en coordenadas polares
xx=rr.*cos(tt);                % Parametrización de las variables xx y yy en función de rr y tt
yy=rr.*sin(tt);
zz=abs((4.*rr.^2)-(6.*rr)+2); 
 subplot(1,2,1);
 surf(xx,yy,zz)
 view(2)
 axis equal
 title('rotacional 2D')
 subplot(1,2,2);
 surf(xx,yy,zz)
 title('rotacional 3D')

Gráficamente comprobamos que tienen más rotacional los puntos exteriores, es decir, tienen más tendencia a rotar los puntos que han sufrido más desplazmiento.

Rotacional del campo vectorial en 2D Y 3D.

5 Estudio de las tensiones debidas a la perturbación

Si consideramos que estamos en un medio elástico lineal, isótropo y homogéneo, podemos definir el tensor de tensiones \(\sigma_{ij}\) a partir de los desplazamientos mediante de la expresión:

                          [math]\sigma_{ij}=\lambda \nabla \cdot \vec u \delta_{ij}+ 2\mu \epsilon_{ij}[/math]

donde [math]\epsilon_{ij}=\frac{(\nabla \vec u + \nabla \vec u^t)}{2}[/math], conocido tensor de deformaciones, es la parte simétrica del tensor gradiente de [math]\vec u[/math] y donde \(\lambda\) y \(\mu\) son conocidos como los coeficientes de Lamé que dependen de las propiedades elásticas de cada material. Tomando \(\lambda=\mu=1\), podemos proceder a calcular el tensor de tensiones

El gradiente se obtiene teniendo en cuenta que: [math]\nabla\vec u=u^i_{.j}\vec g_i \otimes \vec g^j[/math] Donde: [math]u^i_{j}= \vec g^i \nabla\vec u \vec g_j=\frac{\partial u^i}{\partial x^j}+\Gamma^i_{kj}u^k[/math]

En el caso particular de nuestro campo tenemos que [math]u^{\theta}=(1-\rho)^2, u^{\rho}=u^z=0[/math]. Obteniendo los símbolos de Christoffel para un sistema de coordenas cilíndricas resulta

[math]u^{\rho}_{\rho}=0 \qquad u^{\rho}_{\theta}=\Gamma^{\rho}_{\theta\theta}u^{\theta}=-\rho(1-\rho)^2 \qquad u^{\rho}_{z}=0 [/math]

[math]u^{\theta}_{\rho}=\frac{1}{\rho}(1-\rho)^2 \qquad u^{\theta}_{\theta}=0\qquad u^{\theta}_z=0[/math]

[math]u^z_{\rho}=u^z_{\theta}=u^z_{z}=0[/math]

Sabemos que la divergencia [math]\nabla \cdot \vec u =0[/math], por lo que podemos definir la matriz del tensor de tensiones en coordenadas 1-Contra 1-Cova:

                                       [math]\sigma^i_{.j}=\left(\begin{array}{ccc}0 & -\rho^3+2\rho^2+\frac{1}{\rho}-2&0\\-\rho^3+2\rho^2+\frac{1}{\rho}-2&0&0\\0&0&0\end{array}\right)[/math]

Para obtener el tensor en coordenadas 2-Covas multiplicamos su matriz por la izquierda por la matriz de Gram:

                                                   [math]G\sigma^i_{.j}=\sigma_{ij}=\left(\begin{array}{ccc}0 & -\rho^3+2\rho^2+\frac{1}{\rho}-2&0\\-\rho^5+2\rho^4+\rho-2\rho^2&0&0\\0&0&0\end{array}\right)[/math]

5.1 Tensiones normales

Por un lado, las tensiones normales en la dirección de \(\vec g_{\rho}\) se obtienen operando [math]\vec g_{\rho} \cdot \sigma \cdot \vec g_\rho[/math]; obteniendo:

                      [math]\vec g_{\rho} \cdot \sigma \cdot \vec g_\rho =\left(\begin{array}{ccc}1&0&0\end{array}\right)\left(\begin{array}{ccc}0 & -\rho^3+2\rho^2+\frac{1}{\rho}-2&0\\-\rho^5+2\rho^4+\rho-2\rho^2&0&0\\0&0&0\end{array}\right)\left(\begin{array}{ccc}1\\0\\0\end{array}\right)=(0)\vec g_{\rho} + (0)\vec g_{\theta} + (0)\vec g_{z}=0[/math]

Por otro lado, las tensiones normales en la dirección de \(\vec g_\theta/\rho\) se obtienen operando [math]\vec g_\theta/\rho \cdot \sigma \cdot \vec g_\theta/\rho[/math]; obteniendo:

                            [math]\vec g_\theta/\rho \cdot \sigma \cdot \vec g_\theta/\rho =\left(\begin{array}{ccc}0&\frac{1}{\rho}&0\end{array}\right)\left(\begin{array}{ccc}0 & -\rho^3+2\rho^2+\frac{1}{\rho}-2&0\\-\rho^5+2\rho^4+\rho-2\rho^2&0&0\\0&0&0\end{array}\right)\left(\begin{array}{ccc}0\\\frac{1}{\rho}\\0\end{array}\right)=(0)\vec g_{\rho} + (0)\vec g_{\theta} + (0)\vec g_{z}=0 [/math]

Puesto que las tensiones normales en ambos casos son nulas, a la hora de representarlas no obtendríamos ninguna gráfica significativa. Únicamente tendríamos la placa plana sobre la que se esta estudiando nuestro campo, sin ningún tipo de tensión aplicada. Seria la misma gráfica que la obtenida en la comparativa entre las tensiones tangenciales y el rotacional y la divergencia, que se muestra en el apartado 5.2 (para el calulo de la divergencia). Es por ello que no se incluye ningún programa de Matlab/Octave

5.2 Tensiones tangenciales

Las tensiones tangenciales al plano ortogonal respecto de [math]\vec g_{\rho}[/math] vienen dadas por la operación [math]|\sigma \cdot \vec g_\rho-(\vec g_\rho \cdot \sigma \cdot \vec g_\rho) \vec g_\rho|[/math]. Tal y como hemos calculado con anterioridad, el producto matricial [math]|-(\vec g_\rho \cdot \sigma \cdot \vec g_\rho)| [/math] es nulo, por lo que de la expresion inicial solo tenemos que realizar el producto matricial definido por:

                                 [math]|\sigma \cdot \vec g_\rho|=\left(\begin{array}{ccc}0 & -\rho^3+2\rho^2+\frac{1}{\rho}-2&0\\-\rho^5+2\rho^4+\rho-2\rho^2&0&0\\0&0&0\end{array}\right)\left(\begin{array}{ccc}1\\0\\0\end{array}\right)=|(1-\rho^5+2\rho^4+-2\rho+\rho)\vec g_{\rho}| [/math]

El programa de Matlab/Octave que nos permite dibujar las gráficas de las tensiones tangenciales respecto al plano ortogonal a Gp, para poder compararlas con las gráficas del rotacional y de la divergencia, con el objetivo de dar una interpretación física, es:

h=0.1;
rho=1:h:2;                      % mallado del intervalo [1,2] con tamaño h= 0.1
theta=0:h:2*pi+h;               % mallado del intervalo [0,2*pi] con tamaño h=0.1
[rr,tt]=meshgrid(rho,theta);    % matrices de coordenadas u y v
xx=rr.*cos(tt);                 % parametrizacion
yy=rr.*sin(tt);

%Graficas En 2D
figure(1)
subplot(2,3,1),
TensionTangencialGp=abs(-rr.^5+2.*rr.^4+rr-2.*rr.^2);    % Campo escalar de las tensiones tangenciales segun el vector g sub rho
surf(xx,yy,TensionTangencialGp)                          % mallado
axis([-3,3,-3,3])  
title('Tension Tangencial Gp')
view(2)   
subplot(2,3,2)
surf(xx,yy,0*xx)
title('Divergencia del Campo')
axis([-3,3,-3,3]) 
view(2)
subplot(2,3,3)
zz=abs((4.*rr.^2)-(6.*rr)+2); ;
surf(xx,yy,zz)
title('Rotacional del Campo')
view(2)

%Graficas en 3D

subplot(2,3,4)
surf(xx,yy,TensionTangencialGp) 
axis([-3,3,-3,3])  
subplot(2,3,5)
surf(xx,yy,0*xx) 
axis([-3,3,-3,3])  
subplot(2,3,6)
surf(xx,yy,zz) 
axis([-3,3,-3,3])

Las graficas obtenidas al ejecutar el código mostrado se exponen a continuacion:

Visualizacion en 2D y 3D de las tensiones tangenciales respecto al plano ortogonal a Gp, del rotacional y de la divergencia

A su vez, las tensiones tangenciales al plano ortogonal respecto de [math]\vec g_\theta/\rho[/math] se obtienen operando [math]|\sigma \cdot \vec g_\theta/\rho-(\vec g_\theta/\rho \cdot \sigma \cdot \vec g_\theta/\rho) \vec g_\theta/\rho|[/math]; Partiendo de la base de que [math]|-(\vec g_\theta/\rho \cdot \sigma \cdot \vec g_\theta/\rho) |[/math] ya hemos demostrado que tiene un valor nulo, solo tenemos que calcular el valor de [math]|\sigma \cdot \vec g_\theta/\rho|[/math]. Para ello operamos de la siguiente manera:

                                 [math]|\sigma \cdot \vec g_\theta/\rho| =\left(\begin{array}{ccc}0 & -\rho^3+2\rho^2+\frac{1}{\rho}-2&0\\-\rho^5+2\rho^4+\rho-2\rho^2&0&0\\0&0&0\end{array}\right)\left(\begin{array}{ccc}0\\\frac{1}{\rho}\\0\end{array}\right)=|(-\rho^2+2\rho+\frac{1}{\rho^2}-\frac{1}{\rho})\vec g_{\theta}|[/math]

El programa de Matlab/Octave que nos permite dibujar las gráficas de las tensiones tangenciales respecto al plano ortogonal a G0, para poder compararlas con las graficas del rotacional y de la divergencia, con el objetivo de dar una interpretación física, es:

h=0.1;
rho=1:h:2;                   % mallado del intervalo [1,2] con tamaño h= 0.1
theta=0:h:2*pi+h;            % mallado del intervalo [0,2*pi] con tamaño h=0.1
[rr,tt]=meshgrid(rho,theta); % matrices de coordenadas u y v
xx=rr.*cos(tt);              % parametrización
yy=rr.*sin(tt);

%Graficas en 2D
figure (2)
subplot(2,3,1)
TensionTangencialG0=abs(-rr.^2+2.*rr+1./(rr.^2)-2./(rr.^-1));    % Campo escalar de las tensiones tangenciales segun el vector g0/rho
surf(xx,yy,TensionTangencialG0)                                  % mallado
title('Tension Tangencial  G0')
axis([-3,3,-3,3]) 
view(2)       
subplot(2,3,2)
surf(xx,yy,0*xx)
title('Divergencia del Campo')
axis([-3,3,-3,3]) 
view(2)
subplot(2,3,3)
zz=abs((4.*rr.^2)-(6.*rr)+2); 
surf(xx,yy,zz)
title('Rotacional del Campo')
view(2)

%Graficas en 3D
subplot(2,3,4)
surf(xx,yy,TensionTangencialG0) 
axis([-3,3,-3,3])  
subplot(2,3,5)
surf(xx,yy,0*xx) 
axis([-3,3,-3,3])  
subplot(2,3,6)
surf(xx,yy,zz) 
axis([-3,3,-3,3])

Las graficas obtenidas al ejecutar el codigo mostrado son las siguientes:

Visulaizacion en 2D y 3D de las tensiones tangenciales respecto al plano ortogonal a G0, asi como de la divergenci y rotacionl

5.3 Tensión de Von Mises

La tensión de Von Mises se define mediante la expresión:

                              [math]\sigma_{VM} = \sqrt{\frac{(\sigma_1-\sigma_2)^2 + (\sigma_2-\sigma_3)^2 + (\sigma_3-\sigma_1)^2 }{2}}[/math]

donde \(\sigma_1\) , \(\sigma_2\) y \(\sigma_3\) son los autovalores de \(\sigma\), también conocido como tensor de tensiones. La tensión de Von Mises es una magnitud escalar que se suele usar como indicador para saber cuándo un material inicia un comportamiento plástico (y no elástico puro). Para poder representar dicha tension, ejecutaremos el programa de Matlab/Octave que a continuacion se muestra.

h=0.1;
rho=1:h:2;               % mallado del intervalo [1,2] con tamaño h= 0.1
theta=0:h:2*pi+h;            % mallado del intervalo [0,2*pi] con tamaño h=0.1
[rr,tt]=meshgrid(rho,theta);

%tension de von mises
m=length(rho);
s=zeros(3,i);
for i=1:m
  ro=rho(1,i);
  r=[0,-ro^3+2*ro^2+1/ro-2,0;-ro^3+2*ro^2+1/ro-2,0,0;0,0,0];
  lambda=eig(r);
  s(1,i)=lambda(1,1);
  s(2,i)=lambda(2,1);
  s(3,i)=lambda(3,1);
   end
s;
k=zeros(1,m);
for j=1:m
 a=s(1,j);
 b=s(2,j);
 c=s(3,j);
 tvm=sqrt(0.5*((a-b)^2+(b-c)^2+(c-a)^2));
 k(1,j)=tvm;
 end
k;

%dibujo de la tension

plot(rho,k)
title ('Representacion de la Tension Von Mises respecto a rho')

Resulta importante destacar que la tension obtenida varia en funcion de los valores que toma rho, puesto que el tensor de tensiones (calculado al inicio del punto 5) empleado para el calculo de la tension de Von Mises solo depende de dicha variable. Es por este motivo por el cual es suficiente representar los distintos valores que adopta la tension en funcion de la variacion de rho, es decir, pintar una grafica en 2D que nos de la curva que relaciona a ambos valores. Dicha representacion queda aqui mostrada:

Representacion de la Tension de Von Mises respecto a los diferentes valores que adopta la variable rho

6 Masa del cuerpo

La masa de la placa la calculamos con la integral doble (debido a que es una superficie) de la función densidad. Para calcularlo hay que hacer un cambio de variable a polares y su consiguiente jacobiano.

Para llegar a este resultado se usa el método del trapecio, para este caso se ha utilizado el siguiente código:

N1=200; N2=100;
%Número de puntos
r1=1; r2=2; t1=0; t2=pi;
%Extremos
h1=(r2-r1)/N1; h2=(t2-t1)/N2;
r=r1:h1:r2; t=t1:h2:t2;
%Coordenadas de la partición
[rr,tt]=meshgrid(r,t);
%Coordenadas del rectángulo
f=(rr.^3).*cos(tt).*sin(tt).*log(rr.*cos(tt)+2);
%Función
w1=ones(N1+1,1);
w(1)=1/2; w(N1+1)=1/2;
%Vectores de peso
w2=ones(N2+1,1);                 
w(1)=1/2; w(N2+1)=1/2;
%Vectores de peso
masa=h1*h2*w2'*f*w1          
%Resultado de la integral