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

De MateWiki
Saltar a: navegación, buscar
Trabajo realizado por estudiantes
Título Visualización de campos escalares y vectoriales en elasticidad. Grupo 1-C
Asignatura Teoría de Campos
Curso 2018-19
Autores José Luis Abia Pascual, Jorge Ismael Sanchez, Alberto Garces Rodriguez, Ricardo Perez
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

Consideramos una placa rectangular plana (en dimensión 2) que ocupa la región \((x,y) ∈ [0,2]×[-1,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)\)
  • Los desplazamientos [math]\vec u(x,y) [/math] producidos por la acción de una fuerza determinada.

De esta forma, si definimos [math]\vec r_{0}(x,y) [/math] el vector de posición de los puntos de la placa antes de la deformacón, la posición de cada punto \((x,y)\) de la placa después de la deformación viene dada por:

[math]\vec r(x,y) = \vec r_{0}(x,y) + \vec u(x,y) [/math]

Vamos a suponer que la fuerza aplicada sobre la placa ha provocado un desplazamiento de los puntos de la misma dado por el vector de desplazamientos:

[math]\vec u(\rho,\theta) = - \frac{\rho * \cos (\theta)}{10} \vec g_\theta [/math]

Dicho desplazamiento en coordenadas cartesianas seria:

[math]\vec u(x, y) = - \frac{\rho * \cos (\theta)}{10} * (- \rho * \sin (\theta) \vec i + \rho * \cos (\theta) \vec j) = - \frac{x}{10} * (-y \vec i + x \vec j) [/math]

Contenido

1 Dibujo del sólido

Dibujar un mallado que represente los puntos interiores del sólido. Tomar los ejes en el rectángulo \((x,y) ∈ [0,2]×[-1,1]\) y como paso de muestreo [math] h = \frac{1}{10} [/math] para las variables \((x,y)\).

Gráfica de la placa sin deformar
%Defimos los intervalos
%eje x
a=0:1/10:2;
%eje y
b=-1:1/10:1;
%Definimos la placa
[x,y]=meshgrid(a,b);
z=ones(size(x));
%dibujamos rectangulo
hold on
mesh(x,y,0*z)
title('Placa sin deformar');
axis([-1,3,-2,2])
view(0,90)
grid on
hold off


2 Temperatura del sólido

La temperatura del sólido viene dada por la función [math] T(x,y) = |\vec u(x, y)|^2 [/math] en grados centígrados. Dibujar las curvas de nivel y decidir en qué punto o puntos la temperatura es máxima y mínima.

2.1 Cálculo de la temperatura en todo el sólido

Cálculamos la temperatura en todo el sólido en coordenadas cartesianas:

[math] T(x,y) = |\vec u(x, y)|^2 = |- \frac{x}{10} * (-y \vec i + x \vec j)|^2 = \frac{x}{10} * \sqrt{x^2+y^2} [/math]

2.2 Visualización del Campo

Adjuntamos el código para visualizar las distintas temperaturas que hay en la placa. 

%definimos la ecuacion T(x,y)
T=x.*sqrt(x.^2+y.^2)./10;
%definimos curvas de nivel
figure 
hold on
mesh(x,y,0*z)
pcolor(x,y,T)
contour(x,y,T,'k')
shading flat
colorbar
hold off


Gráfica de la variación de la temperatura en toda la placa

2.3 Temperaturas máxima y mínima

De la matriz T creada en el apartado anterior sacariamos los valores maximo y mínimo de la temperatura de la placa. 

%definimos las temperaturas maximas y minimas
ma=max(max(T));
mi=min(min(T));
fprintf('La temperatura máxima es %.2f ºC y la mínima es %.2f ºC\n',ma,mi)


El programa nos devuelve 2 valores, el de máxima temperatura 45ºC que se encontraría en el borde inferior y superior derecho de la placa (La zona más calida del dibujo anterior) y el de mínima temperatura 0ºC que corresponderia al punto (0,0) (ya que es el que tiene el color más frío).

3 Gradiente de la temperatura

Calcular [math]\nabla T [/math] y dibujarlo como campo vectorial. Observar gráficamente que [math]\nabla T [/math] es ortogonal a dichas curvas.

3.1 Cálculo del gradiente de la temperatura

Hacemos el cálculo del gradiente de la temperatura en coordenadas cartesianas:

[math]\nabla T =\frac{\partial T}{\partial x_i} \vec e_i = \frac{\partial (\frac{x}{10} * \sqrt{x^2+y^2})}{\partial x} \vec i + \frac{\partial (\frac{x}{10} * \sqrt{x^2+y^2})}{\partial y} \vec j + \frac{\partial (\frac{x}{10} * \sqrt{x^2+y^2})}{\partial z} \vec k = \frac{2x^2+y^2}{10\sqrt{x^2+y^2}} \vec i + \frac{xy}{10\sqrt{x^2+y^2}} \vec j + 0 \vec k [/math]

3.2 Visualización del gradiente de la temperatura

El cálculo hecho anteriormente para calcular el gradiente lo pasamos a matlab y generamos una representación de dicho campo. 

%Definimos el gradiente
Tx=@(x,y)(2*x.^2+y.^2)./sqrt(x.^2+y.^2)/10;
Ty=@(x,y)x.*y./sqrt(x.^2+y.^2)/10;
%Tk=0*z;
Ti=Tx(x,y);
Tj=Ty(x,y);
%dibujamos el gradiente
figure
hold on
contour(x,y,T)
quiver(x,y,Ti,Tj)
hold off

Podemos observar que los vectores son perpendenculares a las curvas de nivel de la temperatura.

Gráfica del campo gradiente de la placa

4 Variación de la Temperatura

Calcular la variación de temperatura por segundo si, partiendo del punto de coordenadas \((x,y)=(0,1)\) me muevo en la dirección \((1,-1)\) con una velocidad de 2 metros por segundo.

  • Planteamiento:
  1. Partiendo del origen \(P=(0,1)\) nos movemos en la dirección \((1,-1)\), que se representa con el vector unitario [math] \vec e = \frac{\vec i - \vec j}{|\vec i - \vec j|} = \frac{\vec i - \vec j}{\sqrt{2}}[/math] con velocidad \(v=2\)m/s.
  2. De las propiedades del gradiente de un campo escalar:
[math] \frac{\partial T}{\partial t} = \frac{T(\vec r (t))}{\partial t} = \nabla{T(\vec r (t))} · \vec r ' = \nabla{T} · \vec v = \nabla{T} · (v \vec e) = v(\nabla{T} · \vec e)[/math]
y puesto que ya conocemos el \(\nabla{T}\):
[math] \nabla{T}(x,y) = \frac{2x^2+y^2}{10\sqrt{x^2+y^2}} \vec i + \frac{xy}{10\sqrt{x^2+y^2}} \vec j [/math]

Entonces operando se obtiene:

[math] \frac{\partial T}{\partial t} = v·\nabla{T}(x,y)·\vec e = 2*(\frac{2x^2+y^2}{10\sqrt{x^2+y^2}} \vec i + \frac{xy}{10\sqrt{x^2+y^2}} \vec j)* (\frac{\vec i - \vec j}{\sqrt{2}}) = \sqrt{2}*(\frac{2x^2+y^2}{10\sqrt{x^2+y^2}} - \frac{xy}{10\sqrt{x^2+y^2}}) = \frac{\sqrt{2}}{10}(\frac{2x^2+y^2-xy}{\sqrt{x^2+y^2}}) [/math] ºC/s


que es la variación de temperatura por segundo en la dirección indicada.

5 Campo de desplazamientos

Dibujar el campo de desplazamientos en los puntos del mallado del sólido.

5.1 Visualización del Campo de desplazamientos

Estudiarémos como el vector de desplazamientos [math] \vec u(x,y) [/math] deformaría la placa, esto lo interpretaremos con la dirección que cogen los vectores en la placa sin deformar. 

%Discretización de los parámetros de la superficie:
a=0:1/10:2;
b=-1:1/10:1;
%Creación del mallado
[x,y]=meshgrid(a,b);
%Componentes en la dirección de i y de j del campo de desplazamientos 
ux=+x/10.*(y);
uy=-x/10.*(x);
figure
%dibujo del mallado 
mesh(x,y,0*x)
hold on 
%campo de desplazamientos
quiver(x,y,ux,uy)
axis([-1,3,-2,2])
view(0,90)
hold off


Gráfica de los desplazamientos de la placa

6 Estudio del comportamiento del sólido durante el desplazamiento

Dibujar el sólido antes y después del desplazamiento dado por el campo de vectores [math] \vec u(x,y) [/math]. Dibujar ambos en la misma fígura usando el comando subplot.

6.1 Visualización del sólido antes del desplazamiento

Dibujamos la placa sin deformar [math]\vec r(x,y) = \vec r_{0}(x,y) [/math].

Gráfica de antes de la deformación de la placa
%Discretización de los parámetros de la superficie:
a=0:1/10:2;
b=-1:1/10:1;
%Creación del mallado
[x,y]=meshgrid(a,b);
%posicion final 
rx=x+x.*y/10;
ry=y-x.*x/10; 
%representacion de la superficie antes del desplazamiento  
figure
%subplot(3,1,1)
surf(x,y,0*x)
title('antes del desplazamiento')
axis([-1,3,-2,2])
view(2) 
xlabel('Eje X')
ylabel('Eje Y') 
zlabel('Eje Z')


6.2 Visualización del sólido después del desplazamiento

Gráfica de después de la deformación de la placa

Dibujamos la placa después de la deformación [math]\vec r(x,y) = \vec r_{0}(x,y) + \vec u(x,y) [/math]. 

%representacion de la superficie después del desplazamiento
% subplot(3,1,2) 
figure
surf(rx,ry,0*rx)
%axis([-1,3,-2,2])
title('después del desplazamiento')
axis([-1,3,-2,2])
xlabel('Eje X')
ylabel('Eje Y') 
zlabel('Eje Z') 
view(0,90)


6.3 Visualización del sólido antes y después del desplazamiento

Comparamos el antes [math]\vec r(x,y) = \vec r_{0}(x,y) [/math] y el después [math]\vec r(x,y) = \vec r_{0}(x,y) + \vec u(x,y) [/math] de la deformación. 

%comparacion del antes y después
%subplot(3,1,3)
figure
surf(x,y,0*x)
hold on 
surf(rx,ry,0*rx)
axis([-1,3,-2,2])
xlabel('Eje X')
ylabel('Eje Y') 
zlabel('Eje Z') 
title('Comparación antes y después')
view(0,90)
hold off


En esta imagen podemos apreciar como la placa se deforma según el campo de desplazamientos visto en el apartado 5.

Gráfica de comparando el antes y después de la deformación de la placa

7 Divergencia del campo de desplazamientos

Dibujar [math] \nabla \cdot \vec u [/math]. Determinar gráficamente los puntos en los que la divergencia de [math] \vec u(x,y)= \frac{yx}{10} \vec i + \frac{-x^2}{10} \vec j [/math] es máxima, mínima y nula. La divergencia es una medida del cambio de volumen local debido al desplazamiento.

7.1 Cálculo de la Divergencia

Cálculamos la divergencia en coodenadas cartesianas:

[math] \nabla \cdot \vec u =\frac{\partial \vec u_1}{\partial x} + \frac{\partial \vec u_2}{\partial y} + \frac{\partial \vec u_3}{\partial z} =\frac{\partial \frac{yx}{10}}{\partial x} + \frac{\partial \frac{-x^2}{10}}{\partial y} + \frac{\partial 0}{\partial z} =\frac{y}{10} + 0 + 0 =\frac{y}{10} [/math]

7.2 Visualización de la Divergencia

Mediante el siguiente código dibujamos la divergencia y donde su variación es máxima, mínima y nula: 

%Definimos la placa
x=0:.1:2;
y=-1:.1:1;
%Mallamos la placa
[X,Y]=meshgrid(x,y);
%Damos valores a la divergencia
Z=(1/10)*Y;
z=y/10;
%Dibujamos la divergencia
figure
pcolor(X,Y,Z)
shading interp
colorbar
xlabel('Eje X')
ylabel('Eje Y')
title('Variación de la Divergencia a lo largo de la placa')
%Vemos los valores máximos, mínimos y nulos
figure
plot(y,z)
hold on
xlabel('Eje Y')
ylabel('Eje Z')
plot(1,0.1,'xr','markersize',10)
plot(0,0,'xg','markersize',10)
plot(-1,-0.1,'xb','markersize',10)
legend('gráfica de la diverencia','punto máximo','punto nulo','punto mínimo');
hold off


Gráficas de la variación de la divergencia

8 Rotacional del campo de desplazamientos

Calcular [math] \mid \nabla \times \vec u \mid [/math] en todos los puntos del sólido y dibujarlo. ¿Qué puntos sufren un mayor rotacional?

8.1 Cálculo del Rotacional

Cálculamos el rotacional

[math] \nabla \times \vec u=\left| \begin{matrix} \vec i & \vec j & \vec k \\ & & \\ \frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \\ & & \\ u_x & u_y & u_z \end{matrix}\right| =\left| \begin{matrix} \vec i & \vec j & \vec k \\ & & \\ \frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \\ & & \\ \frac{xy}{10} & \frac{-x^2}{10} & 0 \end{matrix}\right| =\left( \frac{\partial 0}{\partial y} - \frac{ \partial \frac{-x^2}{10}}{\partial z}\right)\vec i + \left(\frac{\partial \frac{xy}{10}}{\partial z} - \frac{\partial 0}{\partial x}\right)\vec j + \left(\frac{ \partial \frac{-x^2}{10}}{\partial x} - \frac{ \partial \frac{xy}{10}}{\partial y}\right)\vec k =\left(0 - 0\right)\vec i + \left(0 - 0\right)\vec j + \left(\frac{-2x}{10} - \frac{x}{10}\right)\vec k =-\frac{3x}{10}\vec k [/math]

Cálculamos el módulo del rotacional

[math] \mid \nabla \times \vec u \mid = \mid -\frac{3x}{10}\vec k \mid = \frac{3x}{10} [/math]

8.2 Visualización del Rotacional

Mediante el siguiente código dibujamos el módulo del Rotacional: 

%Definimos la placa
x=0:.1:2;
y=-1:.1:1;
%Mallamos la placa
[X,Y]=meshgrid(x,y);
%Damos valores al módulo del rotacional
Z=3*X/10;
%Dibujamos el rotacional
figure
pcolor(X,Y,Z)
shading interp
colorbar
xlabel('Eje X')
ylabel('Eje Y')
title('Módulo del rotacional')


Gráfica del módulo del rotacional

Como podemos observar en la imagen anterior los puntos que sufren mayor rotacional serian todos los puntos del extremo derecho de la placa

9 Tensiones Normales

9.1 Cálculo del tensor de deformaciones

El tensor de deformaciones [math]\epsilon (\vec u) [/math] está definido como la parte simétrica del tensor gradiente de [math]\vec u [/math], [math]\nabla \vec u [/math] . Por tanto, para calcularlo debemos calcular el gradiente de [math]\vec u [/math] y su traspuesto, [math]\nabla \vec u ^ t [/math].

[math]\vec u(x, y) = \frac{xy}{10} \vec i + \frac{-x^2}{10} \vec j [/math]


[math]\nabla \vec u = \frac{\partial u}{\partial x} \vec i + \frac{\partial u}{\partial y} \vec j = ( \frac{y}{10}\vec i - \frac{2x}{10}\vec j)\cdot\vec i + (\frac{x}{10}\vec i)\cdot\vec j = \frac{y}{10}\vec i \otimes \vec i - \frac{2x}{10}\vec i \otimes \vec j + \frac{x}{10}\vec j \otimes \vec i[/math]


[math]\nabla \vec u ^ t = \frac{y}{10}\vec i \otimes \vec i + \frac{x}{10}\vec i \otimes \vec j - \frac{2x}{10}\vec j \otimes \vec i [/math]


Así, el tensor de deformaciones quedará:


[math]\epsilon (\vec u) = \frac{\nabla \vec u + \nabla \vec u ^ t}{2} = \frac{y}{10}\vec i \otimes \vec i - \frac{x}{20}\vec i \otimes \vec j - \frac{x}{20}\vec j \otimes \vec i [/math]

9.2 Cálculo del tensor de tensiones [math]\sigma_{ij}[/math]

Cuando nos encontramos ante un medio lineal, isótropo y homogéneo, como nuestra placa; podemos definir el tensor de tensiones [math]\sigma_{ij}[/math] mediante la siguiente expresión:


[math]\sigma_{ij} = \lambda \nabla \cdot \vec u \boldsymbol {1} + 2 \mu \epsilon [/math]


En esta expresión aparecen los coeficientes de Lamé, [math]\lambda [/math] y [math] \mu [/math], que determinan el comportamiento elástico del sólido en pequeñas deformaciones, y que dependen del módulo de Young ([math] E [/math]) y del parámetro de Poisson ([math] \upsilon [/math]).


[math] \lambda = \frac {E\upsilon}{(1+\upsilon)(1-2\upsilon)} [/math]

[math] \mu = \frac {E}{(1+\upsilon)} [/math]


Vamos a tomar como valores de estos coeficientes [math]\lambda = \mu = 1 [/math] para el cálculo del tensor de tensiones. Necesitamos la divergencia del campo de desplazamientos, [math]\nabla \cdot \vec u[/math]:


[math] \nabla \cdot \vec u = \frac{\partial u_1}{\partial x} \vec i + \frac{\partial u_2}{\partial y} = \frac {y}{10} [/math]


Y así, nos queda el siguiente campo de tensiones:


[math] \sigma_{ij} = \lambda \nabla \cdot \vec u \boldsymbol {1} + 2 \mu \epsilon = \frac {y}{10} ( \vec i \otimes \vec i + \vec j \otimes \vec j + \vec k \otimes \vec k ) + 2 ( \frac {y}{10} \vec i \otimes \vec i - \frac {x}{20}\vec i \otimes \vec j - \frac {x}{20} \vec j \otimes \vec i ) = \frac {3y}{10} \vec i \otimes \vec i - \frac {x}{10} \vec i \otimes \vec j - \frac {x}{10} \vec j \otimes \vec i + \frac {y}{10} \vec j \otimes \vec j + \frac {y}{10} \vec k \otimes \vec k) [/math]


[math] \sigma_{ij} = \begin{pmatrix} \frac {3y}{10} & - \frac {x}{10} & 0\\ - \frac {x}{10} & \frac {y}{10} & 0\\ 0 & 0 & \frac {y}{10} \end{pmatrix} [/math]

9.3 Cálculo de las tensiones normales

Para calcular las tensiones normales, vamos a proyectar el tensor de tensiones [math]\sigma_{ij} [/math] sobre los ejes.

9.3.1 Tensiones normales en la dirección que marca el eje [math]\vec i [/math]

Las tensiones normales en la dirección que marca el eje [math]\vec i [/math] se obtienen a partir de la expresión:

[math] \vec i \cdot \sigma_{ij} \cdot \vec i = \frac {3y}{10} [/math]

9.3.2 Tensiones normales en la dirección que marca el eje [math]\vec j [/math]

Las tensiones normales en la dirección que marca el eje [math]\vec j [/math] se obtienen a partir de la expresión:

[math] \vec j \cdot \sigma_{ij} \cdot \vec j = \frac {y}{10} [/math]

9.3.3 Tensiones normales en la dirección que marca el eje [math]\vec k [/math]

Las tensiones normales en la dirección que marca el eje [math]\vec k [/math] se obtienen a partir de la expresión:

[math] \vec k \cdot \sigma_{ij} \cdot \vec k = \frac {y}{10} [/math]

9.4 Visualización de las tensiones normales

Mediante el siguiente código dibujamos las tensiones normales en los ejes [math] \vec i, \vec j , \vec k [/math]: 

%Discretización de los parámetros de la superficie:
%Siendo 0.1 el "paso de muestreo"
a=0:1/10:2;
b=-1:1/10:1;
%Creación del mallado
[x,y]=meshgrid(a,b);
%Tensiones normales en las direcciones de los ejes coordenados sobre el mallado 
TNi=y/10*3;
TNj=y/10;
TNk=y/10; 
%Representacion de las tensiones normales 
%Tensiones normales en la dirección del eje i
figure
subplot(1,3,1) 
pcolor(x,y,TNi)
colorbar
xlabel('Eje X')
ylabel('Eje Y') 
zlabel('Eje Z')
title('Tensiones normales en la dirección del eje i')  
%Tensiones normales en la dirección del eje j
subplot(1,3,2) 
pcolor(x,y,TNj)
colorbar
xlabel('Eje X')
ylabel('Eje Y') 
zlabel('Eje Z')
title('Tensiones normales en la dirección del eje j')
%Tensiones normales en la dirección del eje k
subplot(1,3,3)
pcolor(x,y,TNk)
colorbar
xlabel('Eje X')
ylabel('Eje Y') 
zlabel('Eje Z')
title('Tensiones normales en la dirección del eje k')


Representación de las tensiones normales en los ejes [math] \vec i, \vec j , \vec k [/math]

10 La tensión de Von Mises

10.1 Cálculo e Interpretación de la tensión de Von Mises

La tensión de Von Mises es una magnitud escalar que se suele usar como indicador para saber cuando un material inicia un comportamiento plástico. Está definida por la fórmula

[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 [/math], [math] \sigma_2 [/math] y [math] \sigma_3 [/math] son las tensiones principales, que se corresponden con los autovalores de [math] \sigma_{ij} [/math].

10.2 Visualización de la tensión de Von Mises

Donde [math]\sigma_{1}, \sigma_{2}[/math] y [math]\sigma_{3}[/math] son los autovalores (también conocidos como tensiones principales) de [math]\sigma[/math] cuya expresión es:


[math] \sigma_{ij} = \begin{pmatrix} \frac {3y}{10} & - \frac {x}{10} & 0\\ - \frac {x}{10} & \frac {y}{10} & 0\\ 0 & 0 & \frac {y}{10} \end{pmatrix} [/math]


La expresión de dichos autovalores, calculada analiticamente, sería:


[math]\sigma=\begin{pmatrix} \frac{2y + \sqrt{x^2+y^2} }{10} & 0 & 0 \\ 0 & \frac{2y - \sqrt{x^2+y^2} }{10} & 0 \\ 0 & 0 & \frac{y}{10} \end{pmatrix}[/math]


El máximo valor de la tensión de Von Mises se alcanza para [math]x=2[/math] y el valor de la tensión de Von Mises para dicho punto es de \(0.4\). 

%Discretización de los parámetros de la superficie:
% siendo 0.1 el "paso de muestreo"
a=0:1/10:2;
b=-1:1/10:1;
%Creación del mallado
[x,y]=meshgrid(a,b);
%calculo de los autovalores 
TX=@(x,y)[y/10*3 -x/10 0];
TY=@(x,y)[-x/10 y/10 0];
TZ=@(x,y)[0 0 y/10];
%Obtención de autovalores de la tensión de Von Mises
for i=1:length(a)
    for j=1:length(b)
    %Obtención de valores para cada componente del mallado para las funciones de las componentes de la función de tensión de Von Mises    
    U=TX(x(i,j),y(i,j));
    V=TY(x(i,j),y(i,j));
    W=TZ(x(i,j),y(i,j));
    %Creación de un vector cuyas componentes sean los vectores en los que se han almacenado los valores de las funciones estudiadas
    vec=[U; V; W]; 
    %Obtención de los autovalores de la tensión de Von Mises
    auto=eig(vec);
    %Cálculo de la tensión de Von Mises para cada componente del mallado 
    VM=sqrt(((auto(1)-auto(2))^2+(auto(2)-auto(3))^2+((auto(3)-auto(1))^2))*1/2);
    z(i,j)=VM;
    end
end
%dibujamos Tensiones de Von Mises
figure
subplot(1,2,1)
surf(x,y,z)
colorbar
xlabel('Eje X')
ylabel('Eje Y')
zlabel('Eje Z')
title('Tensión de Von Mises')
%Estudio del máximo valor de la tensión de Von Mises
subplot(1,2,2)
hold on
m=z(1,1:length(z));
%Dibujo de la proyección sobre el plano XOZ de la superficie que genera la tensión de Von Mises
plot(a,m,'b');
maxz=max(m);
%Bucle que dibuja los puntos donde la tensión de Von Mises es máxima sobre la proyección
for k=1:length(m)
   if m(k)==maxz
         plot(a(k),maxz,'xr','markersize',10)
     end
end
axis([-0.5 2.5 0 0.5]);
xlabel('Eje X')
ylabel('Eje Z')
legend('Proyección sobre el plano XOZ','Punto máximo')
hold off
fprintf('El máximo valor de la tensión de Von Mises es %f\n',maxz)


Representación de las tensiones de Von Mises

11 Campo de Fuerzas [math]\vec{F}[/math] que actua sobre el sólido

El campo de fuerzas [math]\vec{F}[/math] que actúa sobre la placa (y que son las causantes del desplazamiento observado) se aproxima usando la ecuación de la elasticidad lineal:

[math] \vec{F} = - \nabla \cdot \sigma = - \frac {\partial \sigma_{ji}}{\partial x_j} \cdot \vec e_i [/math]

Dibujar [math]\vec{F}[/math] e interpretar la grafíca.

11.1 Cálculo del Campo de Fuerzas [math]\vec{F}[/math]

El desplazamiento de la placa dado por el vector [math]\vec u(\rho,\theta) = - \frac{\rho * \cos (\theta)}{10} \vec g_\theta = - \frac{x}{10} * (-y \vec i + x \vec j)[/math] está causado por un campo de fuerzas [math]\vec{F}[/math] que actúa sobre el interior de la placa.

Para calcularlo, vamos a utilizar la ecuación de la elasticidad lineal [math]\vec{F} = - \nabla \cdot \sigma_{ij} = - \frac {\partial \sigma_{ji}}{\partial x_j} \cdot \vec e_i [/math]. Es decir, que para calcular las componentes del campo de fuerzas haremos la divergencia de los vectores fila que forman la matriz [math]\sigma_{ij}[/math].


Así, obtenemos las siguientes componentes:

[math] \vec{F_1} = - \frac {\partial \sigma_{j1}}{\partial x} \cdot \vec i = - ( \frac {\partial \frac{3y}{10} }{\partial x} + \frac {\partial \frac{-x}{10} }{\partial y} + \frac {\partial 0 }{\partial z} ) \cdot \vec i = 0 \vec i [/math]

[math] \vec{F_2} = - \frac {\partial \sigma_{j2}}{\partial y} \cdot \vec j = - ( \frac {\partial \frac{-x}{10} }{\partial x} + \frac {\partial \frac{y}{10} }{\partial y} + \frac {\partial 0 }{\partial z} ) \cdot \vec j = (-\frac {1}{10} + \frac {1}{10}) \vec j = 0 \vec j [/math]

[math] \vec{F_3} = - \frac {\partial \sigma_{j3}}{\partial z} \cdot \vec k = - ( \frac {\partial 0}{\partial x} + \frac {\partial 0}{\partial y} + \frac { \partial \frac {y}{10}}{\partial z} ) \cdot \vec k = 0 \vec k [/math]

11.2 Visualización del Campo de Fuerzas [math]\vec{F}[/math]

Mediante el siguiente código dibujamos el campo de fuerzas [math]\vec{F}[/math]: 

%Discretización de los parámetros de la superficie:
% siendo 0.1 el "paso de muestreo"
a=0:1/10:2;
b=-1:1/10:1;
%Creación del mallado
[x,y]=meshgrid(a,b);
%componentes del campo de fuerzas
Fx=zeros(size(x)); 
Fy=zeros(size(x));
%Dibujo del campo de fuerzas 
figure
mesh(x,y,0*x) 
hold on 
quiver(x,y,Fx,Fy,'k')
view(0,90)
xlabel('Eje X')
ylabel('Eje Y')
title('Campo de fuerzas')
hold off


Representación del campo de fuerzas [math]\vec{F}[/math]

11.3 Interpretación del Campo de Fuerzas [math]\vec{F}[/math]

El campo de fuerzas que causa los desplazamientos que estudiamos es [math]\vec{F} = \vec 0 [/math]. Esto, aparentemente, no tiene sentido, ya que los desplazamientos causados por [math]\vec u(x, y) = \frac{xy}{10} \vec i + \frac{-x^2}{10} \vec j [/math] son distintos de [math] \vec 0 [/math].

Sin embargo, es necesario tener en cuenta que además del campo de fuerzas que actúa en el interior de la plaza, también existen unas fuerzas en la frontera (en los bordes de la placa) que pueden causar desplazamientos en la misma, aunque el campo de fuerzas interiores sea [math]\vec{F} = \vec 0 [/math]. Estas fuerzas se pueden calcular multiplicando el tensor de tensiones [math] \sigma_{ij} [/math] por los vectores normales a las fronteras de la placa, es decir, [math]\vec i[/math] y [math]\vec j[/math].

12 Masa total de la placa

Supongamos que la densidad de la placa es [math] d(x,y)=1+e^{\frac{-|x|^2}{(y+1)^{4}}} [/math]. Calcular la masa total aproximando la integral correspondiente numéricamente.

12.1 Cálculo de la masa total de la placa

Calculamos la masa total de la placa:

[math]\int_{A} d(x,y) dA=\int_{0}^{2}\int_{-1}^{1} d(x,y) dydx=\int_{0}^{2}\int_{-1}^{1} 1+e^{\frac{-|x|}{(y+1)^{4}}} dydx=5.7973 [/math]

Como la primitiva es difícil de calcular, vamos a utilizar métodos numéricos, aproximando la integral mediante la regla del trapecio: 

%Extremos de integracion (debido a que el punto y=-1 es singular en la función de densidad, aproximamos a -0.999999999): 
a=0; b=2; c=-0.999999999; d=1; 
%Subintervalos
N1=200; N2=N1; 
h1=(b-a)/N1; h2=(d-c)/N2;
%Mallado de la superficie
x=a:h1:b; 
y=c:h2:d; 
[u,v]=meshgrid(x,y);
%Función de densidad 
d=1+exp(-abs(u).^2./((v+1).^4));
%formula del trapecio 
w1=ones(N1+1,1);
w1(1)=1/2; 
w1(N1+1)=1/2;
w2=ones(N2+1,1);
w2(1)=1/2;
w2(N2+1)=1/2; 
M=h1*h2*(w2')*d*w1;
fprintf('El valor de la masa es %.4f', M)
Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Visualización_de_campos_escalares_y_vectoriales_en_elasticidad_(Grupo_C1)&oldid=42271»