Grupo 1c: visualizacion de campos escalares y vectoriales en elasticidad

Trabajo realizado por estudiantes
Título	Trabajo campos grupo 1C: Visualización de campos escalares y vectoriales en elasticidad.
Asignatura	Teoría de Campos
Curso	2019-20
Autores	Marta Lozano Martinez, Candela Carmen Martin Blanco, Ilenia Maria Morales Perez, Daniel Segura Santos
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) ∈ [-2,2]×[0,4]\). 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(x,y) = -\frac{x^2}{20} \vec i + f(y) \vec j [/math].

Donde f(y) es cierta funcion

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) ∈ [-3,3]×[0,5]\) y como paso de muestreo [math] h = \frac{1}{10} [/math] para las variables \((x,y)\).

Mallado del sólido

%Discretización de los parámetros de la superficie:
% siendo 0.1 el "paso de muestreo"
 x=-2:0.1:2;  
 y=0:0.1:4; 
% creacion del mallado     
 [xx,yy]=meshgrid(x,y); 
% dibujo
 figure(1)
 mesh(xx,yy,0*xx) 
% ejes    
 axis([-3,3,0,5])      
view(2)

2 Visualización campo de la temperatura

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

Campo de la Temperatura

%Discretización de los parámetros de la superficie:
% siendo 0.1 el "paso de muestreo"
x=-2:0.1:2;  
y=0:0.1:4;   
% creación del mallado           
[xx,yy]=meshgrid(x,y);
figure(1)
% campo escalar
T=(xx+2).^2+(yy+2).^2; 
% dibujo
contour(xx,yy,T)       
hold on               
plot(x,x-x,'k','linewidth',1);
plot(x,4+x-x,'k','linewidth',1);
plot(-2+y-y,y,'k','linewidth',1);
plot(2+y-y,y,'k','linewidth',1);
axis([-3,3,0,5])      
view(2)               
colorbar

Se observa que en el punto \((2,4)\) la temperatura es máxima.

3 Cálculo y visualización del gradiente de Temperatura

Calcular el gradiente la temperatura (campo escalar) y dibujarlo como el campo vectorial. Observar gráficamente que el gradiente de la temperatura es ortogonal a dichas curvas. utilizaremos la misma escala de los ejes que el apartado anterior par apreciar el ángulo que forman las curvas de nivel y el gradiente

[math]\nabla T =\frac{\partial T}{\partial x_i} \vec e_i = \frac{\partial T}{\partial x}\vec i + \frac{\partial T}{\partial y} \vec j = (2+2x)\vec i + (2+2y) \vec j [/math]

Como se puede observar en el dibujo el gradiente es ortogonal a las curvas de nivel del campo \(T\) y por tanto el [math]\nabla T[/math] es ortogonal a la superficie del campo \(T\) en todos los puntos de éste.

Campo vectorial: Gradiente de Temperatura

%Discretización de los parámetros de la superficie:
% siendo 0.1 el "paso de muestreo"
x=-2:0.1:2;  
y=0:0.1:4;
% mallado              
[xx,yy]=meshgrid(x,y); 
figure(1)
% temperatura 
T=(xx+2).^2+(yy+2).^2; 
%dibujo
contour(xx,yy,T)      
hold on                
plot(x,x-x,'k','linewidth',1);
plot(x,4+x-x,'k','linewidth',1);
plot(-2+y-y,y,'k','linewidth',1);
plot(2+y-y,y,'k','linewidth',1);
axis([-3,3,0,5])      
view(2)                
colorbar 
hold on
x=-2:0.1:2;     
y=0:0.1:4;           
[xx,yy]=meshgrid(x,y); 
figure(1)
fx=(2.*xx+2);                
fy=2.*yy+2;                 
quiver(xx,yy,fx,fy)     
axis([-3,3,0,5])     
view(2)

4 Campo de desplazamientos

Calcular el campo de desplazamientos [math]\vec u [/math], sabinedo:

• Los puntos situados en el eje y=0 no sufren desplazamientos en la dirección [math]\vec j [/math]
• [math]\nabla ·\vec u = -\frac{x}{10} + \frac{1}{10} [/math]

En primer lugar hayamos la divergencia de :[math]\vec u(x,y) = -\frac{x^2}{20} \vec i + f(y) \vec j [/math] y la igualamos con la condición que nos han dado:

[math]\nabla ·\vec u = \frac{\partial u}{\partial x} + \frac{\partial u}{\partial y} = -\frac{x}{10} + f'(y) [/math]

Igualando: [math] -\frac{x}{10} + f'(y) = -\frac{x}{10} + \frac{1}{10} [/math]
Si despejamos f'(y), nos queda que f'(y) = [math] \frac{1}{10} [/math] y por lo tanto f(y) seria la integral de f'(y) respecto y, lo que nos da que f(y) = [math] \frac{y}{10} [/math]
Es decir: [math]\vec u(x,y) = -\frac{x^2}{20} \vec i + \frac{y}{10} \vec j [/math]

5 Campo de vectores

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

Campo de vectores

%Discretización de los parámetros de la superficie:
% siendo 0.1 el "paso de muestreo"
 x=-2:0.1:2;
 y=0:0.1:4;
 %Creación del mallado
[xx,yy]=meshgrid(x,y);
 %Componentes en la dirección de i y de j del campo de desplazamientos 
 ux=(-xx.^2)/20;
 uy=yy/10;
 figure
 %dibujo del mallado 
mesh(xx,yy,0*xx)
 hold on 
 %campo de desplazamientos
 quiver(xx,yy,ux,uy,'k')
 axis([-3,3,0,5])
 view(0,90)
 hold off

6 Representación del sólido antes y después del desplazamiento

Dibujar el sólido antes y después del desplazamiento dado por el campo de vectores [math]\vec u [/math]. Dibujar ambos en la misma figura usando el comando subplot.

Representación del sólido antes y después del desplazamiento

%Discretización de los parámetros de la superficie:
  % siendo 0.1 el "paso de muestreo"
  x=-2:0.1:2;
  y=0:0.1:4;
  %Creación del mallado
  [xx,yy]=meshgrid(x,y);
  %posicion final 
  rx=((-xx.^2)/20)+xx;
  ry=(yy/10)+yy; 
 %representacion de la superficie antes del desplazamiento  
 figure
 subplot(1,3,1)
 surf(xx,yy,0*xx)
 title('antes del desplazamiento')
 axis([-2.5,2.5,0,5])
 view(2) 
 xlabel('x')
 ylabel('y') 
 zlabel('z')
 %representacion de la superficie después del desplazamiento
 subplot(1,3,2) 
 surf(rx,ry,0*rx)
 axis([-2.4,2.5,0,5])
 title('después del desplazamiento')
 axis([-3,3,0,5])
 xlabel('x')
 ylabel('y') 
 zlabel('z')
  
 %comparacion
 subplot(1,3,3)
 surf(xx,yy,0*xx)
 hold on 
 surf(rx,ry,0*rx)
 axis([-2,2,-1,3])
 xlabel('x')
 ylabel('y') 
 zlabel('z') 
 title('comparacion')
 view(0,90)
 hold off

Como podemos observar en las gráficas adjuntas, el sólido se encuentra posicionado en el eje x entre -2 y 2, mientras que en el eje y entre 0 y 4. Después de sufrir un desplazamiento, la placa se desplaza en sentido izquierdo del eje x debido al símbolo negativo presente en la componente x del desplazamiento [math]\vec u [/math] , mientras que en el eje y se estira un 10% de su longitud situándose la componente y en 4.4. Gracias al último gráfico podemos comparar este desplazamiento.

7 Divergencia máxima, mínima y nula

Dibujar ∇ · [math]\vec u [/math]. Determinar analíticamente los puntos en los que la divergencia de [math]\vec u [/math] es máxima, mínima y nula. La divergencia es una medida del cambio de volumen local debido al desplazamiento. ¿Se puede apreciar esto en la gráfica?

Representación de la divergencia

%Discretización de los parámetros de la superficie:
 % siendo 0.1 el "paso de muestreo"
 x=-2:0.1:2;
 y=0:0.1:4;
 %Creación del mallado
 [xx,yy]=meshgrid(x,y);
 %valor de la divergencia
 divu=-xx/10 +1/10;
 %representacion 
 surf(xx,yy,divu)
 shading interp 
 colorbar 
 title('divergencia') 
 xlabel('x')
 ylabel('y') 
 zlabel('z')
 %Dibujo de la superficie que representa la divergencia y de sus curvas de 
 %nivel visto cuando z tiende a infinito 
 figure
 pcolor(xx,yy,divu)
 shading interp 
 hold on 
 contour(xx,yy,divu,'k')
 xlabel('x')
 ylabel('y') 
 hold off 
 %Dibujo de la proyección sobre el plano YOZ de la superficie que representa
 %la divergencia 
 figure
  plot(yy,divu)
 title('proyección sobre el eje YOZ')
 xlabel('y')
 ylabel('z') 
 grid on
 maxi=max(max(divu))
 mini=min(min(divu))

Siendo la divergencia: -x/10 + 1/10. Hemos obtenido que la divergencia máxima dentro de la región de estudio es 0.3 y en x=-2. Mientras que la divergencia mínima tiene un valor de -0,1 en x=2. Finalmente, se anula la divergencia en x=1, como se aprecia en la gráfica la divergencia abarca desde 0 a 4 en la componente y. En la gráfica no apreciamos cambios de volumen.

8 Cálculo y representación del rotacional

Calcular |∇ × [math]\vec u [/math]| en todos los puntos del sólido y dibujarlo. ¿Qué puntos sufren un mayor 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{-x^2}{20} & \frac{y}{10} & 0 \end{matrix}\right| =\left( \frac{\partial 0}{\partial y} - \frac{ \partial \frac{y}{10}}{\partial z}\right)\vec i + \left(\frac{\partial \frac{-x^2}{20}}{\partial z} - \frac{\partial 0}{\partial x}\right)\vec j + \left(\frac{ \partial \frac{y}{10}}{\partial x} - \frac{ \partial \frac{-x^2}{20}}{\partial y}\right)\vec k =\left(0 - 0\right)\vec i + \left(0 - 0\right)\vec j + \left(0-0\right)\vec k =0 [/math]

Cálculamos el módulo del rotacional:

Como el rotacional es cero su módulo también es nulo. Por lo tanto no hemos podido representarlo.

9 Cálculo y representación de tensiones

Definamos ([math]\epsilon(\vec u) [/math]) = (∇[math]\vec u [/math] + ∇[math]\vec u [/math]t)/2, la parte simétrica del tensor gradiente de [math]\vec u [/math] conocido como tensor de deformaciones. En un medio elástico lineal, isótropo y homogéneo los desplazamientos permiten escribir el tensor de tensiones σij a través de la fórmula σ = λ∇ · [math]\vec u [/math] 1 + 2µ([math]\epsilon [/math]), donde λ y µ son los conocidos como coeficientes de Lamé que dependen de las propiedades elásticas de cada material. A pesar de que los desplazamientos son planos (es decir [math]\vec u [/math] no tiene componente en la dirección de [math]\vec k [/math]) las tensiones no tienen por qué ser planas y puede haber tensiones en la dirección ortogonal al plano de la placa. Tomando λ = µ = 1, dibujar las tensiones normales en la dirección que marca el eje v [math]\vec i [/math], es decir [math]\vec i [/math] · σ ·[math]\vec i [/math], las tensiones normales en la dirección que marca el eje [math]\vec j [/math], es decir [math]\vec j [/math] · σ · [math]\vec j [/math] y las correspondientes al eje [math]\vec k [/math], es decir [math]\vec k [/math] · σ · [math]\vec k [/math].

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{-x^2}{20} \vec i + \frac{y}{10} \vec j [/math]

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

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

De esta forma, el tensor de deformaciones quedará:

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

En un medio elástico lineal, isótropo y homogéneo, los desplazamientos permiten escribir el tensor de tensiones [math]\sigma_{ij}[/math] a través de la fórmula:

[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 dependen de las propiedades elásticas de cada material.

Tomando [math]\lambda = \mu = 1 [/math]. Calculamos 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 {-x}{10} + \frac{1}{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 {-x}{10} + \frac {1}{10} ) ( \vec i \otimes \vec i + \vec j \otimes \vec j + \vec k \otimes \vec k ) + 2 ( \frac {-x}{10} \vec i \otimes \vec i + \frac {1}{10}\vec j \otimes \vec j) = (\frac {-3x}{10} +\frac {1}{10} )\vec i \otimes \vec i + (\frac {3}{10} +\frac {-x}{10})\vec j \otimes \vec j + (\frac {-1}{10} +\frac {-x}{10}) \vec k \otimes \vec k) [/math]

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

9.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 {-3x+1}{10} [/math]

9.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 {3-x}{10} [/math]

9.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 {1-x}{10} [/math]

Tensiones normales en la dirección que marca el eje i, j y k

%Discretización de los parámetros de la superficie:
  % siendo 0.1 el "paso de muestreo"
  x=-2:0.1:2;
  y=0:0.1:4;
  %Creación del mallado
  [xx,yy]=meshgrid(x,y);
   
  %tensiones normales en las direcciones de los ejes coordenados sobre el
  %mallado 
 TNi=(-3.*xx+1)/10;
 TNj=(3-xx)/10;
 TNk=(1-xx)/10; 
 %representacion de las tensiones normales 
 figure
 subplot(1,3,1) 
 surf(xx,yy,TNi)
 axis image
 view(2) 
 colorbar
 xlabel('x')
 ylabel('y') 
 zlabel('z')
 title('Tensiones normales en la dirección del eje i')
  
 subplot(1,3,2) 
 surf(xx,yy,TNj)
 axis image
 view(2) 
 colorbar
 xlabel('x')
 ylabel('y') 
 zlabel('z')
 title('Tensiones normales en la dirección del eje j')
  
 subplot(1,3,3)
 surf(xx,yy,TNk)
 axis image
 view(2) 
 colorbar
 xlabel('x')
 ylabel('y') 
 zlabel('z')
 title('Tensiones normales en la dirección del eje k')

10 Calcular las tensiones respecto al plano ortogonal a

Calcular las tensiones tangenciales respecto al plano ortogonal a [math]\vec i [/math], es decir |σ · [math]\vec i [/math] − ( [math]\vec i [/math]· σ · [math]\vec i [/math]) [math]\vec i [/math]|.

como: [math] \vec i \cdot \sigma_{ij} \cdot \vec i = \frac {-3x+1}{10} [/math]

y ademas: [math] \sigma_{ij} \cdot \vec i =\frac {-3x+1}{10}\vec i [/math]

obtenemos: [math] \sigma_{ij} =0 [/math]

11 La tensión de Von Mises

11.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].

11.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 {-3x + 1}{10} & 0 & 0\\ 0 & \frac {3-x}{10} & 0\\ 0 & 0 & \frac {1-x}{10} \end{pmatrix} [/math]

La expresión de dichos autovalores, calculada analiticamente, sería igual a la anterior ya que sólo tenemos valores en la diagonal.

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.529150\).

%Discretización de los parámetros de la superficie:
  % siendo 0.1 el "paso de muestreo"
  x=-2:0.1:2;
  y=0:0.1:4;
  %Creación del mallado
  [xx,yy]=meshgrid(x,y);
  %calculo de los autovalores 
  A= inline('(-3.*x +1)/10','x','y');
  B= inline('(3-x)/10','x','y');
  C= inline('(1-x)/10','x','y');
  sigma=[];
  for i=1:length(x)
      for j=1:length(y)
          J=A(xx(i,j),yy(i,j));
          K=B(xx(i,j),yy(i,j));
          I=C(xx(i,j),yy(i,j));
          sigma=[J 0 0; 0 K 0 ; 0 0 I ];
          AU=eig(sigma);  %autovalores de sigma 
          VM=sqrt(((AU(1)-AU(2))^2+(AU(2)-AU(3))^2+(AU(3)-AU(1))^2)/2);
          zz(i,j)=VM;
      end 
  end 
 %dibujo 
 figure
 subplot(1,2,1)
 surf(xx,yy,zz)
 colorbar
 xlabel('x')
 ylabel('y') 
 zlabel('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=zz(1,1:length(zz));
 %Dibujo de la proyección sobre el plano XOZ de la superficie que genera la tensión de Von Mises
 plot(x,m,'b');
 maxzz=max(m);
 %Bucle que pinta los puntos donde la tensión de Von Mises es máxima sobre la proyección
 for k=1:length(m);
    if m(k)==maxzz
          plot(x(k),maxzz,'xr','markersize',10)
      end
 end
axis([-3 3 -0 2]);
 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',maxzz)

Representación de las tensiones de Von Mises

12 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.

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

El desplazamiento de la placa dado por el vector [math]\vec u(x,y) = -\frac{x^2}{20} \vec i + \frac{y}{10} \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{-3x-1}{10} }{\partial x} + \frac{\partial 0 }{\partial y} + \frac {\partial 0 }{\partial z} ) \cdot \vec i = {3}/{10}\vec i [/math]

[math] \vec{F_2} = - \frac {\partial \sigma_{j2}}{\partial y} \cdot \vec j = - ( \frac {\partial 0}{\partial x} + \frac {\partial \frac{3-x}{10}}{\partial y} + \frac {\partial 0 }{\partial z} ) \cdot \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{1-x}{10} }{\partial z} ) \cdot \vec k = 0 \vec k [/math]

12.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]:

Campo de vectores

%Discretización de los parámetros de la superficie:
% siendo 0.1 el "paso de muestreo"
 x=-2:0.1:2;      
 y=0:0.1:4;            

%Creación del mallado
 [xx,yy]=meshgrid(x,y); 

%componentes del campo de fuerzas
 fx=(3/10)+0.*xx;                 
 fy=0+0.*yy;   

%Dibujo del campo de fuerzas   
 figure(1)           
 quiver(xx,yy,fx,fy)     
 axis([-3,3,0,4])   
 xlabel('Eje X')
 ylabel('Eje Y')
 title('Campo de fuerzas')
 view(2)

12.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} = \frac{3}{10} \vec i [/math].Esto significa que las fuerzas solo actúan en la dirección i. Si observamos la placa después de la deformación vemos que sufre un alargamiento transversal que hace que adopte una mayor anchura.

13 Masa total de la placa

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

13.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_{-2}^{2}\int_{0}^{4} d(x,y) dydx=\int_{-2}^{2}\int_{0}^{4} 1+e^{\frac{-|x|}{(y+1)^{2}}} dydx= 29.474 [/math]

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

%Discretización de los parámetros de la superficie:
% siendo 0.1 el "paso de muestreo"
 x=-2:0.1:2;
 y=0:0.1:4;

%Creación del mallado
 [xx,yy]=meshgrid(x,y);  

%extremos integracion: 
 a=-2; b=2; c=0; d=4; 
 N1=length(x)-1; N2=N1;
 
%subintervalos 
 h1=(b-a)/N1; h2=(d-c)/N2; 
 u=a:h1:b; 
 v=c:h2:d; 
 [U,V]=meshgrid(u,v);

%densidad 
 densidad=1+exp(-abs(U)./((V+1).^2));

%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; 
 masa=h1*h2*(w2')*densidad*w1;
fprintf('El valor de la masa es %.4f', masa)