Placa Plana (Grupo 26)

De MateWiki
Saltar a: navegación, buscar
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

Visualizámos la placa rectangular dada, en la que dibujaremos los distintos tipos de campos. Para poder representar esta placa hacemos uso del programa Matlab, dibujaremos el mallado determinado por la región [-2;2] x [0;3].
Para poder representar el mallado utilizamos el siguiente código:

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

El gradiente de la temperatura T ([math]\nabla T[/math]) se calcula a partir del campo escalar T, lo que genera como resultado un vector que representa la variación de la temperatura en cada dirección del espacio.

La temperatura viene dada por la siguiente expresión [math] T = (1 - x^4) (\frac {1}{2}- y)[/math], este campo escalar dependerá de las variables x e y. La representación del campo vectorial gradiente indica, que en cada punto del solido, la dirección en la cual la temperatura aumenta mas rápido, y su módulo indicara la rapidez con la que la temperatura aumenta en esa dirección.

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


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 eś 0,5000 en el punto [math](x , y) = (0.000 , 0.000)[/math]

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 dibujaremos el campo vectorial.

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 código Matlab es el siguiente:

%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);


centro

5 Campo de desplazamientos

Estudiaremos si alguno de los puntos del mallado del campo de desplazamientos se quedan fijos. Para poder comprobar esto, dibujamos un mallado que indique el movimiento de cada uno de los puntos que componen a este. El desplazamiento estará indicado por las flechas.


%Campo de desplazamiento
Ux=(mx.*yy)./10;
Uy=(-yy.*mx.^2)/10);
quiver(mx,yy,Ux,Uy,'r');
centro

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.

El código en Matlab es el siguiente:

%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
centro

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]

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

Para la obtención del rotacional, deberemos calcular tres limites, considerando tres curvas situadas en planos perpendiculares. El rotacional de un campo se puede calcular siempre y cuando este sea continuo y diferenciadle en todos sus puntos. Calcularemos el producto vectorial.


[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 haremos 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])


Ejecutando el programa en Matlab, podémoslos determinar que los puntos donde el rtoracional 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

Calculamos la tension de deformación, para ello, hayamos la derivada parcial de componente tal y como se representa en los siguientes cálculos. La formula de la Tension esta compuesta por la divergencia, esa el la razón por la cual obtendremos una matriz multiplicada por un numero.


[math] \vec{u}(x,y) = \frac{xy\vec{i}-yx^2\vec{j}}{10} = ∇•\vec{u} (\frac{xy}{10},\frac{-yx^2}{10}) [/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]

[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] ∇\vec{u} + ∇\vec{u^t}= \begin{pmatrix} \frac{y}{10} & \frac{x}{10} & 0 \\ -2x\frac{y}{10} & \frac{-x^2}{10} & 0 \\ 0 & 0 & 0 \end{pmatrix} + \begin{pmatrix} \frac{y}{10} & -2x\frac{y}{10} & 0 \\ \frac{x}{10} & -\frac{x^2}{10} & 0 \\ 0 & 0 & 0 \end{pmatrix} [/math] = 2є= [math]\begin{pmatrix} \frac{2y}{10} & (-2y+1)\frac{x}{10} & 0 \\ (-2y+1)\frac{x}{10} & -2\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 σ usando la fórmula dada en el enunciado: [math] σ=λ∇·\vec{u}1+2µϵ [/math]; teniendo tanto λ como µ valor 1.

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

center

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);


El resultado de la integral tras ejecutar el código es 19,4333