Placa Plana (Grupo 26)

Trabajo realizado por estudiantes
Título	Placa plana. Grupo 26
Asignatura	Teoría de Campos
Curso	2024-25
Autores	Jorge Muñoz Jimenez Eva Aragon Peña Armando de Tomas Fernandez Antonio Gurría Casas Daniel Galarza Polo
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

Una placa rectangular plana en la región [math](x, y) ∈ [-1 , 1] × [0,f(x)][/math] viene definida en dimension 2. La función [math]f(x)[/math] es la siguiente: [math]f(x) = 2 + x^2 [/math] Supondremos que están definidas dos cantidades físicas.

• La Temperatura
• Los Desplazamientos

La temperatura [math]T(x, y)[/math] viene dada por la ecuación:

[math]T(x, y)=(1-x^4)+(\frac{1}{2}-y)[/math]

Los desplazamientos [math]u(x, y)[/math] producidos por la acción de una fuerza determinada.

Al definir el vector de posición de los puntos de la placa antes de que se produzca cualquier deformación [math]\vec{r_{0}}(x, y)=x\vec{i}+y\vec{j}[/math] , la posición de cada punto [math](x, y)[/math] de la placa después de la deformación vendrá dada mediante la ecuación [math]\vec{r}(x, y)=\vec{r_{0}}(x, y)+\vec{u}(x, y)[/math]. Debido a la aplicación de fuerza sobre la placa, esta sufre un desplazamiento de los puntos, este desplazamiento viene determinado por el vector [math]\vec{u}(x, y) = \frac{xy \vec{i} - yx^2 \vec{j}}{10}[/math].

Haciendo uso del programa Matlab podremos determinar las gráficas de las operaciones calculadas en los siguientes apartados.

1 Mallado

A continuación podemos ver la representación del mallado que contiene a los puntos interiores del sólido realizado con MatLab.

Tomamos como ejes \((x,y) ∈ [−2,2] × [0,3]\) y un paso de muestreo, es decir, el intervalo entre punto y punto, [math]h=\frac{1}{10}[/math] para las variables [math]x[/math] e [math]y[/math].

derecha

% configuración de los ejes
axis equal 
axis([-2,2,0,3])
view(2)
% APARTADO 1- Malla
h=0.1; %paso de muestreo
%definicion de las variables
x=(-1:h:1);
y=(0:h:1);
[mx,my]=meshgrid(x,y);
%Deformacion parabolica de la malla
yy=my.*(mx.^2+2);
%mallado
hold on 
mesh(mx,yy,0.*mx);

2 Temperatura

La temperatura viene dada por la siguiente expresión [math] T = (1 - x^4) (\frac {1}{2}- y)[/math], que depende únicamente de x e y. Las curvas de nivel de este campo escalar vienen representadas en la gráfica en color azul marino.

El gradiente de la temperatura ([math]\nabla T[/math]) se obtiene derivando parcialmente la expresión de T respecto de x e y, y posteriormente, sumándolas. Este se representa en la gráfica mediante flechas de color azul celeste y nos indica la dirección máxima de crecimiento. Como se puede observar en la imagen, las flechas de dicho campo vectorial son perpendiculares a las curvas de nivel de nuestra placa plana.

Tras ejecutar el programa para determinar el punto en el que la temperatura es máxima, se puede determinar que el valor máximo de T es 0,5000 en el punto [math](x , y) = (0.000 , 0.000)[/math].

derecha

%Configuracion de los ejes
axis equal
axis([-2,2,0,3])
view(2)
%Apartado 1 - Malla
h=0.1; %Paso de muestreo
%Definición de variables
x=(-1:h:1);
y=(0:h:1);
[mx,my]=meshgrid(x,y);
%Deformacion parabólica de la malla
yy=my.*(mx.^2+2);
hold on
%Apartado 2 - Temperatura y Gradiente
%Dibujamos las curvas de nivel de la temperatura
T=(1-mx.^4).*(0.5-yy);
contour(mx,yy,T,50,´b´);
%Dibujamos el gradiente
[gx,gy]=gradient(T,h,h);
quiver(mx,yy,gx,gy);
hold off

3 Ley de Fourier

La Ley de Fourier concluye que tras estudiar el flujo de calor entre dos cuerpos, se determina que la diferencia de temperatura entre ambos es directamente proporcional, solo podrá ir del cuerpo mas caliente al cuerpo mas frio, lo que significa que ira en una sola dirección. Para la implementación de esta ley se deben cumplir tres condiciones.

• Sistema isotropo
• Gradiente de temperatura pequeño
• No hay transferencia de calor por convección ni radiación
  

La fórmula de la Ley de Fourier es: [math]\vec{Q}= −κ∇T[/math]

Calculamos el valor de [math]\vec{Q}[/math], donde la conductividades térmica [math]κ[/math] tendrá valor de 1.

A partir de Matlab dibujamos el campo vectorial.

Como en el apartado anterior, las flechas son perpendiculares a las curvas de nivel. Sin embargo, en este caso en dirección opuesta.

derecha

%Configuracion de los ejes
axis equal
axis([-2,2,0,3])
view(2)
%Apartado 1 - Malla
h=0.1; %Paso de muestreo
%Definición de variables
x=(-1:h:1);
y=(0:h:1);
[mx,my]=meshgrid(x,y);
%Deformacion parabólica de la malla
yy=my.*(mx.^2+2);
%Mallado
hold on
%Dibujamos las curvas de nivel de la temperatura
T= (1-mx.^4).*(0.5-yy);
%Dibujamos el gradiente 
[gx,gy]=gradient(T,h,h);
%Calculamos y representamos Q
k=-1;
qx=k.*gx;
qy=k.*gy;
quiver(mx,yy,qx,qy);
hold off

4 Variación de Temperatura

Para estudiar la variación del campo escalar [math]T[/math], haremos uso del gradiente [math]\nabla T[/math] . Se puede observar, por la propia definición del [math]\nabla[/math], que dicho campo vectorial es perpendicular a las curvas de nivel de nuestra placa plana. A medida que nos acercamos a temperaturas más elevadas, el módulo de [math]\vec\nabla[/math] en cada punto va siendo mayor.

El punto de variación máxima de la temperatura se alcanza cuando el módulo del gradiente sea máximo. En nuestra placa plana, ese punto se encuentra en el (-1,1). La dirección de la variación de temperatura viene dada por la flecha roja representada sobre la placa.

El código Matlab es el siguiente:

derecha

%Hay que dibujarlo como un solido
mod=sqrt(gx.^2+gy.^2);
%Configuración de los ejes
axis equal
axis([-2,2,0,3])
view(2)
%Dibujamos la Malla
h=0.1;
%Definicion de variables
x=(-1:h:1);
y=(0:h:1);
[mx,my]=meshgrid(x,y);
%Deformacion parabolica de la malla
yy=my.*(mx.^2+2);
%Calculo del gradiente de la deformación 
[gy,gx]=gradient(yy,h,h); %Gradiente en x,y (orden invertido)
%Calculo de la magnitud del gradiente
mod_grad=sqrt(gx.^2+gy.^2);
%Encontrar el indice dle punto de mayor magnitud del gradiente
[max_grad_value, idx] = max(mod_grad(:)); %Mayor valor del gradiente
[row,col] = inds2sub(size(mod_grad),idx); %Indices del punto
%Coordenadas del punto con mayor gradiente
x_max = mx(row, col);
y_max = my(row,col);
%Direccion del gradiente en el punto de mayor magnitud
gx_max = gx(row, col);
gy_max = gy(row, col); 
%Normalizacion del vector de dirección 
magnitude = sqrt(gx_max^2 + gy_max^2);
direction_vector = [gx_max / magnitude, gy_max / magnitude];
%Graficar la malla
hold on 
mesh(mx, yy, zeros(size(mx))); %Malla
%Dibujar el punto rojo en el lugar donde la magnitud del gradiente es maxima
plot3(x_max, y_max, 0, 'ro', 'MarkerSize', 10, 'Linewidth', 2);

5 Campo de desplazamientos

El campo vectorial de desplazamiento viene dado por [math]\vec{u}(x, y) = \frac{xy \vec{i} - yx^2 \vec{j}}{10}[/math]. Este campo se representa mediante las flechas rojas.

El único punto de la placa que no se ve afectado por dicho desplazamiento es el origen de coordenadas (0,0).

derecha

%Hay que dibujarlo como un solido
mod=sqrt(gx.^2+gy.^2);
%Configuración de los ejes
axis equal
axis([-2,2,0,3])
view(2)
%Dibujamos la Malla
h=0.1;
%Definicion de variables
x=(-1:h:1);
y=(0:h:1);
[mx,my]=meshgrid(x,y);
%Deformacion parabolica de la malla
yy=my.*(mx.^2+2);
%Campo de desplazamiento
Ux=(mx.*yy)./10;
Uy=(-yy.*mx.^2)/10);
quiver(mx,yy,Ux,Uy,'r');

6 Desplazamiento de los vectores

La placa plana sufre el desplazamiento de los vectores, para poder contrastar el desplazamiento ocurrido, representaremos en dos gráficas para poder observar el cambio que sucede.

Los colores representan de más cálido a más frío el desplazamiento, siendo los puntos representados en colores fríos los que más desplazamiento sufren.

center

%Configuracion de los ejes
h=0.1; %paso de muestreo
%Definicion de variables
x=(-1:h:1);
y=(0:h:1);
[mx,my]=meshgrid(x,y);
%Deformacion parabolica de la malla
yy=my.*(mx.^2+2);
%Campo de desplazamiento
Ux=(mx.*yy)./10;
Uy=(-yy.*mx.^2)/10;
%Placa pre-desplazamiento
subplot(1,3,1);
surf(mx,yy,0.*mx);
axis equal
axis([-2,2,0,3])
view(2)
title('Original')
%Placa post-desplazamiento
subplot(1,3,2);
surf(mx,yy,Ux,Uy);
axis equal
axis([-2,2,0,3])
view(2)
title('Desplazado')
%Comparacion
subplot(1,3,3);
hold on
surf(mx,yy,0.*mx);
surf(mx,yy,Ux,Uy);
axis equal
axis([-2,2,0,3])
view(2)
title('Comparativa')
hold off

7 Divergencia

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

FALTA PONER BONITO ESTE PÁRRAFO Tanto el máximo como el mínimo se alcanzan donde el módulo de la divergencia de U sea máxima o mínima respectivamente. El punto donde el módulo de la divergencia de U es nula, será donde se anula la divergencia.

. La divegencia es nula en el punto (0,0). . La divergencia es máxima en el punto (1,1) y con un valor de -1/10- . La divergencia es mínima en el punto (-1,0) y con un valor de -1/5.

Se representa la gráfica en un mapa de colores (barra de color a la derecha). En base a ello, interpretamos el cambio de volumen local debido al desplazamiento:

1. Colores positivos (amarillos/verde claro): Indican las regiones en las que la divergencia es positiva, lo que implica que el volumen local está aumentando debido al desplazamiento.
2. Colores negativos (azules): En estos colores se indican las regiones donde la divergencia es negativa, y por tanto el volumen está disminuyendo en consecuencia del desplazamiento.
3. Zonas de transición (cerca de 0): Indican los puntos donde el cambio de volumen es mínimo o nulo debido al desplazamiento.

derecha

%Paso de muestreo
h=0.1;
%Definicion de variables
x=(-1:h:1);
y=(0:h:1);
[mx,my]=meshgrid(x,y);
%Deformacion parabolica de la malla
yy=my.*(mx.^2+2);
mesh(mx,yy,0.*mx);
%Definimos variable U
Ux=(mx.*yy)./10);
Uy=(-yy.*mx.^2)/10;
%Divergencia de U
div=divergence(mx,yy,Ux,Uy);
maxdiv=max(max(div));
mindiv=min(min(div));
surf(mx,yy,div);
%Configuración de los ejes
axis equal
axis([-2,2,0,3])

8 Rotacional

El rotacional se utiliza para ver cómo afecta el campo de fuerzas a un objeto en su movimiento, pero, en especial, a su movimiento angular, es decir, el rotacional nos indica cómo girará un cuerpo en un determinado punto afectado por el campo de fuerzas. Calcularemos el rotacional del campo como se expone abajo.

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

[math]|∇ × \vec{u}| = \begin{bmatrix} \vec{i} & \vec{j} \\ \frac{\partial}{\partial x} & \frac{\partial}{\partial y} \\ \vec{u_1} & \vec{u_2} \end{bmatrix} [/math] => [math] \begin{bmatrix} \vec{i} & \vec{j} & \vec{k} \\ \frac{y}{10} & \frac{-2yx}{10} & 0 \\ \frac{xy}{10} & \frac{-yx^2}{10} & 0 \end{bmatrix}[/math] [math] = (\frac{y}{10} × \frac{-yx^2}{10}) \vec{k} - (\frac{xy}{10} × \frac{-2yx}{10}) \vec{k} = [/math]

[math] = (\frac{-x^2y^2}{100} - \frac{ 2x^2y^2}{100}) \vec{k} = (\frac{-3x^2y^2}{100}) \vec{k} [/math]

Se observa que el resultado final se encuentra en el vector [math]\vec{k} [/math], esto se debe a que el rotacional sera perpendicular tanto a la componente [math]\vec{i} [/math] y [math]\vec{j} [/math] respectivamente.

A continuación, representaremos el rotacional gráficamente, para ello hacemos uso del programa Matlab.

derecha

%Paso de muestreo
h=0.1;
%Defenicion de variables
x=(-1:h:1);
y=(0:h:1);
[mx,my]=meshgrid(x,y);
%Deformacion parabolica de la malla
yy=my.*(mx.^2+2);
mesh(mx,yy,0.*mx);
%Definimos variable U
Ux=(mx.*yy)./10;
Uy=(-yy.*mx.^2)/10;
%Rotacional de U 
[rot,ang]=curl(mx,yy,Ux,Uy);
surf(mx,yy,rot);
maxrot=max(max(rot));
%Configuracion de los ejes
axis equal
axis([-2,2,0,3])

Tras la ejecución del programa, podemos determinar que los puntos donde el rotacional es máximo son los siguientes: [math] x = -1,00 ; y = 3,00[/math] El resultado final es 0,82.

La deformación parabolica con el puto de mayor gradiente de temperatura, se observa el la siguiente imagen.

center

9 Tension Deformación

Definiendo [math]ε(\vec u)=\frac {\nabla \vec u + \nabla \vec u ^t}{2}[/math], 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 a través de la fórmula [math]σ=λ\nabla·\vec u 1+2με[/math].

Calculamos el gradiente de [math] ∇\vec{u}[/math]

[math] ∇\vec{u} = \begin{pmatrix}\frac{\partial u_1}{\partial x} & \frac{\partial u_1}{\partial y} & \frac{\partial u_1}{\partial z} \\\frac{\partial u_2}{\partial x} & \frac{\partial u_2}{\partial y} & \frac{\partial u_2}{\partial z} \\\frac{\partial u_3}{\partial x} & \frac{\partial u_3}{\partial y} & \frac{\partial u_3}{\partial z}\end{pmatrix} = \begin{pmatrix}\frac{y}{10} & \frac{x}{10} & 0 \\-2x\frac{y}{10} & -\frac{x^2}{10} & 0 \\0 & 0 & 0\end{pmatrix}[/math]

A partir de esto obtendremos también su traspuesta, que quedara de la siguiente manera [math]∇\vec{u^t}= \begin{pmatrix}\frac{y}{10} & -2x\frac{y}{10} & 0 \\\frac{x}{10} & -\frac{x^2}{10} & 0 \\0 & 0 & 0\end{pmatrix} [/math]

[math] є= \frac{∇\vec{u} + ∇\vec{u^t}}{2} = \begin{pmatrix} \frac{y}{5} & (-y+0,5)\frac{x}{5} & 0 \\ (-y+0,5)\frac{x}{5} & -\frac{x^2}{10} & 0 \\ 0 & 0 & 0 \end{pmatrix} [/math]

A continuación, hallamos el valor de σ

[math] σ=\begin{pmatrix} \frac{1}{10}(y-x^2) & 0 & 0 \\ 0 & \frac{1}{10}(y-x^2) & 0 \\ 0 & 0 & \frac{1}{10}(y-x^2) \end{pmatrix} [/math] + [math] \begin{pmatrix} \frac{2y}{10} & (-2y+1)\frac{x}{10} & 0 \\ (-2y+1)\frac{x}{10} & \frac{-2x^2}{10} & 0 \\ 0 & 0 & 0 \end{pmatrix} [/math]

Por lo tanto:

[math]σ=\begin{pmatrix} (3y-x^2)\frac{1}{10} & (-2y+1)\frac{x}{10} & 0 \\ (-2y+1)\frac{x}{10} & (y-3x^2)\frac{1}{10} & 0 \\ 0 & 0 & (y-x^2)\frac{1}{10} \end{pmatrix} [/math]

Todos estos cálculos están representados en la siguiente gráfica, realizada a partir del programa Matlab.

derecha

%Paso de muestreo
h=0.1;
%Defenicion de variables
x=(-1:h:1);
y=(0:h:1);
[mx,my]=meshgrid(x,y);
%Deformacion parabolica de la malla
yy=my.*(mx.^2+2);
mesh(mx,yy,0.*mx);
%Definimos variable U
Ux=(mx.*yy)./10;
Uy=(-yy.*mx.^2)/10;
%Trazado i*tension*i.
t1=(1/10).*(3.*yy-mx.^2);
subplot(1,3,1);
surf(mx,yy,t1);
title('I*tension*i);
%Trazado j*tension*j
t2=(1/10).*(yy-3.*mx.^2)
subplot(1,3,2);
surf(mx,yy,t2);
title('j*tension*j');
%Trazado k*tension*k
t3=(1/10).*(yy-mx.^2);
subplot(1,3,3);
surf(mx,yy,t3);
title('k*tension*k);

10 Tensiones tangenciales

La tension tangencial respecto al plano ortogonal [math] \vec{i} [/math]. Los cálculos son. [math] |σ·\vec{i} - ( \vec{i}·σ·\vec{i})\vec{i}|= | \begin{pmatrix} \frac{1}{10}( 3y-x^2 ) & (\frac{x}{10} - \frac{2xy}{10}) & 0 \\ \frac{x}{10}- \frac{2xy}{10} & - \frac{1}{10} (y - 3x^2) & 0 \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix} - (\frac {1}{10} (3y - x^2)) \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix} | =[/math]

[math] |\begin{pmatrix} \frac{1}{10}(3y - x^2) \\ (\frac{x}{10} - \frac{2xy}{10}) \\ 0 \end{pmatrix} - \begin{pmatrix} \frac{1}{10}(3y - x^2) \\ 0 \\ 0 \end{pmatrix}| = | \begin{pmatrix} 0 \\ (\frac{x}{10} - \frac{2xy}{10}) \\ 0 \end{pmatrix} | = (\frac {x}{10}) - (\frac{2xy}{10})[/math]

derecha

%Paso de muestreo
h=0.1;
%Defenicion de variables
x=(-1:h:1);
y=(0:h:1);
[mx,my]=meshgrid(x,y);
%Deformacion parabolica de la malla
yy=my.*(mx.^2+2);
mesh(mx,yy,0.*mx);
%Definimos variable U
Ux=(mx.*yy)./10;
Uy=(-yy.*mx.^2)/10;
%Tension tangencial
G=(1/10).*(-2.*yy.*mx+mx);
surf(mx,yy,G);

Los colores de la placa plana representan la ondulación de esta, la parte mas verde esta fija, la parte amarilla se eleva hacia arriba por el contrario de la zona azul que va hacia abajo, esto se puede observar mejor en la siguiente imagen representada tridimensionalmente.

11 Tensión de Von Mises

La tensión de Von Mises es un campo escalar que se emplea para analizar cómo reacciona un material específico frente a un esfuerzo, permitiendo diferenciar entre un comportamiento plástico y elástico, así como identificar el origen de un posible fallo. Se calcula a partir de los autovalores de la matriz de tensiones:

[math] σ_{VM} = \sqrt{\frac{(σ_{1}-σ_{2})^2 + (σ_{2}-σ_{3})^2 + (σ_{3}-σ_{1})^2}{2}} [/math]

Para calcularlo, primero tendremos que calcular los autovalores de dicha matriz y luego calcular los valores de Von Mises para cada punto:

%Paso de muestreo
h=0.1;
%Defenicion de variables
x=(-1:h:1);
y=(0:h:1);
[mx,my]=meshgrid(x,y);
%Deformacion parabolica de la malla
yy=my.*(mx.^2+2);
mesh(mx,yy,0.*mx);
%Definimos variable U
Ux=(mx.*yy)./10;
Uy=(-yy.*mx.^2)/10;
%Definimos la matriz de tensiones
sigma = [100 20 30; 20 50 10; 30 10 80];
%Calculamos los autovalores de la matriz de tensiones
eigenvalues = eig(sigma);
%Calculamos la tension de Von Mises
vonMisesStress = sqrt(0.5*((eigenvalues(1)) - eigenvalues (2))^2 +...
%Mostramos la tension de Von Mises 
disposición(['Tensiones de Von Mises:', 
num2str(vonMisesStress)]);
%Dibujamos un grafico en 3D para representar la matriz de tensiones y señalar el máximo 
figure;
scatter3(eigenvalues(1), eigenvalues(2), eigencalues(3), 'filled');
title('Puntos de autovalores de la matriz de tensiones');
xlabel('Autovalor 1');
ylable('Autovalor 2');
zlabel('Autovalor 3');
grid on;
%Identificamos y mostramos el punto de maxima tension de Von Mises
[maxValue, idx] = max([vonMisesStress]);
fprintf('El maximo valir de la tension de Von Mises se alcanza en el punto %d\n', idx);

centro

12 Elasticidad Lineal

Realizamos los cálculos para determinar la eslasticidad lineal, para ello realizaremos la divergencia del campo vectorial de los vectores cuyas componentes son las filas de la matriz [math] σ [/math].

[math] \vec{F} = -∇·σ [/math]

El desplazamiento es el siguiente: [math] \vec{u}(x,y) = (\frac{xy}{2} - yx^2) \vec{i} + (-\frac{y^2}{2})\vec{j}[/math]

La matriz de tensiones es la siguiente: [math]σ (x,y) = \begin{pmatrix} σ{_x}{_x} & σ{_x}{_y}\\ σ{_y}{_x} & σ{_y}{_y}\end{pmatrix}[/math]

[math]∇·σ = \begin{pmatrix} \frac{dσ{_1}{_1}}{dx} & \frac{dσ{_1}{_2}}{dy} & \frac{dσ{_1}{_3}}{dz} \\\frac{dσ{_2}{_1}}{dx} & \frac{dσ{_2}{_2}}{dy} & \frac{dσ{_2}{_3}}{dz} \\\frac{dσ{_3}{_1}}{dx} & \frac{dσ{_3}{_2}}{dy} & \frac{dσ{_3}{_3}}{dz} \end{pmatrix} = \begin{pmatrix} \frac{-2x}{10} + (\frac{-2x}{10}) + 0 \\ (\frac{1}{10}) - (\frac{2y}{10}) + (\frac{1}{10}) + 0 \\ 0 + 0 + 0 \end{pmatrix} = \begin{pmatrix} \frac{-4x}{10} \\ \frac{2}{10} - \frac{2y}{10} \\ 0 \end{pmatrix} [/math]

El resultado final teniendo en cuenta los signos:

[math] \vec{F} = -∇·σ = \begin{pmatrix} \frac{4x}{10} \\ \frac{2}{10} (y - 1) \\ 0 \end{pmatrix} [/math]

13 Masa a partir de la densidad

Para poder calcular la masa total, deberemos realizar la integral correspondiente. La densidad de la placa viene determinada por la función

[math] d(x,y) = (2-|x|)(4-y) [/math]

La integral a calcular es: [math]\displaystyle \int_{-1}^{1} \int_0^{2 + x^2} (2 - |x|)(4 - y) dydx [/math] Para poder realizar los cálculos de la interfaz hemos usado el Programa Matlab

%Paso de muestreo
h=0.1;
%Defenicion de variables
x=(-1:h:1);
y=(0:h:1);
[mx,my]=meshgrid(x,y);
%Deformacion parabolica de la malla
yy=my.*(mx.^2+2);
mesh(mx,yy,0.*mx);
%Definimos variable U
Ux=(mx.*yy)./10;
Uy=(-yy.*mx.^2)/10;
%Definimos los limites y la función f(x)
f=@(x)2 + x.^2;
%Definir la densidad d(x,y)
d=@(x,y)(2 - abs(x)).* (4-y);
%Calcular la masa total con integral 
masa_total = integral2(@(x,y)d(x,y),-1,1,0,f);
%Mostrar el resultado
fprintf('La masa total de la placa es aproximadamente: %.6f\n',masa_total);

Finalmente, tras ejecutar el código obtenemos 19,433 como resultado de la integral.

14 Bibliografía

Para realizar el trabajo nos hemos apoyado en: