Presa Triangular (Grupo 27)

De MateWiki
Saltar a: navegación, buscar
Trabajo realizado por estudiantes
Título Presa Triangular . Grupo 27
Asignatura Teoría de Campos
Curso 2024-25
Autores Rafael Moreno Orellana
Miguel Rubio Arraztio
Victor Jesus Sepulveda Fernandez
Pablo Cortina Gomez
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura
  1. Consideramos la sección transversal de una presa triangular (ver figura 2).
  2. Dibujar la estructura, que es el triángulo de la figura, manteniendo las líneas coordenadas tal y como están en la figura.
  3. Suponer lo siguiente:

• Parametrizar la superficie [math](x, y) ∈ [0, 2] × [0, f(x)][/math] con:[math]f(x) = \min(3, \frac{3}{2}(2 − x))[/math].

• La temperatura viene dada por la función: [math]T(x, y) = \frac{y \cdot x^2}{2}[/math].

• Los desplazamientos se corresponden con el campo: [math]\vec{u}(x, y) = \frac{2(2 − x)y \cdot \vec{i} − y \cdot \vec{j}}{50}[/math].

• Tomar como densidad: [math]d(x, y) = (2 − |x − \frac{1}{2}|)(4 − y)[/math].

1 Representacion del mallado en 2D

Dibujar un mallado que represente los puntos interiores del sólido. Para ello, parametrizar el sólido de manera que las líneas coordenadas sean las mismas que las dibujadas en 1. Tomar los ejes (comando axis) en el rectángulo [math](x, y) ∈ [−2, 2] × [0, 3][/math] y como paso de muestreo [math]h = \frac{1}{10}[/math] para las variables [math]x[/math] e [math]y[/math].

mallado de la presa triangular
clc; clear all;
% Definimos la función f(x)
f = @(x) min(3, 3/2 * (2 - x));

% Paso de discretización y número de puntos
h = 1/10;  
Nx = ((2 - 0) / h);  
Ny = ((3 - 0) / h);  

% Discretizamos los rangos de x e y
x = linspace(0, 2, Nx);  
y = linspace(0, 3, Ny);  

% Malla de puntos (X, Y) 
[X, Y] = meshgrid(x, y);   

% Definimos Z solo para los valores dentro de f(x)
for i = 1:Nx
    Z(Y <= f(x(i)) & X == x(i)) = 0; 
end

% Configuración de la gráfica
figure('Name', 'Mallado de la Figura');  
mesh(X, Y, Z);  
xlabel('Eje X'); ylabel('Eje Y'); 
title('Superficie parametrizada: [0, 2] × [0, f(x)]');
grid on;
view(2);  
axis([-1 3 -1 4]);  
hold on;
plot([0, 2], [0, 0], 'k-', [0, 2], [f(0), f(2)], 'k-', [0, 0], [0, f(0)], 'k-', [2, 2], [0, f(2)], 'k-', 'LineWidth', 2);
hold off;


2 Curvas de Nivel de la Temperatura y [math]\nabla T[/math]

Dibujar las curvas de nivel de la temperatura (comando contour) y decidir en qué punto la temperatura es máxima a partir de la gráfica. Calcular [math]\nabla T[/math] y pintarlo como campo vectorial en la misma gráfica. Observar gráficamente que [math]\nabla T[/math] es ortogonal a dichas curvas. (Nota: es muy importante que la escala de los ejes sea la misma para apreciar el ángulo que forman las curvas de nivel y el gradiente).

2.1 Curvas de Nivel de la Temperatura

Curvas de nivel 2D
% Paso de discretización y número de puntos
h = 1/10; 
x = 0:h:2; 
y = 0:h:3; 

% Malla de puntos (X, Y)
[X, Y] = meshgrid(x, y);

% Función límite de la región válida
f = @(x) min(3, (3/2)*(2 - x));   

% Temperatura en cada punto
T = (Y .* X.^2) / 2;  
T(~region) = NaN;  

% Curvas de nivel
contour3(X, Y, T, 20);  
colorbar; 

% Configuración de la gráfica
axis equal; 
xlabel('Eje X');  
ylabel('Eje Y'); 
zlabel('Temperatura T(x, y)');  
title('Curvas de Nivel de la Temperatura');  
grid on;  
hold off;


2.2 [math]\nabla T[/math] (Gradiente)

Gradiente
Curvas de nivel y Gradiente
% Paso de discretización y número de puntos
h = 1/10;  
x = 0:h:2; 
y = 0:h:3; 

% Malla de puntos (X, Y)
[X, Y] = meshgrid(x, y); 

% Función límite de la región válida
f = @(x) min(3, (3/2)*(2 - x));   

% Temperatura en cada punto
T = (Y .* X.^2) / 2;  
T(~region) = NaN;  

% Gradiente de la temperatura
[Tx, Ty] = gradient(T, h);  

% Curvas de nivel
contour3(X, Y, T, 20); 

% Configuración de la gráfica
hold on;
quiver3(X, Y, T, Tx, Ty, zeros(size(Tx)), 'r', 'LineWidth', 1.5);  
axis equal; 
axis([0, 2, 0, 3, 0, 1]); 
xlabel('Eje X'); 
ylabel('Eje Y'); 
zlabel('Temperatura T(x, y)'); 
title('Curvas de Nivel y Gradiente de Temperatura'); 
colorbar; 
grid on; 
view(2);
hold off;


3 Ley de Fourier

De acuerdo a la Ley de Fourier la energía calorífica [math]\vec{Q} = -\kappa \nabla T[/math] donde [math]\kappa[/math] es la constante de conductividad térmica de la placa que supondremos [math]\kappa = 1[/math]. Calcular [math]\vec{Q}[/math] y dibujarlo como campo vectorial.

Q y Curvas de nivel
Q
% Paso de discretización y número de puntos
h = 1/10;  
x = 0:h:2; 
y = 0:h:3; 

% Malla de puntos (X, Y)
[X, Y] = meshgrid(x, y);  

% Función límite de la región válida
f = @(x) min(3, (3/2)*(2 - x));   

% Temperatura en cada punto
T = (Y .* X.^2) / 2;  
T(~region) = NaN; 

% Gradiente de la temperatura
[Tx, Ty] = gradient(T, h);  

% Flujo de energía calorífica
Qx = -Tx;  
Qy = -Ty;  

% Curvas de nivel
contour3(X, Y, T, 20); 

% Configuración de la gráfica
hold on;
quiver3(X, Y,T,Qx,Qy, zeros(size(Qx)), 'r'); 
axis equal;  
axis([0, 2, 0, 3, 0, 1]);  
xlabel('Eje X'); 
ylabel('Eje Y');  
zlabel('Temperatura T(x, y)'); 
title('Curvas de Nivel y Flujo de Energía');  
colorbar;  
grid on; 
view(2);
hold off;


4 Punto de Máximo Valor de [math]\nabla T[/math] (Gradiente)

Determinar analíticamente (o con el ordenador) en qué punto y en qué dirección la variación de temperatura es mayor y dibujarlo en la gráfica de la placa con un punto rojo.

Máximo gradiente 2D
Máximo gradiente 3D
% Paso de discretización y rango de puntos
h = 1/50;  
x = 0:h:2; 
y = 0:h:3; 

% Malla de puntos (X, Y)
[X, Y] = meshgrid(x, y);

% Función límite para definir la región válida
f = @(x) min(3, 3/2 * (2 - x));
Region = (Y <= f(X));  
X(~Region) = NaN; 
Y(~Region) = NaN;    

% Cálculo de la temperatura y gradiente
T = (Y .* X.^2) / 2; 
gradX = Y .* X; 
gradY = X.^2 / 2; 
GRAD = sqrt(gradX.^2 + gradY.^2);  

% Identificar el punto con el mayor valor de gradiente
[MAXgrad, MAXI] = max(GRAD(:)); 
[MAXgradX, MAXgradY] = deal(X(MAXI), Y(MAXI));  

% Submuestreo para las flechas
XGrid = X(1:5:end, 1:5:end); 
YGrid = Y(1:5:end, 1:5:end); 
GradientX = gradX(1:5:end, 1:5:end); 
GradientY = gradY(1:5:end, 1:5:end); 

% Representación Campo de Temperaturas y Punto con Mayor Gradiente
figure; 
hold on;
quiver3(XGrid, YGrid, zeros(size(XGrid)), GradientX, GradientY, zeros(size(GradientX)), 'b', 'LineWidth', 1);

plot3(MAXgradX, MAXgradY, 0, 'ro', 'MarkerSize', 10, 'LineWidth', 3);  
text(MAXgradX, MAXgradY, 0, sprintf('Max Grad = %.2f', MAXgrad), 'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom');

axis equal;
axis([-1, 3, -1, 4, -1, 1]); 
grid on;
xlabel('Eje X'); 
ylabel('Eje Y'); 
zlabel('Eje Z');
title('Campo de Temperatura y Dirección de Máxima Variación');
view(3); 
hold off;


5 Campo de Desplazamientos

Dibujar el campo de desplazamientos en los puntos del mallado del sólido. ¿Hay algunos puntos que se encuentren fijos?

Desplazamientos 3D
Desplazamientos 2D
% Paso de discretización y número de puntos 
h = 1/10;  
x = 0:h:2; 
y = 0:h:3;

% Malla de puntos (X, Y)
[X, Y] = meshgrid(x, y);

% Función límite de la región válida
f = @(x) min(3, 3/2 * (2 - x));
Region = (Y <= f(X));  
X(~Region) = NaN; 
Y(~Region) = NaN;    
Z = zeros(size(X)); 

% Campo de Desplazamientos
ux = (2 * (2 - X) .* Y) / 50; 
uy = -Y / 50; 

% Superficie representada como malla
mesh(X, Y, Z, 'EdgeColor', 'b'); 
hold on;
quiver3(X, Y, Z, ux, uy, zeros(size(ux)), 1.5, 'g');  
axis equal;  
axis([-0.5, 2.5, -0.5, 3.5, -0.5, 0.5]); 
xlabel('Eje X');  
ylabel('Eje Y');  
zlabel('Eje Z');  
title('Campo de Desplazamientos'); 
view(3);  
grid on;  
hold off;


6 Sólido Antes y Después del Desplazamiento

Dibujar el sólido antes y después del desplazamiento dado por el campo de vectores [math]\vec{u}[/math]. Dibujar ambos en la misma figura usando el comando subplot para poder compararlos.

Desplazamientos de la placa
% Paso de discretización y número de puntos 
h = 1/10;  
x = 0:h:2; 
y = 0:h:3;

% Malla de puntos (X, Y)
[X, Y] = meshgrid(x, y);

% Función límite de la región válida
f = @(x) min(3, 3/2 * (2 - x));
Region = (Y <= f(X));  
X(~Region) = NaN; 
Y(~Region) = NaN;    

% Campo de Desplazamientos
ux = (2 * (2 - X) .* Y) / 50; 
uy = -Y / 50; 

% Sólido antes del desplazamiento
subplot(1, 2, 1);
mesh(X, Y, zeros(size(X)));
view(2);
title("Sólido antes del desplazamiento");
axis([-1, 3, -1, 4]);
xlabel('Eje X'); 
ylabel('Eje Y');
grid on;

% Sólido después del desplazamiento
subplot(1, 2, 2);
mesh(X + ux, Y + uy, zeros(size(X)));
view(2);
title("Sólido después del desplazamiento"); 
axis([-1, 3, -1, 4]);
xlabel('Eje X'); 
ylabel('Eje Y');
grid on;


7 Divergencia

Dibujar [math]\nabla \cdot \vec{u}[/math] en t = 0. Determinar analíticamente los puntos en los que la divergencia de [math]\vec{u}[/math] es máxima, mínima y nula. La divergencia es una medida del cambio de volumen local debido al desplazamiento. ¿Se puede apreciar esto en la gráfica?


Para calcular la divergencia de un campo vectorial lo primero es definir los conceptos: La divergencia de un campo vectorial queda definida por:

[math]\vec{u} = \frac{\partial u_x}{\partial x} + \frac{\partial u_y}{\partial y} [/math]


Teniendo en cuenta que nuestro campo vectorial es:

[math] \vec u=\frac {1}{50}[(2(2-x)y)(\vec i) -y(\vec j)][/math]

Podemos definir cada uno de los términos de la divergencia como:

[math]u_x(x, y) = \frac {1}{50}[(2(2-x)y)][/math]

[math]u_y(x, y) = \frac {-y}{50}[/math]

Una vez que tenemos los términos definidos tenemos que calcular las derivadas parciales para obtener la divergengia.

[math]\frac{\partial u_x}{\partial x} = \frac{\partial (xy)}{\partial x} = \frac {-2y}{50}[/math]

[math]\frac{\partial u_y}{\partial y} = \frac{\partial (-yx^2)}{\partial y} = \frac {-1}{50}[/math]

Tras calcular las derivadas parciales podemos definir la divergencia: [math]\nabla \cdot \vec{u}(x, y) = \frac {-2y-1}{50}[/math]

Una vez que tenemos la divergencia calculada utilizamos Matlab para obtener su representación gráfica:

Divergencia 2D y 3D
% Paso de discretización y número de puntos 
h = 1/10;  
x = 0:h:2; 
y = 0:h:3;

% Malla de puntos (X, Y)
[X, Y] = meshgrid(x, y);

% Función límite de la región válida
f = @(x) min(3, 3/2 * (2 - x));
Region = (Y <= f(X));  
X(~Region) = NaN; 
Y(~Region) = NaN;   

% Divergencia
div_u = (-2 .* Y - 1) ./ 50;

% Configuración de la gráfica
figure;
surf(X, Y, div_u); 
colorbar;
title('Divergencia del campo u(x, y)');
xlabel('x');
ylabel('y');
zlabel('z');
view(3);


8 [math]|\nabla \times \vec{u}|[/math]

Calcular [math]|\nabla \times \vec{u}|[/math] en todos los puntos del sólido y dibujarlo. ¿Qué puntos sufren un mayor rotacional?

Sabiendo que [math] \vec u=\frac {1}{50}[(2(2-x)y)(\vec i) -y(\vec j)][/math]

Y que; por tanto: [math]u_x(x, y) = \frac {1}{50}[(2(2-x)y)][/math]

[math]u_y(x, y) = \frac {-y}{50}[/math]


[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_1 & u_2 & u_3 \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 {1}{50}[(2(2-x)y)] & \frac {-y}{50} & 0 \end{matrix}\right| =\left( \frac{\partial 0}{\partial y} - \frac{ \partial \frac {-y}{50}}{\partial z}\right)\vec i + \left(\frac{\partial \frac {1}{50}[(2(2-x)y)]}{\partial z} - \frac{\partial 0}{\partial x}\right)\vec j + \left(\frac{ \partial \frac {-y}{50}}{\partial x} - \frac{ \partial \frac {1}{50}[(2(2-x)y)]}{\partial y}\right)\vec k =0\vec i + 0\vec j + \left(\frac {2-x}{25}\right)\vec k [/math]


[math]|\nabla \times \vec{u}|[/math] = [math] \frac {2-x}{25} [/math]

Rotacional 3D Y 2D
% Paso de discretización y número de puntos
h = 1/10;  
x = 0:h:2; 
y = 0:h:3; 

% Malla de puntos (X, Y)
[X, Y] = meshgrid(x, y);  

% Función límite de la región válida
f = @(x) min(3, (3/2)*(2 - x));  
Region = (Y <= f(X));  
X(~Region) = NaN; 
Y(~Region) = NaN;    

%Módulo del rotacional
ROT = (2 - X) / 25;

% Configuración de la gráfica
figure
surf(X,Y,ROT)
view(3)
title("Rotacional 3-D en t=0")
colorbar
xlabel('Eje X')
ylabel('Eje Y')
zlabel('Eje Z')


9 Apartado 9

Definamos [math]\epsilon(\mathbf{u}) = \frac{\nabla \mathbf{u} + (\nabla \mathbf{u})^\top}{2}[/math], la parte simétrica del tensor gradiente de [math]\mathbf{u}[/math], conocido como el tensor de deformaciones. En un medio elástico lineal, isótropo y homogéneo, los desplazamientos permiten escribir el tensor de tensiones [math]\sigma_{ij}[/math] a través de la fórmula:

[math]\sigma = \lambda (\nabla \cdot \mathbf{u}) \mathbf{1} + 2\mu \epsilon,[/math]

donde [math]\mathbf{1}[/math] es el tensor identidad en el conjunto de vectores libres del espacio [math]\mathbb{R}^3[/math] y [math]\lambda[/math], [math]\mu[/math] son los conocidos como coeficientes de Lamé que dependen de las propiedades elásticas de cada material.

A pesar de que los desplazamientos son planos (es decir, [math]\mathbf{u}[/math] no tiene componente en la dirección de [math]\mathbf{k}[/math]), las tensiones no tienen por qué ser planas y puede haber tensiones en la dirección ortogonal al plano de la placa. Tomando [math]\lambda = 1[/math] y [math]\mu = 1[/math], dibujar: Las tensiones normales en la dirección que marca el eje [math]\mathbf{i}[/math], es decir, [math]\mathbf{i} \cdot \sigma \cdot \mathbf{i}[/math]. Las tensiones normales en la dirección que marca el eje [math]\mathbf{j}[/math], es decir, [math]\mathbf{j} \cdot \sigma \cdot \mathbf{j}[/math]. Las tensiones correspondientes al eje [math]\mathbf{k}[/math], es decir, [math]\mathbf{k} \cdot \sigma \cdot \mathbf{k}[/math] (dibujar las que no son nulas).

10 Tensiones Perpendiculares

Calcular las tensiones tangenciales respecto al plano ortogonal a [math]\vec{i}[/math], es decir,[math]|\sigma \cdot \vec{i} − (\vec{i} \cdot \sigma \cdot \vec{i}) \vec{i}|[/math]. Dibujar sólo las que no son nulas.

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 los autovalores de [math]\sigma[/math] (también conocidos como tensiones principales). Se trata de una magnitud escalar que se suele usar como indicador para saber cuándo un material sufre una mayor tensión y puede iniciar un comportamiento plástico (y no elástico puro). Pintar la tensión de Von Mises y señalar en qué punto se alcanza el mayor valor. (Para calcular autovalores con OCTAVE/MatLab usar el comando eig.m).

Tensión Von Mises y punto máximo
% Paso de discretización y rango de puntos
h = 1/10;  
x = 0:h:2; 
y = 0:h:3; 

% Malla de puntos (X, Y)
[X, Y] = meshgrid(x, y);

% Función límite para definir la región válida
f = @(x) min(3, 3/2 * (2 - x));
Region = (Y <= f(X));  
X(~Region) = NaN; 
Y(~Region) = NaN; 

%Autovalores
AV_1 = zeros(length(y), length(x)); 
AV_2 = zeros(length(y), length(x));
AV_3 = zeros(length(y), length(x));

for i = 1:length(x)
    for j = 1:length(y)
        T_ij = [(-3*Y(j,i))/25 - 1/25, (4*(2 - X(j,i)))/25, 0; 
                (4*(2 - X(j,i)))/25, (-Y(j,i) - 2)/25, 0;
                0, 0, 0];
        AVs = eig(T_ij);
        AV_1(j, i) = AVs(1);
        AV_2(j, i) = AVs(2);
        AV_3(j, i) = AVs(3);
    end
end

%Tensión de Von Mises
VM = sqrt(((AV_1 - AV_2).^2 + (AV_2 - AV_3).^2 + (AV_3 - AV_1).^2) / 2);

% Configuración de la gráfica
figure;
surf(X, Y, VM);
xlabel('x');
ylabel('y');
zlabel('Von Mises');
title('Tensión de Von Mises');

% Identificar el punto con el mayor valor
[max_tension, idx] = max(VM(:)); 
[row, col] = ind2sub(size(VM), idx);
max_x = X(row, col);
max_y = Y(row, col);
hold on;

%Dibujar el punto de mayor valor
scatter3(max_x, max_y, max_tension, 100, 'r', 'filled');
disp(['Valor máximo de la tensión: ', num2str(max_tension)]);
disp(['Coordenadas del valor máximo: (x, y) = (', num2str(max_x), ', ', num2str(max_y), ')']);


El punto de mayor valor, como podemos ver en la gráfica se encuentra en (x,y) = (0,3).

12 Campo de fuerzas [math]\vec{F}[/math]

El campo de fuerzas [math]\mathbf{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]\mathbf{F} = -\nabla \cdot \sigma[/math], donde [math]\nabla \cdot \sigma[/math] es el campo vectorial que se obtiene al hacer la divergencia de los vectores cuyas componentes son las filas de la matriz [math]\sigma[/math]. Calcular [math]\mathbf{F}[/math] y dibujarlo. 

% Paso de discretización y rango de puntos
h = 1/10;  
x = 0:h:2; 
y = 0:h:3; 

% Malla de puntos (X, Y)
[X, Y] = meshgrid(x, y);

% Función límite para definir la región válida
f = @(x) min(3, 3/2 * (2 - x));
Region = (Y <= f(X));  
X(~Region) = NaN; 
Y(~Region) = NaN;  

%Definición de F
F1 = (-3*Y)/25 - 1/25 + (4*(2 - X))/25;
F2 = (4*(2 - X))/25 + (-Y - 2)/25;
[div_x1, div_y_1] = gradient(F1, x, y);
[div_2x, div_2y] = gradient(F2, x, y);
divergencia = div_x1 + div_2y;
F_x = -divergencia;
F_y = -divergencia;
Fuerza = sqrt(F_x.^2 + F_y.^2);
max_F = max(Fuerza(:));

% Configuración de la gráfica
figure;
quiver(X, Y, F_x, F_y);
xlabel('x');
ylabel('y');
title('Fuerza en la región triangular');


13 Cálculo de Masa por Integración

Supongamos que la densidad de la placa es [math]d(x, y) = (2 - |x|)(4 - y)[/math]. Calcular la masa total aproximando la integral correspondiente numéricamente:

[math]M = \iint_{\text{Región}} d(x, y) \, dx \, dy[/math]. 

%Establecemos valores límite
x_min = 0;
x_max = 2;
y_min = 0;
y_max = 3;

%Cálculo de masa
d = @(x, y) (2 - abs(x - 1/2)) .* (4 - y);
masatotal = integral2(d, x_min, x_max, y_min, y_max);
disp(['La masa total es: ', num2str(masatotal)]);


El valor de la masa obtenido por integración 20,625 unidades de masa.

Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Presa_Triangular_(Grupo_27)&oldid=79297»