Campo de temperaturas y deformaciones en placa plana. (Grupo 11)

Trabajo realizado por estudiantes
Título	Campos de temperatura y deformaciones en placa plana(Grupo 11)
Asignatura	Teoría de Campos
Curso	2024-25
Autores	Luis Ignacio Quintana Sierra Eduardo López Reyes Rodrigo Ladrón de Guevara Solera Ángel Martín Cruz
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

Visualización de campos escalares y vectoriales en elasticidad. Consideramos una placa rectangular plana (en dimensión 2) que ocupa la región [math](x, y) ∈ [-1, 1]×[0, f(x)][/math]. Donde [math] f(x)= 2+x^2 [/math].

En ella vamos a suponer que tenemos definidas dos cantidades físicas: la temperatura [math]T(x, y)[/math], que viene dada por,

[math]T(x, y) = (1-x^4)(1/2-y),[/math] 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)= x \vec{i} + y \vec{j} [/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_{d}}(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 ondulatorio de los puntos de la misma dado por el vector:

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

1 Dibujo del mallado que representa los puntos interiores del sólido.

Resultado de ejecución

Esta gráfica muestra el mallado de la placa y el código utilizado en MatLab para obtenerlo. Un breve resumen del funcionamiento del código seria:

• La primera línea del código utiliza algo básico en Matlab, que es el uso del clear y el clc para que se borren todas las variables anteriormente usadas y no causen confusión
• En las siguientes líneas de código discretizamos las variables x e y.
• y parametrizándola entre las curvas y=0 e [math] y=2+x^2 [/math] mediante combinación lineal según la distancia a la que nos encontremos.
• Ya en las últimas líneas mostramos los resultados en forma de gráfica.

clc,clear,
h=1/10;
x = -1:h:1;
y1 = 2 + x.^2;  
y2 = zeros(size(x));  
figure;
plot(x, y1, 'b', 'LineWidth', 1); 
hold on;
plot(x, y2, 'b', 'LineWidth', 1); 
plot([1,1],[0,3],'b','LineWidth', 1);plot([-1,-1],[0,3],'b','LineWidth', 1);
ncurvas = 10; 
for k = 1:ncurvas
    alpha = k / ncurvas; 
    yy = (1 - alpha) * y1 + alpha * y2; 
    plot(x, yy); 
end
for i=1:length(x)
    plot([x(i) x(i)],[0 y1(i)])
end
axis([-2 2 0 3]);
xlabel('X');
ylabel('Y');
title('Mallado Placa Plana');
grid on;
hold off;

2 Curvas de nivel de la temperatura

Resultado de ejecución

Resultado de ejecución

La siguiente gráfica representa las curvas de nivel de la temperatura.

El gradiente de un campo vectorial es el siguiente: [math]\nabla T(x,y) =\frac{d∂}{dx} +\frac{d∂}{dy}[/math]
siendo el campo de temperatura: [math]T(x, y) = (1-x^4)(1/2-y)[/math] El gradiente será: [math]\nabla T(x,y) = -4x^3(0.5-y)\vec{i} + (1-x^4)\vec{j}[/math]

El código usado para la representación gráfica del gradiente y de la distribución de temperatura es:

clc,clear,
h=1/10;
x=-1:h:1;
f=@(x) 2+x.^2;
[X,Y]=meshgrid(x,linspace(0,1,length(x)));
for i=1:length(x)
Y(:,1)=linspace(0,f(x(i)),length(x));
end
T=(1-X.^4).*(1/2-Y);
figure;
surf(X, Y, T);
hold on
contour3(X,Y,T);
xlabel('X');
ylabel('Y');
zlabel('Temperatura(x, y)');
title('Distribución de Temperatura y Gradiente');
colorbar;
view(2);
[Tx, Ty]=gradient(T,h,h);
quiver(X,Y,Tx, Ty,'m',LineWidth=1.2);
hold off

3 Cálculo de energía calorífica con la Ley de Fourier

De acuerdo a la Ley de Fourier la energía calorífica [math] \vec{Q} [/math] viaja de acuerdo a la fórmula [math] \vec{Q}=−κ∇T [/math], donde [math] κ [/math] es la constante de conductividad térmica de la placa que supondremos [math] κ=1 [/math]. Para poder calcular la energía calorífica, primero debemos de calcular el gradiente de la temperatura, que viene dada por: [math]T(x, y) = (1-x^4)(1/2-y)[/math]

El gradiente un campo escalar se obtiene derivando parcialmente respecto de cada coordenada, por lo que el gradiente del campo es:

[math] ∇T(x,y)=(-2x^3+4yx^3)\vec{i} + (x^4-1)\vec{j}[/math] , y por lo tanto, la energía calorífica [math] \vec{Q} [/math] es: [math] \vec{Q}= (2x^3-4yx^3)\vec{i} + (-x^4+1)\vec{j}[/math]

Una vez calculado, se procede a dibujarlo utilizando Matlab:

Resultado de ejecución

clc,clear,
h = 1/10;
x = -1:h:1;  
f = @(x) 2 + x.^2;

[X, Y] = meshgrid(x, linspace(0, 1, length(x)));

for i = 1:length(x)
    Y(:, i) = linspace(0, f(x(i)), length(x));
end
T = (1 - X.^4) .* (1/2 - Y);
figure;
[Tx, Ty]=gradient(T,h,h);
Qx=-Tx; Qy=-Ty;
contour3(X,Y,T,30);  
hold on;
quiver3(X,Y,T,Qx,Qy,zeros(size(X)),'m');  
axis equal  
axis([-2, 2, 0, 3, 0, 1])
xlabel('Eje X')
ylabel('Eje Y')  
zlabel('Temperatura T(x, y)')
title('Energía Calorífica')
colorbar
grid on
hold off
view(2)

4 Punto y dirección de la máxima variación de temperatura

A continuación se adjunta el resultado gráfico, así como el código utilizado en Matlab:

Resultado de ejecución

% Parámetros iniciales
h = 1/50;
x = -1:h:1;
f = @(x) 2 + x.^2;
[X, Y] = meshgrid(x, linspace(0, 1, length(x)));

% Definición de la región válida
for i = 1:length(x)
    Y(:, i) = linspace(0, f(x(i)), length(x));
end

% Temperatura
Temperatura = (1 - X.^4) .* (0.5 - Y);

% Cálculo de derivadas parciales
ParcialX = X.^3 .* (4 .* Y - 2);  
ParcialY = X.^4 - 1;

% Magnitud del gradiente
gradiente = hypot(ParcialX, ParcialY);

% Ubicación del máximo gradiente
[MaxGradiente, MaxIdx] = max(gradiente(:));
[MaxRow, MaxCol] = ind2sub(size(gradiente), MaxIdx);
PtoMaximo = [X(MaxRow, MaxCol), Y(MaxRow, MaxCol)];

% Submuestreo para graficar
step = 5;
X2 = X(1:step:end, 1:step:end);
Y2 = Y(1:step:end, 1:step:end);
ParcialX2 = ParcialX(1:step:end, 1:step:end);
ParcialY2 = ParcialY(1:step:end, 1:step:end);

% Graficación
figure;
hold on;

% Flechas del gradiente
quiver3(X2, Y2, zeros(size(X2)), ParcialX2, ParcialY2, zeros(size(ParcialX2)), 'black', 'LineWidth', 1, 'MaxHeadSize', 1);

% Punto de máximo gradiente con flecha de dirección
plot3(PtoMaximo(1), PtoMaximo(2), 0, 'ro', 'MarkerSize', 1, 'LineWidth', 5); 
text(PtoMaximo(1), PtoMaximo(2), 0, sprintf(' Max Grad: %.2f', MaxGradiente), 'FontSize', 10, 'Color', 'k', 'VerticalAlignment', 'bottom');
quiver3(PtoMaximo(1), PtoMaximo(2), 0, ParcialX(MaxRow, MaxCol)/40, ParcialY(MaxRow, MaxCol), 0, 'r' , 'LineWidth', 2, 'MaxHeadSize', 1);

% Configuración de los ejes
axis equal;
axis([-1.5, 1.5, -0.5, 3.5, -1, 1]);
xlabel('Eje de abcisas');
ylabel('Eje de ordenadas');
zlabel('Eje de aplicadas');
title('Gradiente y Dirección del Máximo Incremento');
grid on;
hold off;

5 Representación gráfica del desplazamiento del sólido

A continuación se incluye la representación del sólido previa al movimiento y después del movimiento (en [math]t = 0[/math]), así como una figura que contiene la representación del sólido previa y tras el desplazamiento, para lo que se ha utilizado el comando subplot. El movimiento viene impuesto por el vector de desplazamientos [math]\vec{u}(x,y)=\frac{xy\vec{i}-yx^2\vec{j}}{10}[/math]. Se puede observar que los puntos con x=0 e y=0 se mantienen fijos a lo largo de todo el desplazamiento. Se adjunta el código de Matlab a continuación:

Resultado de ejecución

clc,clear,
h = 1/10;
x = -1:h:1;  
f = @(x) 2 + x.^2;
[X, Y] = meshgrid(x, linspace(0, 1, length(x)));

for i = 1:length(x)
    Y(:, i) = linspace(0, f(x(i)), length(x));
end
U_x = (X .* Y) / 10; 
U_y = (-Y .* X.^2) / 10;  

figure;
subplot(2,2,1)
quiver(X, Y, U_x, U_y);
axis equal;
xlabel('X');
ylabel('Y');
title('Campo de desplazamientos en la placa');

subplot(2,2,2)
y1 = 2 + x.^2;  
y2 = zeros(size(x));  
hold on
plot(x, y1, 'b-', 'LineWidth', 1); 
plot(x, y2, 'b-', 'LineWidth', 1); 
plot([1,1],[0,3],'b','LineWidth', 1);plot([-1,-1],[0,3],'b','LineWidth', 1);
ncurvas = 10; 
for k = 1:ncurvas
    alpha = k / ncurvas; 
    yy = (1 - alpha) * y1 + alpha * y2; 
    plot(x, yy); 
end
for i=1:length(x)
    plot([x(i) x(i)],[0 y1(i)])
end
axis([-2 2 0 3]);
xlabel('X');
ylabel('Y');
title('Mallado Placa Plana');
grid on;
hold off

subplot(2,2,3)
Xf=X+U_x; Yf=Y+U_y;
mesh(X,Y,X.*0,'EdgeColor','red');
hold on
mesh(Xf,Yf,Xf.*0);
hold off
title('Antes y después del desplazamiento');
xlabel=('X');
ylabel=('Y');
view(2)

6 Divergencia [math]∇·\vec{u}[/math]

Para calcular la divergencia usaremos la siguiente formula: [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} [/math] usada para calcular la divergencia en campos escalares.

Siendo en [math]t=0[/math]: [math]\vec u = \frac{xy}{10}\vec{i} + \frac{-yx^2}{10}\vec{j} + 0\vec{k}[/math]

[math]\ \nabla \cdot \vec u = \frac{y-x^2}{10}[/math]

En la gráfica se puede observar el cambio de volumen debido al desplazamiento, que además va variando dependiendo del punto en el que lo estudiemos.

Resultado de ejecución

clc,clear,
h = 1/10;
x = -1:h:1;  
f = @(x) 2 + x.^2;

[X, Y] = meshgrid(x, linspace(0, 1, length(x)));

for i = 1:length(x)
    Y(:, i) = linspace(0, f(x(i)), length(x));
end
Div=(Y-X.^2)/10;
surf(X,Y,Div);
title('Divergencia en t=0');
xlabel('X');
ylabel('Y');
zlabel('Divergencia')
colorbar

7 Rotacional [math]\left | ∇ \times \vec{u} \right |[/math]

Para calcular el rotacional usaremos la siguiente formula: [math]∇ × \vec{u}= \begin{vmatrix}\vec{i} & \vec{j}& \vec{k}\\ \frac{\partial }{\partial x} & \frac{\partial }{\partial y}& \frac{\partial }{\partial z}\\ u_{1}& u_{2}& u_{3} \end{vmatrix}[/math] Usada para calcular el rotacional en campos escalares.

Siendo en [math]t=0[/math]: [math]\vec u = \frac{xy}{10}\vec{i} + \frac{-yx^2}{10}\vec{j} + 0\vec{k}; [/math] [math]∇ × \vec{u}[/math] = [math]\begin{vmatrix} \vec{i} & \vec{j} & \vec{k}\\ \frac{\partial }{\partial x} & \frac{\partial }{\partial y} &\frac{\partial }{\partial z} \\ \frac{xy}{10} & \frac{-yx^2}{10} & 0\end{vmatrix} = \frac{-2yx}{10} \vec{k} - \frac{x}{10} \vec{k} = \frac{x(-2y-1)}{10}\vec{k}; [/math]

Buscamos el modulo: [math]|∇ × \vec{u}|= \frac{x(2y+1)}{10}[/math]

En la grafica se puede apreciar que los puntos de mayor rotacional son los pertenecientes a las rectas [math] y=0; y=6; y=12 [/math], representadas en amarillo.

Resultado de ejecución

clc,clear,
h = 1/10;
x = -1:h:1;  
f = @(x) 2 + x.^2;

[X, Y] = meshgrid(x, linspace(0, 1, length(x)));

for i = 1:length(x)
    Y(:, i) = linspace(0, f(x(i)), length(x));
end
Rot=(X.*(2*Y+1))/10;
surf(X,Y,Rot);
title('Rotacional');
xlabel('X');
ylabel('Y');
zlabel('Rotacional');
colorbar

8 Tensor de tensiones

Para la realización de este apartado introduciremos dos nuevos conceptos: el tensor de deformaciones, Ԑ, y el tensor de tensiones, σ, definidos a continuación,
[math]Ԑ(\vec{u}) = (∇\vec{u} + ∇\vec{u}t)/2[/math] [math]σ=λ∇·\vec{u}1 + 2μԐ[/math],

donde Ԑ será la parte simétrica del tensor [math]∇·\vec{u}[/math]; 1 es el tensor identidad en [math]R^3[/math], y ([math]λ[/math], [math]µ[/math]) son los llamados coeficientes de Lamé, que determinan las características elásticas de cada material. Para este caso tomaremos que [math]λ=µ=1[/math]. Se hallaran las tensiones normales en la dirección de los ejes cartesianos, y se graficarán en caso de ser no nulas.

Recordando el vector [math]\vec{u}=(\frac{xy}{10},\frac{-x^2y}{10}, 0)[/math] previamente definido, el primer paso será calcular su gradiente:

[math] ∇\vec{u}=\begin{pmatrix}\frac{\partial u_{1}}{\partial x} & \frac{\partial u_{1}}{\partial y} & \frac{\partial u_{1}}{\partial t} \\ \frac{\partial u_{2}}{\partial x} & \frac{\partial u_{2}}{\partial y} & \frac{\partial u_{2}}{\partial t}\\ \frac{\partial u_{3}}{\partial x} & \frac{\partial u_{3}}{\partial y} & \frac{\partial u_{3}}{\partial t}\end{pmatrix} = \begin{pmatrix}\frac{y}{10} & \frac{x}{10} & 0\\\frac{-2xy}{10} & \frac{-x^2}{10} & 0\\ 0 & 0 & 0\end{pmatrix}[/math]; [math] ∇\vec{u}t=\begin{pmatrix}\frac{y}{10} & \frac{-2xy}{10} & 0\\ \frac{x}{10} & \frac{-x^2}{10} & 0\\0 & 0 & 0\end{pmatrix} [/math]
Con estos resultados, calculamos el tensor de deformaciones:

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

Con la divergencia del campo hallada previamente, [math]∇·\vec{u}=\frac{y-x^2}{10}[/math], definimos el tensor de tensores, para luego calcular las tensiones normales a los ejes.

[math]σ = λ∇·\vec{u}1 + 2 μԐ=\begin{pmatrix}\frac{y-x^2}{10} & 0 & 0\\0 & \frac{y-x^2}{10} & 0\\ 0 & 0 & \frac{y-x^2}{10}\end{pmatrix} + \begin{pmatrix}\frac{y}{5} & \frac{x(1-2y)}{10} & 0\\\frac{x(1-2y)}{10} & \frac{-x^2}{5} & 0\\ 0 & 0 & 0\end{pmatrix} = \begin{pmatrix}\frac{3y-x^2}{10} & \frac{x(1-2y)}{10} & 0\\\frac{x(1-2y)}{10} & \frac{y-3x^2}{10} & 0\\ 0 & 0 & \frac{y-x^2}{10}\end{pmatrix}[/math]

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

Atendiendo al resultado observamos como el valor de las tensiones obtenidas es cero. Esto implica la inexistencia de dichas tensiones, y la imposibilidad de trazarlas en un gráfico.

Resultado de ejecución

Resultado de ejecución

Resultado de ejecución

Código Matlab

clc,clear,
h = 1/10;
x = -1:h:1;  
f = @(x) 2 + x.^2;

[X, Y] = meshgrid(x, linspace(0, 1, length(x)));

for i = 1:length(x)
    Y(:, i) = linspace(0, f(x(i)), length(x));
end
Teni=(3*Y-X.^2)/10;
Tenj=(Y-3*X.^2)/10;
Tenk=(Y-X.^2)/10;
figure
surf(X,Y,Teni)
axis image
colorbar
xlabel('X')
ylabel('Y') 
zlabel('Tensión dirección i')
title('Tensiones normales eje i')

figure
surf(X,Y,Tenj)
axis image
colorbar
xlabel('X')
ylabel('Y') 
zlabel('Tensión dirección j')
title('Tensiones normales eje j')

figure
surf(X,Y,Tenk)
axis image
colorbar
xlabel('X')
ylabel('Y') 
zlabel('Tensión dirección k')
title('Tensiones normales eje k')

9 Tensiones tangenciales al plano ortogonal a [math]\vec{i}[/math]

Buscaremos analizar las tensiones que sufre la placa con respecto a la dirección especificada, en t = 0. Para ello, definiremos esta tensión como [math] |σ·\vec{i}−(\vec{i}·σ·\vec{i})\vec{i}|[/math]. Recordando el resultado obtenido anteriormente, [math]\vec{i}·σ·\vec{i}= \frac{3y-x^2}{10}\vec{i}[/math], concluimos que la tensión buscada es igual a [math]σ·\vec{i}[/math]. Continuando el desarrollo,

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

Haciendo la operación queda por tanto: [math] |σ·\vec{i}−(\vec{i}·σ·\vec{i})\vec{i}| = |\begin{pmatrix} \frac{3y-x^2}{10} \\ \frac{x(1-2y)}{10} \\ 0 \end{pmatrix} - \begin{pmatrix} \frac{3y-x^2}{10} \\ 0 \\ 0 \end{pmatrix}| = |\frac{x(1-2y)}{10}|[/math]

Resultado de ejecución

Resultado de ejecución

clc,clear,
h = 1/10;
x = -1:h:1;  
f = @(x) 2 + x.^2;

[X, Y] = meshgrid(x, linspace(0, 1, length(x)));

for i = 1:length(x)
    Y(:, i) = linspace(0, f(x(i)), length(x));
end
Teni=(X-2.*X.*Y)/10;
figure
surf(X,Y,Teni)
axis image
colorbar
xlabel('X')
ylabel('Y') 
zlabel('Tensión')
title('Tensiones tangenciales al plano ortogonal a i')

10 Tensión de Von Mises

La Tensión de Von Mises es una magnitud física escalar usada en campos como la ingeniería estructural, calculada a partir de los valores de las tensiones principales de cada punto del espacio, la cual indica la tensión a aplicar a cada punto de un material para que éste inicie su comportamiento plástico.

Sea la Tensión de Von Mises [math] \sigma_{VM} [/math], calculada para un punto P de un sólido deformable, y sean las tensiones principales del tensor tensión para dicho punto [math]\sigma_1, \sigma_2, \sigma_3[/math], correspondientes a los autovalores de [math] \sigma_{ij} [/math], entonces se comprueba que la Tensión de Von Mises viene dada por la siguiente expresión:

[math] \sigma_{VM} = \sqrt{ \frac { ( \sigma_1 - \sigma_2 )^2 + ( \sigma_2 - \sigma_3 )^2 + ( \sigma_3 - \sigma_1 )^2 }{2} }[/math]

Repetimos el código como en anteriores apartados, añadiendo esta vez la función de Von Mises ([math]VonMises[/math]).

Resultado de ejecución

clc,clear,
h = 1/10;
X = -1:h:1;  
f = @(x) 2 + x.^2;

[x, y] = meshgrid(X, linspace(0, 1, length(X)));

% Generar la matriz y en función de X
for i = 1:length(X)
    y(:, i) = linspace(0, f(X(i)), length(X));
end

Ux=(x.*y)/10;
Uy=(-y.*x.^2)/10;

[gradUxx, gradUxy]=gradient(Ux,h,h);
[gradUyx, gradUyy]=gradient(Uy,h,h);

% Tensor de deformaciones ε (simétrica)
exx=(gradUxx+gradUxx)/2;
eyy=(gradUyy+gradUyy)/2;
exy=(gradUxy+gradUyx)/2;

% La fórmula del tensor de tensiones σ = ∇·u + 2ε
divU=gradUxx+gradUyy;  % Divergente de u
sigmaxx=divU+2*exx;
sigmayy=divU+2*eyy;
sigmaxy=divU+2*exy;

% Tensión de Von Mises para cada punto en la malla
sigmavm=sqrt(((sigmaxx-sigmayy).^2+(sigmayy-sigmaxy).^2+(sigmaxy-sigmaxx).^2)/2);

% Graficar la tensión de Von Mises
figure;
surf(x, y, sigmavm);
xlabel('X');
ylabel('Y');
zlabel('Tensión de Von Mises');
title('Distribución de la Tensión de Von Mises');
colorbar;

% Encontrar el punto con el valor máximo de la tensión de Von Mises
[maxval,idx]=max(sigmavm(:));
[maxx,maxy]=ind2sub(size(sigmavm),idx);

% Señalar el punto máximo en el gráfico
hold on;
plot3(x(maxx,maxy),y(maxx,maxy),maxval, 'ro', 'MarkerFaceColor', 'r');

% Mostrar el valor máximo de la tensión de Von Mises
disp(['El valor máximo de la Tensión de Von Mises se alcanza en X = ',num2str(x(maxx,maxy)), ' y Y = ',num2str(y(maxx,maxy))]);
disp(['Valor máximo de Von Mises: ',num2str(maxval)]);

11 Campo de fuerzas que actúa sobre la placa

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}=-∇· σ[/math]

donde [math]∇ · σ[/math] es el campo vectorial que se obtiene al hacer la divergencia de los vectores cuyas componentes son las filas de la matriz [math]σ[/math].

Calculando [math] - ∇ · σ[/math] queda [math] - ∇ · σ = \begin{pmatrix}\frac{2x}{5} \\ \frac{y}{5} - \frac{1}{5} \\ 0\end{pmatrix} [/math]

Esta matriz de fuerzas nos dice que las fuerzas en la dirección de [math]\vec{i}[/math] la fuerza es proporcional a x, en la dirección de [math]\vec{j}[/math] la fuerza depende linealmente de y, y en la dirección de [math]\vec{k}[/math] no hay fuerza.

Resultado de ejecución

El código de Matlab es:

clc,clear,
h = 1/10;
x = -1:h:1;  
f = @(x) 2 + x.^2;
[X, Y] = meshgrid(x, linspace(0, 1, length(x)));
for i = 1:length(x)
    Y(:, i) = linspace(0, f(x(i)), length(x));
end
Fx = 2*X/5;
Fy = (Y-1)/5;
Fz = zeros(size(X));

figure;
quiver(X, Y, Fx, Fy, 'AutoScale', 'on', 'AutoScaleFactor', 0.5);
xlabel('X');
ylabel('Y');
title('Fuerzas sobre la placa plana');
axis equal;
grid on;

12 Masa total de la placa

Para este último apartado se estudiará la masa total de la placa a partir de su función de densidad [math]d(x,y)=(2-|x|)(4-y)[/math]. Sólo tendremos que integrar la densidad entre la región en la que se encuentra nuestra placa [math][-1,1] [/math] x [math][0,f(x)] [/math] donde [math] f(x)=2+x^2[/math]:

[math]\int_{-1}^{1} \int_{0}^{2+x^2} (2-|x|)(4-y) \, \mathrm{d}y \, \mathrm{d}x[/math]

Primero resolvemos la integral interna que depende de y:

[math]\int_{0}^{2+x^2} (4-y) \, \mathrm{d}y \ = \left[4y - \frac{y^2}{2}\right]_{0}^{2+x^2} = (4(2+x^2)-\frac{(2+x^2)^2}{2}) = 6+2x^2-\frac{x^4}{2}[/math]

Ahora podemos calcular la integral externa que depende de x que es:

[math]\int_{-1}^{1} (2-|x|)(6+2x^2-\frac{x^4}{2}) \, \mathrm{d}x[/math]

Al ser el integrando simétrico al eje x=0, reescribimos la integral:

[math]2\int_{0}^{1} (x^4-4x^2-12)(x-2) \, \mathrm{d}x = 2\int_{0}^{1} \left[12-6x+4x^2-2x^3-x^4-\frac{x^5}{2}\right] \, \mathrm{d}x[/math]

Integrando esto nos queda:

[math]2 \left[12x-3x^2+\frac{4}{3}x^3-\frac{1}{2}x^4-\frac{x^5}{5}-\frac{x^6}{12}\right]_{0}^{1} = 2\left[12-3+\frac{4}{3}-\frac{1}{2}-\frac{1}{5}-\frac{1}{12}\right] [/math]

La masa total será [math] M=2*(9.55)=19.1[/math] unidades de masa aproximadamente.

[math]\boxed {M = 19.1 u}[/math]

Si la calculamos con el ordenador nos sale 19.4333 u de masa.