Diferencia entre revisiones de «Placa Plana (Grupo 26)»

Revisión actual del 00:26 10 dic 2024

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:

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

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]

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.

• 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 va en la dirección del vector [math]\vec{k} [/math], esto se debe a que el rotacional será perpendicular tanto a la componente [math]\vec{i} [/math] y [math]\vec{j} [/math] respectivamente porque [math]\vec{u} [/math] solo tiene componentes [math]\vec{i} [/math] y [math]\vec{j} [/math].

A continuación, representamos 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]

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.

center

%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

Hallamos las tensiones tangenciales respecto al plano ortogonal a [math] \vec{i} [/math]. Los cálculos son los siguientes: [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]

Los colores de la placa plana representan la ondulación de esta, la zona turquesa se mantiene prácticamente fija, la zona amarilla se eleva positivamente y la zona azul se curva hacia abajo.

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

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:

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

derecha

% Calculo de integral
% Paso de muestreo
h = 0.1;
% Definición de variables
x = -1:h:1;
y = 0:h:1;
[mx, my] = meshgrid(x, y);
% Deformación parabólica de la malla
yy = my .* (mx.^2 + 2);
% Inicializar matriz para la tensión de Von Mises
sigma_v = zeros(size(mx));
% Bucle para calcular la tensión en cada punto
for i = 1:size(mx, 1)
    for j = 1:size(mx, 2)
        % Coordenadas actuales
        xi = mx(i, j);
        yi = my(i, j);
        % Tensor de tensiones (ejemplo simbólico, ajusta según el problema)
        sigma = [2*xi, xi*yi; xi*yi, 3*yi]; % Define tu tensor aqui
        % Calcular autovalores del tensor de tensiones
        lambda = eig(sigma);
        % Calcular la tensión de Von Mises
        sigma_v(i, j) = sqrt(lambda(1)^2 - lambda(1)*lambda(2) + lambda(2)^2);
    end
end
% Encontrar el valor máximo de la tensión de Von Mises y su posición
[max_sigma_v, idx] = max(sigma_v(:));
[max_row, max_col] = ind2sub(size(sigma_v), idx);
max_x = mx(max_row, max_col);
max_y = my(max_row, max_col);
% Mostrar resultados
fprintf('El valor máximo de la tensión de Von Mises es %.4f\n', max_sigma_v);
fprintf('Se alcanza en el punto (%.4f, %.4f)\n', max_x, max_y);
% Graficar la tensión de Von Mises
figure;
surf(mx, my, sigma_v, 'EdgeColor', 'none');
colorbar;
title('Tensión de Von Mises');
xlabel('x');
ylabel('y');
zlabel('\sigma_v');
view(3);

El valor máximo de la tensión de Von Mises es 4.6904. se alcanza en el punto (-1.0000, 1.0000).

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

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]

derecha

% Calculo de integral
% Paso de muestreo
h = 0.1;
% Definición de variables
x = (-1:h:1); % Rango de x
y = (0:h:1);  % Rango de y
[mx, my] = meshgrid(x, y); % Crear mallas de las coordenadas
% Deformación parabólica de la malla
yy = my .* (mx.^2 + 2);
mesh(mx,yy,0.*mx);
% Inicialización de la matriz sigma como una estructura tridimensional
% Matriz 3D: [filas, columnas, filas_sigma, columnas_sigma]
sigma = zeros([size(mx), 3, 3]); 
% Definición directa de los componentes de sigma
sigma(:, :, 1, 1) = (3 .* my - mx.^2) .* (1/10);   % Elemento (1,1)
sigma(:, :, 1, 2) = (-2 .* my + 1) .* (1/10);      % Elemento (1,2)
sigma(:, :, 1, 3) = 0;                            % Elemento (1,3)
sigma(:, :, 2, 1) = sigma(:, :, 1, 2);            % Elemento (2,1) 
sigma(:, :, 2, 2) = (my - 3 .* mx.^2) .* (1/10);  % Elemento (2,2)
sigma(:, :, 2, 3) = 0;                            % Elemento (2,3)
sigma(:, :, 3, 1) = 0;                            % Elemento (3,1)
sigma(:, :, 3, 2) = 0;                            % Elemento (3,2)
sigma(:, :, 3, 3) = (my - mx.^2) .* (1/10);       % Elemento (3,3)
% Visualización de la matriz sigma en el centro de la malla
disp('Matriz sigma en el centro de la malla:');
disp(sigma(round(end/2), round(end/2), :, :));
% Calcular derivadas parciales
[dsigma11_dx, dsigma11_dy] = gradient(sigma(:, :, 1, 1), h, h); % Derivadas de sigma11
[dsigma12_d]x, dsigma12_dy] = gradient(sigma(:, :, 1, 2), h, h); % Derivadas de sigma12
[dsigma21_dx, dsigma21_dy] = gradient(sigma(:, :, 2, 1), h, h); % Derivadas de sigma21
[dsigma22_dx, dsigma22_dy] = gradient(sigma(:, :, 2, 2), h, h); % Derivadas de sigma22
% Divergencia de cada componente (vector de divergencia)
div_sigma_x = dsigma11_dx + dsigma12_dy; % Componente x de div(sigma)
div_sigma_y = dsigma21_dx + dsigma22_dy; % Componente y de div(sigma)
% Multiplicar la divergencia por -1
Fx = -div_sigma_x;
Fy = -div_sigma_y;
% Graficar el resultado como un campo vectorial
hold on 
quiver(mx, my, Fx, Fy, 'AutoScale', 'on');
xlabel('x');
ylabel('y');
title('Campo F');
grid on;
% Resultados
F = (cat(3, Fx, Fy));
% Visualización de la divergencia en un punto
disp('Divergencia de sigma en el centro de la malla:');
disp(squeeze(F(round(end/2), round(end/2), :)));
view(2);

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

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

Tras ejecutar el código, obtenemos que la masa total es 19,433.

14 Bibliografía

Para realizar el trabajo nos hemos apoyado en las siguientes aplicaciones, páginas web y documentos:

• Matlab - MathWorks
• Apuntes propios de clase
• Trabajos de años anteriores (https://mat.caminos.upm.es/wiki/Temperatura_en_una_placa_semicircular_en_2-D_(Grupo_19A) )
• Chat GPT