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

De MateWiki
Saltar a: navegación, buscar
Trabajo realizado por estudiantes
Título Visualización de campos escalares y vectoriales en elasticidad. Grupo 8C
Asignatura Teoría de Campos
Curso 2019-20
Autores Cruceta Fraile, Luis Miguel
Przybylski, Olivier Rafał
Rodelgo Gómez, Alejandro
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)\), 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 deformació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 una cierta función.

1 DISEÑO DE LA PLACA

Dibujamos un mallado que represente los puntos interiores del sólido, que al ser en 2D es una placa. Tomamos los ejes en el punto \((0, 0)\) del rectángulo \((x,y) ∈ [-2,2]×[0,4]\) y como paso de muestreo [math] h = \frac{1}{10} [/math] para las variables \((x,y)\).

Gráfica de la placa inicial
%Defimos los intervalos
%Lados
x = -2:1/10:2;
y = 0:1/10:4;
%Definimos la placa
[Mx,My] = meshgrid(x,y);
Mz = ones(size(Mx));
%Plotear placa
hold on
mesh (Mx,My,0*Mz)
title ('Placa inicial');
axis ([-3,3,-1,5])
view (0,90)
grid on
hold off


2 TEMPERATURA DE LA PLACA

La temperatura de la placa viene dada por la función [math] T(x,y) = (x+2)^2 + (y+2)^2 [/math] en grados centígrados.

2.1 Visualización del campo de temperaturas

A continuación ponemos el código de Matlab para ver el campo:

Temperatura de la placa
%Definimos la ecuación T(x,y)
T = (Mx+2).^2+(My+2).^2;
%Definimos curvas de nivel
figure 
hold on
mesh (Mx,My,0*Mz)
pcolor (Mx,My,T)
contour (Mx,My,T,'k')
shading flat
colorbar
hold off

2.2 Temperatura máxima y mínima

La temperatura máxima viene dada en el punto \((-2, 0)\) y es \(52,00 ºC\), y la temperatura mínima es \(4,00 ºC\) en el punto \((2, 4)\), ambos puntos en extremos opuestos en la placa.

3 CAMPO VECTORIAL GRADIENTE DE T \(( \nabla T)\)

Primero debemos hallar el gradiente de T \(( \nabla T)\) para expresarlo como campo vectorial, y después lo visualizamos con Matlab.

3.1 Cálculo del gradiente de T

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

3.2 Visualización del campo vectorial Gradiente de T

A continuación ponemos el código necesario para dicha visualización en Matlab.

Gradiente de la temperatura en la placa
%Definimos el gradiente
Tx = @(Mx,My) 2.*Mx+4;
Ty = @(Mx,My) 2.*My+4;
%Lo ponemos en ejes cartesianos
Ti = Tx(Mx,My);
Tj = Ty(Mx,My);
%Dibujamos el gradiente
figure
hold on
contour (Mx,My,T)
quiver (Mx,My,Ti,Tj)
hold off

4 CAMPO DE DESPLAZAMIENTOS [math]\vec u [/math]

Consideramos el campo [math]\vec u [/math] de desplazamientos sabiendo que:

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

Analicemos los dos puntos por separado.
En el primer punto, se nos indica que la función [math]\ \vec u(x,y) = \vec u_{x} (x,y) \vec i + \vec u_{y} (x,y) \vec j [/math], en su sentido más general, no sufre desplazamiento en la dirección [math]\vec j [/math] en el eje X, es decir, [math]\ \vec u (x,0) = \vec u_{x} (x,0) \vec i + 0 \vec j [/math].

En el segundo punto, se nos indica la divergencia del campo. De acuerdo a la definición de divergencia, [math]\ \nabla·\vec u =\frac{\partial \vec u_{x}}{\partial x} + \frac{\partial \vec u_{y}}{\partial y} [/math]. Es decir, en nuestro caso, [math]\ \frac{\partial \vec u_{x}}{\partial x} + \frac{\partial \vec u_{y}}{\partial y} = -\frac {x}{10} + \frac {1}{10}[/math]. Por tanto, de acuerdo al enunciado, [math]\ \vec u (x,y) = -\frac{x^2}{20} \vec i + f(x,y) \vec j[/math] . La divergencia de nuestro campo es, en este caso, [math]\ \nabla·\vec u =\frac{x}{10} + \frac{\partial f}{\partial y}[/math]. Es decir, [math]\ \nabla·\vec u =\frac{x}{10} + \frac{\partial f}{\partial y}=-\frac {x}{10} + \frac {1}{10} [/math]. Por tanto, [math]\ \frac{\partial f}{\partial y} = \frac {1}{10}[/math]. De aquí se puede obtener que [math]\ u_{y} (x,y) = \frac {1}{10}y + k(x)[/math].
De esta manera, nos queda que el campo de desplazamientos es [math]\ \vec u (x,y) = -\frac{x^2}{20} \vec i + \frac {1}{10}y + k(x) \vec j [/math]. Aplicando la condición de contorno [math]\ \vec u (x,y) = -\frac{x^2}{20} \vec i + (\frac {1}{10}*0 + k(x)) \vec j [/math] da que \(k(x)=0\) porque de otra manera el término [math]\ u_{y} (x,0) ≠ 0[/math]. Por ello, el campo de desplazamientos a considerar para ejercicios posteriores es: [math]\ u (x,y)=-\frac{1}{20} x^2 \vec i + \frac {1}{10} y \vec j [/math].

5 CAMPO DE VECTORES EN LOS PUNTOS DEL MALLADO

Disponiendo de la fórmula del campo de desplazamientos del apartado anterior, [math]\ u (x,y)=-\frac{1}{20} x^2 \vec i + \frac {1}{10} y \vec j [/math], podemos emplear un mallado de la placa sencillo donde en cada punto representemos con un vector el campo de desplazamientos.
El código de Matlab es el siguiente:

Gráfica de los desplazamientos de la placa
%Componentes en la dirección de i y de j del campo de desplazamientos: 
 ux=(-1/20).*(Mx).^2;
 uy=-1/10.*(My);
 figure
 %Dibujo de la placa: 
 mesh(Mx,My,0*Mx)
 hold on 
 %Dibujo del campo de desplazamientos:
 quiver(Mx,My,ux,uy)
 %Ajuste de ejes y ventana de visualización: 
 axis equal 
 axis([-3,3,-1,5]) 
 view(2) 
 hold off

Como se puede apreciar en la imagen, el campo no presenta desplazamientos en el sentido del vector [math]\vec j [/math] en el eje \(y=0\), tal y como se indicaba en el enunciado.

6 PLACA ANTES Y DESPUÉS DEL DESPLAZAMIENTO DADO POR [math]\vec u [/math]

Gracias al comando subplot, podemos dibujar la situación de la placa antes y después del desplazamiento. La representación de la placa antes del desplazamiento es trivial, ya que basta con realizar un mallado de la placa y representarlo en el plano de coordenadas. Para obtener la placa una vez aplicado el campo de desplazamientos, se toma la placa en posición inicial y se aplica a cada punto del mallado su desplazamiento correspondiente, de manera que [math]P(x,y) = P_{0}(x,y) + \vec u(x,y)[/math], siendo [math]\vec u(x,y) )=-\frac{1}{20} x^2 \vec i + \frac {1}{10} y \vec j [/math].
El código de Matlab utilizado para obtener dicho campo es el siguiente:

Desplazamientos de la placa antes y después
%Posición de la placa deformada:
rx=Mx-(1/20).*[Mx.^2];
ry=My+(1/10).*My; 
%Representación de la placa antes del desplazamiento: 
subplot(1,2,1)
surf(Mx,My,0*Mx)
title('Posición original')
axis equal
axis([-3,3,-1,5])
view(2) 
xlabel('Eje X')
ylabel('Eje Y') 
zlabel('Eje Z')
%Representación de la placa después del desplazamiento:
subplot(1,2,2)
surf(rx,ry,0*rx)
title('Después del desplazamiento')
axis equal
axis([-3,3,-1,5])
xlabel('Eje X')
ylabel('Eje Y') 
zlabel('Eje Z') 
view(2)

7 DIVERGENCIA DEL CAMPO DE DESPLAZAMIENTOS

A continuación calcularemos la divergencia del campo de desplazamientos y determinaremos cuando la divergencia será máxima, mínima y nula.

7.1 Cálculo de la divergencia

Para el cálculo de la divergencia en coordenadas cartesianas, se procede a derivar el campo de desplazamientos hallado anteriormente:[math]\ \vec u (x, y)= \frac{-x^2}{20} \vec i + \frac {y}{10}\vec j[/math]

[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{-x^2}{20})}{\partial x} + \frac{\partial (\frac{y}{10})}{\partial y} + \frac{\partial 0}{\partial z} = - \frac{x}{10} + \frac{1}{10} [/math]

7.2 Visualización de la divergencia

Mediante el siguiente código se puede representar la divergencia:

Variación de la divergencia y su máximo, mínimo y punto nulo.
%Ponemos la divergencia
Mz = -(1/10)*Mx+1/10;
z = -x/10+1/10;
%Dibujamos la divergencia
figure
subplot(1,2,1)
pcolor(Mx,My,Mz)
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
subplot(1,2,2)
plot(y,z)
hold on
title('Máximo, mínimo y nulo')
xlabel('Eje Y')
ylabel('Eje Z')
plot(0,0.3,'xr','markersize',10)
plot(3,0,'xg','markersize',10)
plot(4,-0.1,'xb','markersize',10)
legend('Diverencia','Máximo','Nulo','Mínimo');
hold off

Como se puede observar en la representación, la divergencia adopta la forma de plano inclinado en el cual el máximo queda determinado por la intersección con el plano \(z=0.3\), el mínimo por el plano \(z=-0.1\) y el nulo por el plano \(z=0\). El cambio de volumen se puede apreciar gracias al color, según nos desplazamos en la dirección del eje X o Y.

8 ROTACIONAL DEL CAMPO DE DESPLAZAMIENTOS

Debemos calcular el rotacional en todos los puntos del sólido e indicar qué puntos sufren un mayor 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 = \vec 0 [/math]
Por tanto, se comprueba que el rotacional en todos los puntos de la placa es nulo, lo que implica que se trata de un campo vectorial conservativo.

9 TENSIONES NORMALES

9.1 Tensor de deformaciones [math]\epsilon (\vec u) [/math]

Para calcular este tensor, que se define como la parte simétrica del tensor gradiente de [math]\ \vec u, \nabla \vec u[/math], debemos sacar el gradiente de [math]\ \vec u[/math] y su traspuesto [math]\ \vec u^t[/math].

Tomando como campo de desplazamientos [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)\vec i + (\frac{1}{10}\vec j)\vec j = -\frac{x}{10}\vec i \otimes \vec i + \frac{1}{10}\vec j \otimes \vec j[/math] Como los tensores diada tienen el mismo vector dos a dos, el gradiente y su traspuesto es el mismo. Por lo tanto, el tensor [math]\epsilon (\vec u) [/math] coincide con el [math]\ \nabla \vec u[/math].

9.2 Tensor tensiones [math] \sigma_{ij}[/math]

Una vez hallado el tensor de deformaciones, se puede hallar el tensor de tensiones mediante la expresión [math]\sigma_{ij} = \lambda \nabla \cdot \vec u \boldsymbol {1} + 2 \mu \epsilon [/math]. La divergencia [math]\ \nabla · \vec u[/math] ha sido hallada previamente: [math] \nabla·\vec u =-\frac {x}{10} + \frac {1}{10}[/math].

[math] \sigma_{ij} = \lambda \nabla · \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 {1-3x}{10} \vec i \otimes \vec i - \frac {3-x}{10} \vec j \otimes \vec j + \frac {1-x}{10} \vec k \otimes \vec k) [/math]


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

9.3 Cálculo de tensiones normales

Para calcular las tensiones normales en las direcciones de los ejes, basta con tomar los elementos de la diagonal:
[math] \vec i · \sigma_{ij} · \vec i = \frac {1-3x}{10} [/math] en la dirección que marca [math] \vec i[/math].
[math] \vec j · \sigma_{ij} · \vec j = \frac {3-x}{10} [/math] en la dirección que marca [math] \vec j[/math].
[math] \vec k · \sigma_{ij} · \vec k = \frac {1-x}{10} [/math] en la dirección que marca [math] \vec k[/math].

9.4 Visualización de las tensiones normales

Estas tensiones normales se pueden dibujar en Matlab mediante el siguiente código: 

TNi=(1-3*Mx)/10;
TNj=(3-Mx)/10;
TNk=(1-Mx)/10; 
figure
subplot(1,3,1) 
 pcolor(Mx,My,TNi)
 colorbar
 xlabel('Eje X')
 ylabel('Eje Y') 
 zlabel('Eje Z')
 title('Tensiones normales en la dirección del eje i')  
subplot(1,3,2) 
 pcolor(Mx,My,TNj)
 colorbar
 xlabel('Eje X')
 ylabel('Eje Y') 
 zlabel('Eje Z')
 title('Tensiones normales en la dirección del eje j')
subplot(1,3,3)
 pcolor(Mx,My,TNk)
 colorbar
 xlabel('Eje X')
 ylabel('Eje Y') 
 zlabel('Eje Z')
 title('Tensiones normales en la dirección del eje k')
Tensiones normales en las direcciones de los ejes

10 TENSIONES TANGENCIALES RESPECTO A [math]\ \vec i [/math]

En este apartado calculamos la tensión tangencial respecto al plano ortogonal a [math]\ \vec i [/math], que viene dada por [math]\ \left| \sigma ·\vec i -(\vec i·\sigma·\vec i)\vec i \right| [/math], siendo [math]\ \sigma_{ij} = \lambda \nabla · \vec u \boldsymbol {1} + 2 \mu \epsilon[/math].
Para ello, se procede matricialmente según la expresión anterior:
[math]\ \left| \sigma ·\vec i -(\vec i·\sigma·\vec i)\vec i \right| = \left| \begin{pmatrix} \frac {1-3x}{10} & 0 & 0\\ 0 & \frac {3-x}{10} & 0\\ 0 & 0 & \frac {1-x}{10} \end{pmatrix} · \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix} - (\frac{1-3X}{10}) · \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix} \right|= \left| \begin{pmatrix} \frac{1-3X}{10} \\ 0 \\ 0 \end{pmatrix} - \begin{pmatrix} \frac{1-3X}{10} \\ 0 \\ 0 \end{pmatrix}\right| = \vec 0[/math]
De acuerdo con el resultado obtenido, podemos afirmar que no hay tensiones tangenciales respecto al plano ortogonal a [math]\ \vec i[/math], ni a ningún plano ortogonal a [math]\ \vec j[/math] o [math]\ \vec k[/math].

11 TENSIÓN DE VON MISES

La tensión de Von Mises se define 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]. Esta tensión de Von Mises nos indica cuando un material inicia un comportamiento plástico. Es una magnitud escalar.

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

[math]\sigma_{1}, \sigma_{2}[/math] y [math]\sigma_{3}[/math] son los autovalores de [math]\sigma[/math]:
[math] \sigma_{ij} = \begin{pmatrix} \frac {1-3x}{10} & 0 & 0\\ 0 & \frac {3-x}{10} & 0\\ 0 & 0 & \frac {1-x}{10} \end{pmatrix} [/math]
Al ser estas tensiones en la diagonal principal, estas son los autovalores. 

%Cálculo de los autovalores 
TX=@(Mx,My)[My/10*3 -Mx/10 0];
TY=@(Mx,My)[-Mx/10 My/10 0];
TZ=@(Mx,My)[0 0 My/10];
%Autovalores de la tensión de Von Mises
for i=1:length(x)
    for j=1:length(y)
    %Valores para cada componente del mallado para la función de tensión de Von Mises    
    U=TX(Mx(i,j),My(i,j));
    V=TY(Mx(i,j),My(i,j));
    W=TZ(Mx(i,j),My(i,j));
    %Vector con los valores obtenidos
    vec=[U; V; W]; 
    %Obtención de los autovalores de la tensión de Von Mises 
    auto=eig(vec);
    %Cálculo de Von Mises para cada componente
    VM=sqrt(((auto(1)-auto(2))^2+(auto(2)-auto(3))^2+((auto(3)-auto(1))^2))*1/2);
    Mz(i,j)=VM;
    end
end
%Representamos Von Mises
figure
subplot(1,2,1)
surf(Mx,My,Mz)
colorbar
xlabel('Eje X')
ylabel('Eje Y')
zlabel('Eje Z')
title('Tensión de Von Mises')
%El máximo valor de la tensión de Von Mises
subplot(1,2,2)
hold on
m=Mz(1,1:length(Mz));
%Proyección sobre el plano XOZ de la superficie dse la tensión de Von Mises
plot(x,m,'b');
maxz=max(m);
%Dibujar el máximo 
for k=1:length(m)
   if m(k)==maxz
         plot(x(k),maxz,'xr','markersize',10)
     end
end
axis([-2.5 2.5 0 0.5]);
xlabel('Eje X')
ylabel('Eje Z')
legend('Proyección sobre el plano XOZ','Punto máximo')
hold off


Tensiones de Von Mises

Según el gráfico, la máxima tensión de Von Mises se alcanza en [math]\ x=2 [/math] y [math]\ x=-2 [/math]

12 CAMPO DE FUERZAS [math]\ \vec F [/math] SOBRE LA PLACA

El desplazamiento de la placa expresado mediante el vector [math]\vec u(x, y) = -\frac{x^2}{20} \vec i + f(y) \vec j, [/math] está causado por un campo de fuerzas [math]\ \vec F [/math] que actúa en el interior de la placa.

12.1 Cálculo del campo de fuerzas [math]\ \vec F [/math]

Para calcularlo, vamos a utilizar la ecuación de la elasticidad lineal: [math]\ \vec F = \nabla ·\sigma_{ij}= -\frac {\partial \sigma_{ij}}{\partial x_j} · \vec e_i[/math]
De esta manera obtenemos las siguientes componentes aplicando la divergencia a los vectores fila de la matriz [math]\ \sigma_{ij}[/math].
[math]\ \vec F_i = -\frac {\partial \sigma_{j1}}{\partial x_j} · \vec e_i = -\left( \frac {\partial \frac {1-3x}{10}}{\partial x} + \frac {\partial 0}{\partial y} + \frac {\partial 0}{\partial z} \right) · \vec i = \frac {-3}{10} \vec i[/math]

[math]\ \vec F_j = -\frac {\partial \sigma_{j1}}{\partial x_j} · \vec e_j = -\left( \frac {\partial 0}{\partial x} + \frac {\partial \frac {3-x}{10}}{\partial y} + \frac {\partial 0}{\partial z} \right) · \vec j = 0 \vec j[/math]

[math]\ \vec F_k = -\frac {\partial \sigma_{j1}}{\partial x_j} · \vec e_k = -\left( \frac {\partial 0}{\partial x} + \frac {\partial 0}{\partial y} + \frac {\partial \frac {1-x}{10}}{\partial z} \right) · \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]:

Fx=inline('Mx.*(-3./10)','Mx');% Definimos el campo de fuerzas F y sus componentes
Fy=0.*My;
U=Fx(Mx);
V=Fy;
quiver(Mx,My,U,V)                          % Dibujamos el campo de fuerzas
title('CAMPO DE FUERZAS')                % Título del gráfico
xlabel('EJE X')                          % Título del eje x
ylabel('EJE Y')                          % Título del eje y
axis([-3,3,-1,5])% Establecimiento de valores mínimos y máximos en los ejes representados
Campo de fuerza [math]\ \vec F [/math] sobre la placa

12.3 Interpretación del campo de fuerzas [math]\ \vec F [/math]

Esto significa que las fuerzas solo actúan en la dirección de [math]\ \vec i [/math]. Si observamos la placa después de la deformación vemos que sufre un estrechamiento transversal que hace que adopte una forma más alargada. Esto se debe a que las fuerzas que actúan en sentido contrario a los dos lados de la placa, la comprime, dando lugar a la deformación citada anteriormente.


13 MASA TOTAL DE LA PLACA

Para este ejercicio hemos decidido aplicar el método del trapecio extendido a dos dimensiones. Primero, hemos calculado la función densidad, y la hemos representado en [math]R^3[/math] así como su proyección sobre la placa para ver cómo varía a lo largo y ancho de la placa.
La densidad de la placa es: [math] d(x,y) = 1 + e^{-|x|/(y+1)^2} [/math]

13.1 Distribución de la densidad

Al disponer ya de la función densidad, bastó con realizar un mallado discreto y una interpolación entre los puntos del mallado para obtener la superficie correspondiente a la función Densidad, que atribuye a cada punto la altura igual a su densidad.
El código Matlab utilizado fue el siguiente:

Función densidad de la placa
%Función densidad:
d=inline('(1+exp(-abs(x)./(y+1).^2))','x','y');
%Representación de la densidad en la 3ª dimensión
Z=d(Mx,My);
subplot(1,2,1)
surf(Mx,My,Z)
axis equal
shading flat
title('Superficie de la función Densidad');
xlabel('Eje X')
ylabel('Eje Y')
zlabel('Eje Z')
%Proyección de la función densidad sobre la placa:
subplot(1,2,2)
pcolor(X,Y,Z)
shading flat
colorbar
axis equal
axis([-3 3 -1 5]);
title('Proyección de la función Densidad sobre la placa');
xlabel('Eje X')
ylabel('Eje Y')


13.2 Cálculo de la masa

En cuanto al cálculo de la masa, primero tuvimos que hacer una disección de la placa en elementos uniformes. En este caso, para la facilitación del computado, se optó por cuadrados. Por el método del trapecio extendido a las dos dimensiones, se van calculando las masas de cada elemento como si de volúmenes se tratase, situándose la densidad en la tercera dimensión. En los extremos de la placa es necesario ponderar los vértices puesto que sólo forman parte de uno o a lo sumo dos elementos a la vez, mientras que los vértices interiores forman parte de cuatro elementos simultáneamente.
Por ello, el código Matlab utilizado para obtener el resultado correspondiente a la masa es el siguiente:

Visualización del mallado de vértices ponderado
%Mallado discreto de la placa:
pm=1/10;
x=-2:pm:2;
y=0:pm:4;
%Mallado de vértices:
N1=length(x)-1; 
N2=N1;
%Límites de integración:
a=-2; 
b=2; 
c=0; 
d=4;
%Intervalos de paso para la creación del mallado:
h1=(b-a)/N1; 
h2=(d-c)/N2;
X=a:h1:b; 
Y=c:h2:d;
%Creación del mallado a partir de los intervalos de paso:
[XX,YY]=meshgrid(X,Y);   
%Función densidad aplicada a cada vértice:
D=(1+exp((-abs(XX))./(YY+1).^2)); 
surf(XX,YY,D)
axis equal
%Coeficientes de ponderación de vértices:
P1=ones(N1+1,1);                 
P1(1)=1/2; 
P1(N1+1)=1/2;
P2=ones(N2+1,1);                 
P2(1)=1/2;  
P2(N2+1)=1/2;
%Cálculo de la masa como producto matricial de volúmenes:
M=h1*h2*rot90(P2)*D*P1;

De aquí se obtiene que la placa tiene 29.4742 kg de masa.

Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Visualización_de_campos_escalares_y_vectoriales_en_elasticidad_(Grupo_C8)&oldid=49284»