Representación y Análisis de Campos Físicos en una Columna Recta. Grupo 13

Trabajo realizado por estudiantes
Título	Representación y Análisis de Campos Físicos en una Columna Recta. Grupo 13
Asignatura	Teoría de Campos
Curso	2024-25
Autores	César Abraham Vélez Rebollo Javier Martínez Hidalgo Héctor López de los Mozos Pérez Sandra Fuzhen Rodríguez Ibáñez
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

Estudiaremos una sección transversal de una columna recta que ocupa la región, en coordenadas cartesianas, [math][-1, 1] \times [0, 10][/math]. La variación de la temperatura (T) en cada punto de la columna viene definida por la función en coordenadas cilíndricas:

[math]T(\rho, \theta) = \sin(2 \pi \rho)[/math]

Los desplazamientos que sufre por la acción de una fuerza externa determinada, con el vector de posición

[math]\rho \hat{\mathbf{e}}_{\rho} + z \hat{\mathbf{e}}_z[/math]

son los siguientes:

[math]\vec{u}(\rho, \theta) = \rho \sin(\theta) \sin\left(\frac{2 \pi \rho}{50}\right) \, \hat{\mathbf{e}}_{\theta}[/math]

Tomaremos como función de densidad de los puntos de la columna la siguiente expresión:

[math]d(\rho, \theta) = (2 - \rho)\left(4 - \cos\left(4\left(\theta + \frac{\pi}{2}\right)\right)\right)[/math]

1 Mallado de los puntos interiores del sólido (Apartado 1)

Comenzaremos dibujando el mallado de nuestra columna; y para ello usaremos el siguiente código de Matlab:

%MALLADO DE LA COLUMNA
h=1/10;
x=-1:h:1;
y=0:h:10;
[Mx,My]=meshgrid(x,y);
mesh(Mx,My,Mx.*00);
title("Mallado de la columna")
view(2)
axis equal
axis tight
grid off
hold on
plot([-1,1],[0,0],'b')
plot([1,1],[0,10],'b')
plot([-1,1],[10,10],'b')
plot([-1,-1],[0,10],'b')
hold off

Mallado de los puntos de la columna

2 Curvas de nivel y gradiente de la temperatura [math]\nabla T[/math] (Apartado 2)

La temperatura viene definida por la siguiente función:

[math]T(\rho, \theta) = \sin(2 \pi \rho)[/math]

Pasada a coordenadas cartesianas sería:

[math]T(x,y) = \sin(2 \pi \sqrt{x^2 + y^2})[/math]

Con el siguiente código, se obtienen las curvas de nivel:

%GRADIENTE Y CURVAS DE NIVEL DE T
x=-1:h:1;
y=0:h:10;
figure
[Mx,My]=meshgrid(x,y);
mesh(Mx,My,Mx.*00);
Mz=sin(2*pi*sqrt((Mx.^2+My.^2))); %la función está pasada a cartesianas
contour(Mx,My,Mz)
colorbar
axis equal
axis tight
xlabel('Eje X')
ylabel('Eje Y')
title('Curvas de nivel de T')
figure
Mi=(2*pi*Mx.*cos(2*pi*sqrt(Mx.^2+My.^2)))./(sqrt(Mx.^2+My.^2));
Mj=(2*pi*My.*cos(2*pi*sqrt(Mx.^2+My.^2)))./(sqrt(Mx.^2+8My.^2));
quiver(Mx,My,-Mi,-Mj)
axis equal
axis tight
xlabel('Eje X')
ylabel('Eje Y')
title('Gradiente de T')
[m,n]=size(Mz);
maxi=0;
for i=1:m
    for j=1:n
        if Mz(i,j)>maxi
            maxi=Mz(i,j);
            pos1=i;
            pos2=j;
        end
    end
end
fprintf('La temperatura máxima es %f y se alcanza en las coordenadas %.2f, %.2f\n',maxi,Mx(pos1,pos2),My(pos1,pos2))

Curvas de nivel de la temperatura T

Gradiente de la temperatura [math]\nabla T[/math]

El valor del máximo, almacenado en la variable maxi, es 0.999989 y se alcanza en el punto de coordenadas (-0.90,8.20)

3 Flujo de energía calórica [math]\vec{Q}[/math] (Apartado 3)

Una vez tengamos el gradiente, es muy fácil calcular el flujo de energía calórica, puesto que el coeficiente de conductividad térmica de la placa es igual a 1:

[math]\vec{Q} = -k \nabla T[/math]

El gradiente de T lo obtenemos de derivar la función de T, en cartesianas, respecto de x e y. Nos queda como:

[math]\nabla T = \frac{2 \pi x \cos \left( 2 \pi \sqrt{x^2 + y^2} \right)}{\sqrt{x^2 + y^2}} \hat{i} + \frac{2 \pi y \cos \left( 2 \pi \sqrt{x^2 + y^2} \right)}{\sqrt{x^2 + y^2}} \hat{j}[/math]

Con el siguiente programa hallamos su representación:

%CAMPO VECTORIAL Q
h=1/10;
x=-1:h:1;
y=0:h:10;
[Mx,My]=meshgrid(x,y);
mesh(Mx,My,Mx.*00);
Mi=(2*pi*Mx.*cos(2*pi*sqrt(Mx.^2+My.^2)))./(sqrt(Mx.^2+My.^2));
Mj=(2*pi*My.*cos(2*pi*sqrt(Mx.^2+My.^2)))./(sqrt(Mx.^2+My.^2));
quiver(Mx,My,-Mi,-Mj)
hold on
Mz=sin(2*pi*sqrt(Mx.^2+My.^2)); 
contour(Mx,My,Mz)
axis equal
axis tight
title('Campo vectorial Q')
xlabel('Eje X')
ylabel('Eje Y')
hold off

Representación del campo vectorial [math]\vec{Q}[/math]

4 Puntos y dirección de la variación mayor de temperatura (Apartado 4)

En este apartado se ha determinado analíticamente en qué punto y en qué dirección la variación de temperatura es mayor, y se ha representado gráficamente en la placa mediante un conjunto de puntos rojos, que en este caso forman un semicírculo.

Cálculo del Gradiente de Temperatura

La temperatura está definida como:

[math]T(ρ,θ)=sin⁡(2πρ),[/math]

donde:

• [math] ρ=\sqrt{x^2 + y^2} [/math] es la distancia radial desde el origen.

• La temperatura no depende de [math]θ[/math], lo que simplifica los cálculos del gradiente.

El gradiente de temperatura se expresa como:

[math]\nabla T = \frac{\partial T}{\partial \rho} \hat{e}_\rho = 2\pi \cos(2\pi \rho) \hat{e}_\rho.[/math]

La magnitud del gradiente, que indica la variación de temperatura, es:

[math]|\nabla T| = |2\pi \cos(2\pi \rho)|.[/math]

La máxima variación ocurre cuando:

[math]|\cos(2\pi \rho)| = 1.[/math]

Esto se da para valores de [math]ρ=n/2, \in \mathbb{Z}.[/math] Dentro del dominio considerado [math](x∈[−1,1], y∈[0,10])[/math], el primer valor relevante es:

[math]\rho = 0.5.[/math]

Por tanto, los puntos donde la variación de temperatura es máxima forman un semicírculo de radio 0.5 centrado en el origen.

Dirección de Variación Máxima

El gradiente de temperatura apunta siempre en la dirección radial [math]\hat{e}_\rho[/math]. Dado que se considera únicamente el rango positivo de [math]cos⁡(2πρ)[/math], la dirección de máxima variación de temperatura corresponde a un vector alejándose del origen.

% Dominio de la placa
h = 0.1; % Paso de muestreo
x = -1:h:1; % Eje x
y = 0:h:10; % Eje y
[X, Y] = meshgrid(x, y);

% Cálculo de la temperatura
Rho = sqrt(X.^2 + Y.^2); % Distancia radial
T = sin(2*pi*Rho); % Temperatura en la placa

% Gradiente de temperatura
dT_drho = 2*pi*cos(2*pi*Rho);

% Puntos de máxima variación de temperatura (semicírculo para rho = 0.5)
rho_max = 0.5; % Radio del semicírculo
theta = linspace(0, pi, 100); % Valores de theta en [0, pi]
points_x = rho_max * cos(theta); 
points_y = rho_max * sin(theta); 

% Graficar la placa
figure;
contourf(X, Y, T, 20, 'LineColor', 'none'); % Contorno de temperatura
colorbar;
hold on;

% Representar gradiente con quiver
[U, V] = gradient(T, h); % Cálculo del gradiente en x, y
quiver(X, Y, U, V, 'k'); % Vectores de gradiente

% Dibujar el semicírculo de puntos de máxima variación
plot(points_x, points_y, 'ro', 'MarkerSize', 5, 'MarkerFaceColor', 'r'); % Puntos del semicírculo

% Personalización de la gráfica
title('Campo de Temperatura y Máxima Variación en la Columna');
xlabel('x');
ylabel('y');
axis equal;
xlim([-1, 1]);
ylim([0, 10]);
grid on;
hold off;

5 Campo de Desplazamientos en los puntos del mallado del sólido. (Apartado 5)

En este apartado se ha representado el campo de desplazamientos en los puntos del mallado del sólido y se ha identificado si existen puntos que permanecen fijos.

Cálculo del Campo de Desplazamientos.

El desplazamiento está definido en coordenadas cilíndricas como:

[math]\vec{u}(\rho, \theta) = \rho \sin(\theta) \sin\left(\frac{2\pi \rho}{50}\right) \hat{e}_\theta ,[/math]

donde:

• [math]\rho = \sqrt{x^2 + y^2}[/math] representa la distancia radial desde el origen.

• [math]\theta = \arctan\left(\frac{y}{x}\right)[/math] es el ángulo polar.

Para su representación en el plano cartesiano, el desplazamiento tangencial [math]\vec{u}_\theta[/math] se ha descompuesto en componentes [math]u_x[/math] y [math]u_y[/math] como:

[math] u_x = -u_\theta \sin(\theta),u_y = u_\theta \cos(\theta) [/math]

Se ha generado un mallado uniforme en el dominio [math][−1,1]×[0,10][/math], y el desplazamiento en cada punto del mallado se ha calculado y representado gráficamente.

Identificación de Puntos Fijos

Los puntos del sólido que permanecen fijos cuando [math]\vec{u} = 0[/math] , las condiciones bajo las cuales el desplazamiento se anula son las siguientes:

1. La magnitud del desplazamiento es nula en:

[math]\rho = 0 \quad \text{(el origen)} \quad \text{y cuando} \quad \sin(\theta) = 0.[/math]

• Para [math]\rho = 0[/math], el desplazamiento es nulo debido a la dependencia radial del término.

• Para [math]\sin(\theta) = 0[/math], el desplazamiento es nulo en las direcciones [math]\theta = 0[/math] y [math]\theta = \pi[/math], que corresponden al eje [math]x[/math].

Por tanto, los puntos fijos identificados son:

• El origen [math](x = 0, y = 0).[/math]

• Los puntos a lo largo del eje [math]x (y=0 [/math] o [math] y=10).[/math]

% Mallado del sólido
h = 0.1; % Paso del mallado
x = -1:h:1; % Eje x
y = 0:h:10; % Eje y
[X, Y] = meshgrid(x, y);

% Convertir a coordenadas cilíndricas
Rho = sqrt(X.^2 + Y.^2); % Distancia radial
Theta = atan2(Y, X); % Ángulo polar

% Cálculo del desplazamiento tangencial
U_theta = Rho .* sin(Theta) .* sin(2*pi*Rho/50);

% Convertir desplazamiento a coordenadas cartesianas
U_x = -U_theta .* sin(Theta); % Componente en x
U_y = U_theta .* cos(Theta); % Componente en y

% Dibujar el rectángulo (borde de la placa)
figure;
hold on;
rectangle('Position', [-1, 0, 2, 10], 'EdgeColor', 'b', 'LineWidth', 2); % Rectángulo

% Dibujar el campo de desplazamientos (flechas ampliadas)
scale = 2; % Escala de ampliación de las flechas
quiver(X, Y, U_x, U_y, scale, 'r'); % Campo de vectores con escala ajustada

% Personalización del gráfico
title('Campo de desplazamientos y bordes del sólido (flechas ampliadas)');
xlabel('x');
ylabel('y');
axis equal;
grid on;
hold off;

6 Sólido Antes y Después del Desplazamiento dado por el Campo de Vectores [math]\vec{u}[/math](Apartado 6)

En este apartado, se ha representado el sólido antes y después de aplicar el desplazamiento definido por el campo de vectores [math]\vec{u}[/math].

Cálculo del Desplazamiento

El desplazamiento está definido en coordenadas cilíndricas como:

[math]\vec{u}(\rho, \theta) = \rho \sin(\theta) \sin\left(\frac{2\pi \rho}{50}\right) \hat{e}_\theta,[/math]

donde:

• [math]\rho = \sqrt{x^2 + y^2}[/math]: es la distancia radial desde el origen.

• [math]\theta = \arctan\left(\frac{y}{x}\right)[/math]: es el ángulo polar.

• [math]\hat{e}_\theta[/math]: es la dirección tangencial en coordenadas cilíndricas.

Para representar este desplazamiento en coordenadas cartesianas, se descompone [math]\vec{u}(\rho, \theta)[/math] en las componentes [math]u_x[/math] y [math]u_y[/math] :

[math]u_x = -u_\theta \sin(\theta), \quad u_y = u_\theta \cos(\theta),[/math]

donde:

[math]u_\theta = \rho \sin(\theta) \sin\left(\frac{2\pi \rho}{50}\right).[/math]

De este modo, las nuevas posiciones de los puntos del sólido tras el desplazamiento son:

[math]x' = x + u_x, \quad y' = y + u_y.[/math]

Se ha generado un mallado uniforme en el dominio[math] [−1,1]×[0,10][/math], y las posiciones desplazadas de cada punto del mallado se han calculado.

% Paso del mallado
h = 0.1;
x = -1:h:1; % Eje x
y = 0:h:10; % Eje y
[X, Y] = meshgrid(x, y);

% Convertir a coordenadas cilíndricas
Rho = sqrt(X.^2 + Y.^2); % Distancia radial
Theta = atan2(Y, X); % Ángulo polar

% Cálculo del desplazamiento tangencial
U_theta = Rho .* sin(Theta) .* sin(2*pi*Rho/50);

% Convertir desplazamiento a componentes cartesianas
U_x = -U_theta .* sin(Theta);
U_y = U_theta .* cos(Theta);

% Nuevas posiciones tras el desplazamiento
X_new = X + U_x;
Y_new = Y + U_y;

% Crear la figura con subplot
figure;

% Subplot 1: Sólido inicial
subplot(1, 2, 1); 
mesh(X, Y, 0*X); % Representación del sólido inicial
title("Sólido inicial")
view(2)
axis equal
ylim([0, 11]);
grid off
grid on
hold on
% Dibujar bordes del sólido inicial
plot([-1, 1], [0, 0], 'b', 'LineWidth', 2); 
plot([1, 1], [0, 10], 'b', 'LineWidth', 2); 
plot([-1, 1], [10, 10], 'b', 'LineWidth', 2); 
plot([-1, -1], [0, 10], 'b', 'LineWidth', 2); 
hold off

% Subplot 2: Sólido desplazado y sin desplazar
subplot(1, 2, 2); % 1 fila, 2 columnas, segunda posición
hold on;
% Representar sólido sin desplazar (líneas azules)
plot(X, Y, 'k.', 'MarkerSize', 3); % Sólido sin desplazar como puntos negros
% Representar sólido desplazado con flechas
quiver(X, Y, U_x, U_y, 1.5, 'r', 'LineWidth', 1.2); % Flechas de desplazamiento en rojo
% Dibujar bordes del sólido sin desplazar
plot([-1, 1], [0, 0], 'b', 'LineWidth', 2); 
plot([1, 1], [0, 10], 'b', 'LineWidth', 2); 
plot([-1, 1], [10, 10], 'b', 'LineWidth', 2); 
plot([-1, -1], [0, 10], 'b', 'LineWidth', 2); 
% Personalización del gráfico
title('Sólido desplazado y sin desplazar');
xlabel('x');
ylabel('y');
ylim([0, 11]);
axis equal;
grid on;
hold off;

% Configurar título general
sgtitle('Comparación: Sólido Inicial y Sólido Desplazado');

7 Análisis y representación de [math]\nabla \cdot \vec{u}[/math] (Apartado 7)

La divergencia de un campo vectorial en cilíndricas viene definida por:

[math]\nabla\cdot\vec{u}=\frac{1}{ρ}\cdot[\frac{\partial}{\partial\rho}(\rho·\vec{ u_ρ })+\frac{\partial}{\partial{\theta}}({\vec u_\theta })+\frac{\partial}{\partial{z}}({\rho·\vec u_z })][/math]

La divergencia es un escalar que representa la diferencia entre el flujo saliente y el entrante de un campo vectorial. Particularizando para nuestro campo vectorial proporcionado en el enunciado [math]\vec{u}(ρ, θ) = \rho\cdot sin(\theta)\cdot sin(\frac{2\pi\rho}{50}) \vec{e_\theta} [/math], se opera y resulta:

[math]\nabla\cdot\vec{u}=\frac{1}{ρ}\cdot[\frac{\partial}{\partial\rho}(\rho·0)+\frac{\partial}{\partial{\theta}}({\rho\cdot sin(\theta)\cdot sin(\frac{2\pi\rho}{50})\vec u_\theta })+\frac{\partial}{\partial{z}}({\rho·0})][/math] [math]=\frac{1}{ρ}\cdot[\rho\cdot sin(\frac{2\pi\rho}{50})\cdot cos(\theta)][/math] [math]= sin(\frac{2\pi\rho}{50})\cdot cos(\theta)[/math]

Divergencia del campo u representada en 2D

Divergencia del campo u representada en 3D

h=2/10; %Coordenadas r y t
r=-1:h:1;
t=0:h:10;
%% Cálculo de coordenadas
[U,V]=meshgrid(r,t); %Generación de la retícula
X=U.*cos(V); %parametrización
Y=U.*sin(V);
mesh(U,V,0*U)
%% Divergencia
figure(2)
view(2)
axis([-1,1,0,10]);
div=sin(2*pi.*V/50).*cos(Y);
surf(U,V,div)
colorbar
axis([-1,1,0,10])
axis equal
title("Divergencia de la placa rectangular plana")
xlabel("Eje X")       
ylabel("Eje Y")
maximo=max(max(div))
minimo=min(min(div))

Se puede apreciar el cambio de volumen local debido al desplazamiento sobre todo en la gráfica 3D. Además, gracias a los últimos dos comandos podemos analizar los puntos donde la divergencia es máxima y mínima siendo:

[math]Divergencia_{máxima} = 0.9511[/math]
[math]Divergencia_{mínima} = 0[/math]

8 Análisis y representación de [math]|\nabla \times \vec{u}|[/math] (Apartado 8)

El rotacional de un campo vectorial en cilíndricas viene dado por:

[math]\nabla×\vec u(ρ,θ) = \frac{1}{ρ}\left|\begin{matrix} \vec e_ρ & \vec e_θ & \vec e_z \\ \frac{\partial}{\partial ρ} & \frac{\partial}{\partial θ} & \frac{\partial}{\partial z} \\ \vec u_ρ & \vec ρu_θ & \vec u_z \end{matrix}\right|[/math]

El resultado (que se trata de un vector), muestra la tendencia de un campo vectorial a inducir rotación alrededor de un punto. Particularizando a nuestro campo proporcionado por el enunciado [math]\vec{u}=sin(\theta)\cdot sin(\frac{2\pi\rho}{50})\vec{u_θ}[/math] resulta:

[math]=\frac{1}{ρ}\left|\begin{matrix} \vec e_ρ & \vec e_θ & \vec e_z \\ \frac{\partial}{\partial ρ} & \frac{\partial}{\partial θ} & \frac{\partial}{\partial z} \\ 0 & \rho\cdot [sin(\theta)\cdot sin(\frac{2\pi\rho}{50})]\vec{u_θ} & 0 \end{matrix}\right| = [(2\cdot sin(\theta)\cdot sin(\frac{2\pi\rho}{50}))+(\rho \cdot sin(\theta)\cdot \frac{2\pi\rho}{50} \cdot cos(\frac{2\pi\rho}{50}))]\vec{e_z}[/math].

Rotacional del campo u representado en 2D

Rotacional del campo u representado en 3D

h=2/10; %Coordenadas r y t
r=-1:h:1;
t=0:h:10;
%% Cálculo de coordenadas
[U,V]=meshgrid(r,t); %Generación de la retícula
X=U.*cos(V); %parametrización
Y=U.*sin(V);
mesh(U,V,0*U)
figure(2)
view(2)
axis([-1,1,0,10]);
axis equal
rot=(sin(V).*2.*sin(2*pi*U/50))+(sin(V).*U.*(2*pi/50).*cos(2*pi*U/50));
surf(U,V,rot)
colorbar
axis([-1,1,0,10])
axis equal
view(3)
[m,n]=size(rot);
maxi=0;
for i=1:m
     for j=1:n
         if rot(i,j)>maxi
             maxi=rot(i,j);
             pos1=i;
             pos2=j;
         end
     end
end
fprintf('El rotacional máximo es  %f y se alcanza en las coordenadas %.2f, %.2f\n',maxi,U(pos1,pos2),V(pos1,pos2))

El punto de mayor rotacional se puede sacar de forma analítica gracias al último comando. Siendo el resultado que nos da Matlab en la ventana de comandos: El rotacional máximo es 0.375179 y se alcanza en las coordenadas 1.00, 1.60 en [math] \rho ,\theta [/math] respectivamente.

9 Tensor de tensiones(Apartado 9)

En un medio elástico lineal, isótropo y homogéneo los desplazamientos permiten escribir el tensor de tensiones a través de [math]\sigma=\lambda\nabla\cdot\vec{u}1+2\mu\epsilon[/math]

Donde 1 es el tensor identidad, [math]\lambda=\mu=1[/math] , [math]\epsilon[/math] es el tensor deformaciones que viene dado por [math]\epsilon(\vec{u})=(\nabla\vec{u}+\nabla\vec{u}^t)/2[/math] y [math]\nabla\cdot\vec{u}[/math] es la divergencia del campo [math]\vec{u}[/math].

[math]\sigma=\begin{pmatrix} \sin(\frac{2\pi\rho}{50})\cdot cos(\theta) & 0 & 0\\ 0 & \sin(\frac{2\pi\rho}{50})\cdot cos(\theta) & 0\\ 0 & 0 & \sin(\frac{2\pi\rho}{50})\cdot cos(\theta) \end{pmatrix}+ 2 \begin{pmatrix} 0 & \frac{sin(\theta)}{2}\cdot (sin(\frac{2\pi\rho}{50})+\rho \cdot \frac{2\pi\rho}{50} \cdot (cos(\frac{2\pi\rho}{50})) & 0\\ \frac{sin(\theta)}{2}\cdot (sin(\frac{2\pi\rho}{50})+\rho \cdot \frac{2\pi\rho}{50} \cdot (cos(\frac{2\pi\rho}{50})) & sin(\frac{2\pi\rho}{50}) \cdot cos(\theta) & 0\\ 0 & 0 & 0 \end{pmatrix}[/math]
La matriz [math]\epsilon[/math] se saca tras operar lo siguiente: (siendo [math]\nabla\vec{u}[/math] el gradiente en cilíndricas del campo proporcionado) [math]= \frac{1}{2} \cdot [\begin{pmatrix}0 & 0 & 0\\ sin(\theta)\cdot (sin(\frac{2\pi\rho}{50})+\rho \cdot \frac{2\pi\rho}{50} \cdot (cos(\frac{2\pi\rho}{50})) & sin(\frac{2\pi\rho}{50}) \cdot cos(\theta) & 0\\ 0 & 0 & 0 \end{pmatrix}+\begin{pmatrix} 0 & sin(\theta)\cdot (sin(\frac{2\pi\rho}{50})+\rho \cdot \frac{2\pi\rho}{50} \cdot (cos(\frac{2\pi\rho}{50})) & 0\\ 0 & sin(\frac{2\pi\rho}{50}) \cdot cos(\theta) & 0\\ 0 & 0 & 0 \end{pmatrix}][/math]

Resultando que la matriz de tensiones es:

[math]\sigma=\lambda\nabla\cdot\vec{u}1+2\mu\epsilon=\begin{pmatrix}\sin(\frac{2\pi\rho}{50})\cdot cos(\theta) & sin(\theta)\cdot (sin(\frac{2\pi\rho}{50})+\rho \cdot \frac{2\pi\rho}{50} \cdot (cos(\frac{2\pi\rho}{50})) & 0\\ sin(\theta)\cdot (sin(\frac{2\pi\rho}{50})+\rho \cdot \frac{2\pi\rho}{50} \cdot (cos(\frac{2\pi\rho}{50})) & 3\cdot sin(\frac{2\pi\rho}{50})\cdot cos(\theta) & 0\\ 0 & 0 & \sin(\frac{2\pi\rho}{50})\cdot cos(\theta) \end{pmatrix}[/math]

A pesar de que los desplazamientos son planos, las tensiones no tienen por qué ser planas y puede haber tensiones en dirección ortogonal a la superficie de la placa. Para sacarlas, habiendo calculado la matriz [math]\sigma [/math] sólo debemos multiplicar por dicha dirección antes y después de la matriz:

Tensiones normales en la dirección [math]\vec{e_\rho}[/math]

[math]\vec{e_\rho}·\sigma·\vec{e_\rho}=(1,0,0)\cdot \begin{pmatrix}\sin(\frac{2\pi\rho}{50})\cdot cos(\theta) & sin(\theta)\cdot (sin(\frac{2\pi\rho}{50})+\rho \cdot \frac{2\pi\rho}{50} \cdot (cos(\frac{2\pi\rho}{50})) & 0\\ sin(\theta)\cdot (sin(\frac{2\pi\rho}{50})+\rho \cdot \frac{2\pi\rho}{50} \cdot (cos(\frac{2\pi\rho}{50})) & 3\cdot sin(\frac{2\pi\rho}{50})\cdot cos(\theta) & 0\\ 0 & 0 & \sin(\frac{2\pi\rho}{50})\cdot cos(\theta) \end{pmatrix} \cdot (1,0,0) [/math]
[math]\vec{e_\rho}·\sigma·\vec{e_\rho}=sin(\frac{2\pi\rho}{50})\cdot cos(\theta)[/math]

Tensiones normales en la dirección [math]\vec{e_\theta}[/math]

[math]\vec{e_\theta}·\sigma·\vec{e_\theta}=(0,1,0)\cdot \begin{pmatrix}\sin(\frac{2\pi\rho}{50})\cdot cos(\theta) & sin(\theta)\cdot (sin(\frac{2\pi\rho}{50})+\rho \cdot \frac{2\pi\rho}{50} \cdot (cos(\frac{2\pi\rho}{50})) & 0\\ sin(\theta)\cdot (sin(\frac{2\pi\rho}{50})+\rho \cdot \frac{2\pi\rho}{50} \cdot (cos(\frac{2\pi\rho}{50})) & 3\cdot sin(\frac{2\pi\rho}{50})\cdot cos(\theta) & 0\\ 0 & 0 & \sin(\frac{2\pi\rho}{50})\cdot cos(\theta) \end{pmatrix} \cdot (0,1,0) [/math]
[math]\vec{e_\theta}·\sigma·\vec{e_\theta}=3\cdot sin(\frac{2\pi\rho}{50})\cdot cos(\theta)[/math]

Tensiones normales en la dirección [math]\vec{e_z}[/math]

[math]\vec{e_z}·\sigma·\vec{e_z}=(0,0,1)\cdot \begin{pmatrix}\sin(\frac{2\pi\rho}{50})\cdot cos(\theta) & sin(\theta)\cdot (sin(\frac{2\pi\rho}{50})+\rho \cdot \frac{2\pi\rho}{50} \cdot (cos(\frac{2\pi\rho}{50})) & 0\\ sin(\theta)\cdot (sin(\frac{2\pi\rho}{50})+\rho \cdot \frac{2\pi\rho}{50} \cdot (cos(\frac{2\pi\rho}{50})) & 3\cdot sin(\frac{2\pi\rho}{50})\cdot cos(\theta) & 0\\ 0 & 0 & \sin(\frac{2\pi\rho}{50})\cdot cos(\theta) \end{pmatrix} \cdot (0,0,1) [/math]
[math]\vec{e_z}·\sigma·\vec{e_z}=sin(\frac{2\pi\rho}{50})\cdot cos(\theta)[/math]

Tensiones normales a cada dirección principal

h=2/10; %Coordenadas r y t
r=-1:h:1;
t=0:h:10;
%% Cálculo de coordenadas
[U,V]=meshgrid(r,t); %Generación de la retícula
X=U.*cos(V); %parametrización
Y=U.*sin(V);
mesh(U,V,0*U)
rho_sigma_rho=sin(2*pi*U/50).*cos(V);
theta_sigma_tetha=3*sin(2*pi*U/50).*cos(V);
z_sigma_z=sin(2*pi*U/50).*cos(V);
figure(2)
subplot(2,2,1)
surf(U,V,rho_sigma_rho)
view(3)
axis equal 
axis([-1,1,0,10]);
colorbar
title('i·\sigma·i')
subplot(2,2,2)
surf(U,V,theta_sigma_tetha)
view(3)
axis equal 
axis([-1,1,0,10])
colorbar
title('j·\sigma·j')
subplot(2,2,3.5)
surf(U,V,z_sigma_z)
view(3)
axis equal 
axis([-1,1,0,10])
colorbar
title('k·\sigma·k')

10 Tensiones tangeciales respecto al plano ortogonal a [math]\vec{i}[/math](Apartado 10)

Dichas tensiones vienen definidas por la expresión:

[math]|\sigma\cdot\vec{e_\rho}-(\vec{e_\rho}\cdot\sigma\cdot\vec{e_\rho})\vec{e_\rho}|[/math]

Puesto que [math](\vec{e_\rho}\cdot\sigma\cdot\vec{e_\rho})[/math] está calculado en el apartado anterior ([math]\vec{e_\rho}·\sigma·\vec{e_\rho}=sin(\frac{2\pi\rho}{50})\cdot cos(\theta)[/math]), sólo nos queda la primera parte: </math>.

[math]\sigma \cdot\vec{e_\rho}=\begin{pmatrix}\sin(\frac{2\pi\rho}{50})\cdot cos(\theta) & sin(\theta)\cdot (sin(\frac{2\pi\rho}{50})+\rho \cdot \frac{2\pi\rho}{50} \cdot (cos(\frac{2\pi\rho}{50})) & 0\\ sin(\theta)\cdot (sin(\frac{2\pi\rho}{50})+\rho \cdot \frac{2\pi\rho}{50} \cdot (cos(\frac{2\pi\rho}{50})) & 3\cdot sin(\frac{2\pi\rho}{50})\cdot cos(\theta) & 0\\ 0 & 0 & \sin(\frac{2\pi\rho}{50})\cdot cos(\theta) \end{pmatrix} \cdot (1,0,0)'[/math]
[math]=(sin(\frac{2\pi\rho}{50})\cdot cos(\theta),sin(\theta)\cdot (sin(\frac{2\pi\rho}{50})+\rho \cdot \frac{2\pi\rho}{50} \cdot (cos(\frac{2\pi\rho}{50})),0) [/math]

Resolviendo la resta de vectores:

[math]=(sin(\frac{2\pi\rho}{50})\cdot cos(\theta),sin(\theta)\cdot (sin(\frac{2\pi\rho}{50})+\rho \cdot \frac{2\pi\rho}{50} \cdot (cos(\frac{2\pi\rho}{50})),0)-(sin(\frac{2\pi\rho}{50})\cdot cos(\theta),0,0) [/math]
[math]|\sigma\cdot\vec{e_\rho}-(\vec{e_\rho}\cdot\sigma\cdot\vec{e_\rho})\vec{e_\rho}|=(0,sin(\theta)\cdot (sin(\frac{2\pi\rho}{50})+\rho \cdot \frac{2\pi\rho}{50} \cdot (cos(\frac{2\pi\rho}{50})),0)[/math]

Al ser las tensiones tangenciales verdaderamente pedidas son las del plano ortogonal a [math]\vec{i}[/math], pasamos a cartesianas nuestro resultado multiplicando por [math]\cdot(-sen(\theta))\vec{i}+cos(\theta)\vec{j})[/math] y cogemos únicamente la componente [math]\vec{i}[/math], que es la pedida.

[math]=[sin(\theta)\cdot (sin(\frac{2\pi\rho}{50})+\rho \cdot \frac{2\pi\rho}{50} \cdot (cos(\frac{2\pi\rho}{50}))]\cdot(-sen(\theta))\vec{i}+cos(\theta)\vec{j})[/math]

Tensiones tangenciales a la dirección principal pedida

%% Datos y Variables
h=2/10; %Coordenadas r y t
r=-1:h:1;
t=0:h:10;
%% Cálculo de coordenadas
[U,V]=meshgrid(r,t); %Generación de la retícula
X=U.*cos(V); %parametrización
Y=U.*sin(V);
mesh(U,V,0*U)


%%Tensiones tangenciales.
tan1=(sin(V).*(sin(2*pi*U/50)+U.*(2*pi.*U/50).*cos(2*pi.*U/50))).*(-sin(V));
tan2=(sin(V).*(sin(2*pi*U/50)+U.*(2*pi.*U/50).*cos(2*pi.*U/50))).*(-cos(V));
mod=sqrt(tan1.^2+tan2.^2);
figure(11)
view (3)
surf(U,V,mod)
axis equal 
axis([-1,1,0,10])
title('módulo de tensiones tangenciales respecto al plano ortogonal a i');

11 Tensión de Von Mises (Apartado 11)

La tensión de Von Mises, \( \sigma_{VM} \), se calcula utilizando la fórmula:

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

Donde:

- \( \sigma_1, \sigma_2, \sigma_3 \) son los autovalores (tensiones principales) del tensor de tensiones \( \sigma \).

- \( \sigma \) es una matriz \( 3 \times 3 \) que contiene las tensiones normales y de corte en las direcciones principales.

1. Definir el tensor de tensiones \( \sigma \):

- El tensor \( \sigma \) es una matriz que representa las tensiones en un punto. Se da o se calcula en función de las condiciones del problema. Aquí asumimos un ejemplo general:

[math]\sigma = \begin{bmatrix} \sigma_{xx} & \tau_{xy} & \tau_{xz} \\ \tau_{xy} & \sigma_{yy} & \tau_{yz} \\ \tau_{xz} & \tau_{yz} & \sigma_{zz} \end{bmatrix}. [/math]

2. Calcular los autovalores (\( \sigma_1, \sigma_2, \sigma_3 \)):

- Usamos el comando eig.m para calcular los autovalores de \( \sigma \).

3. Calcular \( \sigma_{VM} \):

- Sustituimos los autovalores en la fórmula de la tensión de Von Mises.

4. Identificar el punto de máxima tensión de Von Mises:

- Iteramos sobre los puntos del dominio, calculamos \( \sigma_{VM} \) en cada uno, y localizamos el valor máximo.

5. Visualizar la tensión de Von Mises:

- Graficar la tensión de Von Mises en el dominio.

- Marcar el punto de máxima tensión en el gráfico

Tensión de Von Mises en el dominio

>> % Dominio del sólido
h = 0.1; 
x = -1:h:1; 
y = 0:h:10; 
[X, Y] = meshgrid(x, y);

% Tensor de tensiones sigma (ejemplo simplificado)
sigma_xx = X.^2 + Y.^2; 
sigma_yy = X + Y;       
sigma_zz = X.^2 - Y.^2; 
tau_xy = X .* Y;        
tau_xz = zeros(size(X)); 
tau_yz = zeros(size(X)); 

% Matriz de tensiones
sigma_vm = zeros(size(X)); 

% Cálculo de σVM en cada punto
for i = 1:numel(X)
    % Tensor de tensiones en un punto
    sigma = [sigma_xx(i), tau_xy(i), tau_xz(i);
             tau_xy(i), sigma_yy(i), tau_yz(i);
             tau_xz(i), tau_yz(i), sigma_zz(i)];
    
    % Autovalores del tensor de tensiones
    eigenvalues = eig(sigma);
    sigma1 = eigenvalues(1);
    sigma2 = eigenvalues(2);
    sigma3 = eigenvalues(3);
    
    % Tensión de Von Mises
    sigma_vm(i) = sqrt(0.5 * ((sigma1 - sigma2)^2 + (sigma2 - sigma3)^2 + (sigma3 - sigma1)^2));
end

% Identificar el punto de máxima tensión de Von Mises
[max_vm, idx_max] = max(sigma_vm(:));
[x_max, y_max] = ind2sub(size(X), idx_max);

% Graficar tensión de Von Mises
figure;
contourf(X, Y, sigma_vm, 20, 'LineColor', 'none'); 
colorbar;
hold on;

% Marcar el punto de máxima tensión
plot(X(x_max, y_max), Y(x_max, y_max), 'ro', 'MarkerSize', 10, 'MarkerFaceColor', 'r');
title('Tensión de Von Mises y Punto de Máxima Tensión');
xlabel('x');
ylabel('y');
axis equal;
grid on;

12 Campo de fuerzas [math]\vec{F}[/math] que actúa sobre la placa(Apartado 12)

El campo de fuerzas \( \vec{F} \) que actúa sobre la placa se calcula como:

[math]\vec{F} = -\nabla \cdot \sigma[/math]

donde:

- \( \sigma \) es el tensor de tensiones, que en el Apartado 9 está definido como:

[math]\sigma = \lambda (\nabla \cdot \vec{u}) \mathbf{I} + 2\mu \varepsilon[/math]

con:

- \( \varepsilon = \frac{1}{2} (\nabla \vec{u} + (\nabla \vec{u})^T) \), el tensor de deformaciones.

- \( \mathbf{I} \), el tensor identidad.

- \( \lambda = 1 \) y \( \mu = 1 \), los coeficientes de Lamé dados.

La divergencia del tensor de tensiones \( \nabla \cdot \sigma \) se calcula fila a fila para obtener las componentes del campo de fuerzas.

Para resolver:

1. Definir el campo de desplazamientos \( \vec{u} = (u_x, u_y) \):

- \( \vec{u} \) es el campo que causa las deformaciones.

2. Calcular el tensor de deformaciones \( \varepsilon \):

- Calculamos las componentes de \( \varepsilon \) en términos de las derivadas parciales de \( u_x \) y \( u_y \).

3. Calcular el tensor de tensiones \( \sigma \):

- Sustituimos \( \nabla \cdot \vec{u} \) y \( \varepsilon \) en la expresión de \( \sigma \).

4. Divergencia del tensor de tensiones \( \nabla \cdot \sigma \):

- Calculamos las derivadas parciales de las filas de \( \sigma \).

5. Campo de fuerzas \( \vec{F} \):

- Calculamos \( -\nabla \cdot \sigma \).

6. Representación gráfica:

- Graficamos el campo de fuerzas \( \vec{F} \) sobre el dominio de la placa.

Campo de fuerzas [math]\vec{F}[/math] sobre el dominio de la placa

% Dominio de la placa
h = 0.1; 
x = -1:h:1; 
y = 0:h:10;
[X, Y] = meshgrid(x, y);

% Desplazamiento u
u_x = -X .* sin(Y); 
u_y = Y .* cos(X);  

% Gradiente del desplazamiento
[du_x_dx, du_x_dy] = gradient(u_x, h); 
[du_y_dx, du_y_dy] = gradient(u_y, h); 

% Tensor de deformaciones ε
epsilon_xx = du_x_dx; 
epsilon_yy = du_y_dy; 
epsilon_xy = 0.5 * (du_x_dy + du_y_dx);

% Divergencia de u
div_u = du_x_dx + du_y_dy;

% Tensor de tensiones σ
sigma_xx = div_u + 2 * epsilon_xx; 
sigma_yy = div_u + 2 * epsilon_yy; 
sigma_xy = 2 * epsilon_xy;         

% Divergencia del tensor de tensiones ∇·σ
[dsigma_xx_dx, ~] = gradient(sigma_xx, h);
[~, dsigma_xy_dy] = gradient(sigma_xy, h); 
[dsigma_xy_dx, ~] = gradient(sigma_xy, h); 
[~, dsigma_yy_dy] = gradient(sigma_yy, h); 

% Componentes del campo de fuerzas F
F_x = -(dsigma_xx_dx + dsigma_xy_dy); 
F_y = -(dsigma_xy_dx + dsigma_yy_dy); 

% Representación gráfica del campo de fuerzas
figure;
quiver(X, Y, F_x, F_y, 'b'); 
title('Campo de Fuerzas \bfF');
xlabel('x');
ylabel('y');
axis equal;
grid on;

13 Masa total(Apartado 13)

La masa total del sólido se calcula como:

[math]\text{Masa total} = \int_V d(x, y) \, dA[/math]

donde:

- \( d(x, y) = (2 - |x|)(4 - y) \) es la densidad de la placa.

- \( dA = dx \, dy \) es el diferencial de área en el plano.

El dominio de integración está definido en el problema como una placa rectangular:

[math]x \in [-1, 1], \quad y \in [0, 10][/math]

Aproximación numérica

La integral será aproximada usando una cuadrícula discreta en el dominio \((x, y)\), con un paso \(h\).

La suma de las densidades en cada punto multiplicada por el diferencial de área \((h^2)\) dará una estimación de la masa total.

Lo que da una masa total de 31.31

% Dominio de la placa
h = 0.1; 
x = -1:h:1; 
y = 0:h:10; 
[X, Y] = meshgrid(x, y);

% Densidad d(x, y)
density = (2 - abs(X)) .* (4 - Y);

% Cálculo de la masa total (suma numérica)
mass = sum(density(:)) * h^2;

% Mostrar el resultado
disp(['Masa total del sólido: ', num2str(mass)]);

% Representación gráfica de la densidad
figure;
contourf(X, Y, density, 20, 'LineColor', 'none');
colorbar;
title('Densidad d(x, y) sobre el sólido');
xlabel('x');
ylabel('y');
axis equal;
grid on;