Ondas y campos vectoriales en un arco (Grupo 32)

Trabajo realizado por estudiantes
Título	Ondas y campos vectoriales en un arco. Grupo 32
Asignatura	Teoría de Campos
Curso	2025-26
Autores	Isabella Eugenia Acosta Montoya Macarena Gil Andrés Guillermo Polo Toledo Marta de la Quintana Zubiría Liam O'Hea Kith
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 . Definición de las características del arco

1.1 El arco

La figura sobre la que se va a trabajar en este artículo es un arco o semi-anillo. Dicho arco está comprendido entre un radio mínimo de 1 y un radio máximo de 2. Por lo tanto, el dominio del arco queda comprendido dentro del siguiente intervalo:[math]\;1\le \sqrt[]{y^{2}+x^{2}}\le 2\;[/math], pero, al tratarse de un semi-anillo, queda restringido al semiplano superior, es decir, para [math]\;y\ge 0[/math]

Sabiendo las condiciones anteriores, se pueden expresar, también, las coordenadas cilíndricas que describen el arco. El radio queda comprendido entre 1 y 2, el ángulo limitado al semiplano superior y la altura libre:

• [math]\rho\in \left[ 1,2 \right][/math]
• [math]\theta\in \left[ 0,\pi \right][/math]

Por otra parte, se va a suponer que sobre el arco se ha aplicado una fuerza que ha provocado una vibración, de manera que los desplazamientos vienen dados por una onda, que depende de tres variables: dos espaciales y una temporal. Sin embargo, se va a poner el foco en las ondas transversales que vienen descritas, en este caso, por[math]\overrightarrow{u}(\rho,\theta)=\frac{1}{5}(\rho-1)\rho\;\overrightarrow{e}_{\rho}[/math]

1.2 Mallado del arco

A continuación, se ha generado el mallado de los puntos interiores del sólido. La imagen resultante se adjunta junto con el código, realizado en Matlab. Se ha tomado como paso de muestreo [math]h=\frac{1}{10}[/math]

Mallado del arco realizado a través de Matlab:

Dibujo del mallado del arco mediante MatLab.

% Paso de muestreo
h = 0.1;

% Valor de u (radios)
u = 1:h:2;

% Calculo de lo puntos que caben aproximadamente
puntos = round(pi/h) + 1;

% Vector v (ángulos) de 0 a pi
v = linspace(0, pi, puntos);

% Matrices de coordenadas polares
[rho, theta] = meshgrid(u, v);

% Conversión a Cartesianas
xx = rho .* cos(theta);
yy = rho .* sin(theta);


plot(xx, yy, 'r'); 

hold on %Se mantiene el gráfico para añadir las líneas radiales

plot(xx', yy', 'r');

hold off

% Ajustes finales visuales
axis equal
axis([-3, 3, -0.5, 2.5]);
title('Arco II');

2 . La temperatura

2.1 Representación de la temperatura del arco

La temperatura del sólido es una magnitud escalar, por lo que no tiene dirección ni sentido. En este caso, la temperatura del arco proviene de un foco de calor muy concentrado en los puntos que están a distancia 1 del origen. Dicha temperatura es una magnitud conocida y se distribuye a lo largo del arco siguiendo la siguiente función:[math]\;T(x,y)=(x-y)^{2} [/math]

Se va a ilustrar como se distribuye la temperatura utilizando curvas de nivel, que muestran los puntos que se encuentran a igual temperatura.
Se incluye una barra de colores que permite interpretar fácilmente los valores de temperatura a lo largo del sólido, asociando un color o escala de colores a una temperatura determinada.

Representación de la temperatura del arco

% Paso de muestreo
h = 0.1;
% Valor de u (radios)
u = 1:h:2;
% Cálculo de cuántos puntos caben aproximadamente
puntos = round(pi/h) + 1;
% Vector v (ángulos) de 0 a pi
v = linspace(0, pi, puntos);
% Matrices de coordenadas polares
[rho, theta] = meshgrid(u, v);
% Conversión a Cartesianas
xx = rho .* cos(theta);
yy = rho .* sin(theta);
% Expresión de la temperatura
T = (xx-yy).^2;

% Representación de la temperatura 
figure
contourf(xx,yy,T,10)
colorbar
title('Temperatura')
xlabel('x'); ylabel('y')
axis equal
axis([-3,3,-3,3])

2.2 Máximos y mínimos

Adicionalmente, se ha escrito el código para obtener la temperatura máxima y la temperatura mínima del arco.

% Encontrar el valor máximo y mínimo
% Puntos críticos
maxT = max(T(:));
minT = min(T(:));
fprintf('La temperatura máxima es %.2f ºC y la mínima es %.2f ºC\n',maxT,minT);

Este estudio de máximos y de mínimos también se puede realizar de manera anaítica.
Como se menciona en el principio, la temperatura viene dada por la función[math]\; T(x,y)=(x-y)^{2} [/math]. Se trata de una función cuadrática, por lo que siempre será mayor o igual que cero.
Para averiguar los valores mínimo y máximo se van a calcular y evaluar los puntos críticos de la función. Para hacer el estudio se va a hacer el paso de coordenadas cartesianas a coordenadas cilíndricas para facilitar los cálculos.

• [math]x=\rho cos\theta[/math]
• [math]y=\rho sen\theta[/math]

Para el cálculo de los puntos críticos, se van a igualar las derivadas parciales a cero, es decir:[math]\;\frac{\partial T}{\partial \rho}=0\;[/math] y [math]\;\frac{\partial T}{\partial \theta}=0[/math]:

• [math]\frac{\partial T}{\partial \rho}=2\rho(cos\theta-sen\theta)^{2}[/math]

Se iguala la expresión a cero: [math]\;2\rho(cos\theta-sen\theta)^{2}=0[/math]. Como [math]\rho\gt0\;[/math] la ecuación solo se satisface si[math]\;cos\theta=sen\theta\;[/math], es decir, si [math]\;\theta=\frac{\pi}{4}\;[/math]. Sustituyendo el valor obtenido en la función temperatura se observa que se obtiene la temperatura mínima, que es [math]T=0[/math]. Por lo tanto, en este caso, es independiente el valor de [math]\;rho\;[/math] para hallar el valor mínimo de la temperatura.

• [math]\frac{\partial T}{\partial \theta}=-2\rho^{2}(cos\theta^{2}-sen\theta^{2})=-2\rho^{2}cos2\theta[/math]

Se iguala la expresión a cero:[math]-2\rho^{2}cos2\theta=0[/math]. Como [math]\rho\gt0\;[/math] la ecuación solo se satisface si [math]\;cos2\theta=0[/math], es decir, si [math]\;\theta=\frac{3\pi}{4}\;[/math]o si [math]\;\theta=\frac{\pi}{4}\;[/math], que se descarta como candidato a máximo por el estudio anterior. Sustituyendo en la función temperatura: [math]\;T(\rho, \frac{3\pi}{4}) = \rho^2 (\cos(\frac{3\pi}{4}) - \sin(\frac{3\pi}{4}))^2 = \rho^2 (-\frac{\sqrt{2}}{2} - \frac{\sqrt{2}}{2})^2 = \rho^2 (-\sqrt{2})^2[/math]

El término [math]\;\rho^2\;[/math] representa la dependencia de la temperatura con la distancia al origen. Para maximizar[math]\;T[/math], el radio tiene que ser el máximo, es decir, 2. Sustituyendo en la función temperatura:[math]\;T(2,\frac{3\pi}{4})=2^2 (-\sqrt{2})^2=8\;[/math]. De esta manera, se obtendría el valor máximo.

Estos resultados se pueden apreciar de manera visual con la imagen que se ha adjuntado. Se puede observar que los tonos más fríos se concentran a lo largo de una línea recta: [math]\;x=y[/math] ([math]\theta = \pi/4[/math]). Esto demuestra que la temperatura mínima no depende de la distancia al origen [math]\;\rho\;[/math], sino únicamente de la condición angular[math]\;\cos\theta = \sin\theta[/math]. Por otra parte, los valores más altos de la temperatura están representados por colores más cálidos, que están situados en el borde exterior del arco ([math]\rho=2[/math]) en específico sobre la línea [math]y=-x[/math], es decir para [math]\;\theta = 3\pi/4[/math].

3 Cálculo y dibujo del gradiente de la Temperatura

El gradiente térmico es un vector que muestra cómo cambia la temperatura por unidad de distancia, indicando la dirección en que aumenta más rápidamente. Al estar expresada la temperatura en coordenadas cartesianas, se calcula el gradiente de la temperatura derivando parcialmente de x e y, obteniendo:

• [math]\nabla T(x,y)=(2(x-y),-2(x-y))[/math].

Al representarlo gráficamente, se observa que el gradiente [math]\nabla T[/math] es ortogonal a las curvas de nivel.

Representación de la temperatura del arco

% Paso de muestreo
h = 0.1;
% Valor de u (radios)
u = 1:h:2;
% Cálculo de cuántos puntos caben aproximadamente
puntos = round(pi/h) + 1;
% Vector v (ángulos) de 0 a pi
v = linspace(0, pi, puntos);
% Matrices de coordenadas polares
[rho, theta] = meshgrid(u, v);
% Conversión a Cartesianas
xx = rho .* cos(theta);
yy = rho .* sin(theta);
% Expresión de la temperatura
T = (xx - yy).^2;

% Gradiente
Tx =  2.*(xx-yy);     
Ty = -2.*(xx-yy);   

% Representación del gradiente
figure
contourf(xx,yy,T,20)    
hold on

% Curvas de nivel
contour(xx, yy, T, 20, 'k', 'LineWidth', 1.0)

quiver(xx,yy,Tx,Ty,'k')
hold off

axis equal
colorbar
title('Gradiente térmico del arco')
xlabel('x')
ylabel('y')

4 Dibujo del campo de vectores en los puntos del mallado del sólido

Texto en negrita

6 . Divergencia de [math]\overrightarrow{u}[/math]

La divergencia de un campo vectorial[math]\overrightarrow{u}[/math] definido en coordenadas cilíndricas se obtiene a partir de la siguiente expresión:

• [math]\nabla \cdot \overrightarrow{u}=\frac{1}{\rho}\left( \frac{\partial (\rho u_{\rho})}{\partial \rho}+\frac{\partial (u_{\theta})}{\partial \theta}+\frac{\partial(\rho u_{z}) }{\partial z} \right)[/math]

El campo vectorial con el que se va a operar es:[math]\;\overrightarrow{u}=\frac{1}{5}(\rho-1)\rho\;\overrightarrow{e}_{\rho}[/math]
Primeramente, se van a determinar las componentes conocidas, es decir,[math]\;u_{\rho}\;;\; u_{\theta}\;;\;u_{z}[/math] para así poder sustituir los valores en [math]\nabla \cdot \overrightarrow{u}[/math]

• [math]\;u_{\rho}=\frac{1}{5}(\rho^{2}-\rho)\;;\; u_{\theta}=0\;;\;u_{z}=0[/math]

• [math]\nabla \cdot \overrightarrow{u}=\frac{1}{\rho}\left( \frac{\partial (\rho u_{\rho})}{\partial \rho}+0+0 \right)=\frac{1}{\rho}\frac{\partial }{\partial \rho}\left( \rho\cdot \frac{1}{5}\cdot (\rho^{2}-\rho) \right)[/math]

Se multiplica[math]\;\;\overrightarrow{u}_{\rho}\;\;[/math] por[math]\;\;\rho\;\;[/math]y se deriva el término respecto de[math]\;\;\rho\;\;[/math], obteniendo:

• [math]\frac{\partial }{\partial \rho}\left( \frac{1}{5}(\rho^{3}-\rho^{2}) \right)=\frac{1}{5}(3\rho^{2}-2\rho)[/math]

Finalmente, se obtiene la divergencia:

• [math]\nabla \cdot \overrightarrow{u}=\frac{1}{\rho}\cdot\frac{1}{5}(3\rho^{2}-2\rho)=\frac{1}{5}(3\rho-2)[/math]

A continuación se adjunta el código desarrollado en Matlab para la representación de la divergencia.

Representación de la divergencia de [math]\overrightarrow{u}[/math] a través de matlab:

Dibujo de la divergencia en el arco mediante MatLab.

%Intervalos de las variables
w=10; 
p=80;
u=linspace(1,2,w);
v=linspace(0,pi,p); 

%mallado
[U,V] = meshgrid(u,v);
X = U.*cos(V);
Y = U.*sin(V);

%divergencia
DIVu = (3.*U - 2)./5;

%gráfica
surf(X,Y,DIVu)
view(2)
axis equal;
colorbar;
title('Divergencia en el arco');
axis([-3, 3, -0.5, 2.5]);

Observando la representación realizada en Matlab se puede apreciar que la variación es radial, los colores cambian progresivamente del centro hacia afuera. Esto sucede porque la fórmula de la divergencia depende linealmente de [math] \rho [/math].
También existe una simetría angular, es decir, si se fija un radio el color y valor de la divergencia será el mismo sin importar el ángulo [math] \theta [/math], esto sucede porque la función de la divergencia no depende de [math] \theta [/math].
En resumen, la imagen muestra que en esta región del espacio el campo vectorial [math]\overrightarrow{u}[/math] aumenta a medida que aumenta el radio sin importar la dirección.

7 . Rotacional de [math]\overrightarrow{u}[/math]

El rotacional de un campo vectorial [math]\overrightarrow{u}[/math] definido en coordenadas cilíndricas se calcula de manera genérica a partir de la siguiente expresión:

• [math] \nabla \times \overrightarrow{u}=\frac{1}{\rho} \begin{vmatrix} \overrightarrow{e_{\rho}} & \overrightarrow{e_{\theta}} & \overrightarrow{e_{z}}\\ \frac{\partial }{\partial \rho}& \frac{\partial }{\partial \theta}& \frac{\partial }{\partial z} \\ u_{\rho}& \rho u_{\theta} & u_{z} \end{vmatrix}[/math]

Realizando el cálculo para el campo vectorial [math]\overrightarrow{u}=\frac{1}{5}(\rho-1)\rho\overrightarrow{e}_{\rho}[/math]

• [math]\nabla \times \overrightarrow{u}=\frac{1}{\rho}·\begin{vmatrix} \overrightarrow{e_{\rho}} &\overrightarrow{e_{\theta}} & \overrightarrow{e_{z}}\\ \frac{2\rho-1}{5}& 0& 0 \\ \frac{\rho^{2}-\rho}{5} & 0 & 0 \end{vmatrix}=0\overrightarrow{e_{\rho}}+0\overrightarrow{e_{\theta}}+0\overrightarrow{e_{z}} [/math]

Como se puede observar el vector rotacional es el vector nulo, por lo que el campo es conservativo. Esto supone que el campo se comporta como un campo radial, es decir, no tiene tendencia a girar, sus líneas apuntan hacia o directamente desde un punto, pero no giran alrededor del centro.
No se representará el rotacional de [math]\overrightarrow{u}[/math] en matlab porque es 0.

8 . Tensor deformaciones

La parte simétrica del tensor gradiente es el tensor deformaciones, que viene descrito por la siguiente expresión:

• [math]\epsilon (\vec u) = \frac{\nabla \vec u + \nabla \vec u ^ t}{2}[/math]

En un medio elástico lineal, isótropo y homogéneo los desplazamientos permiten describir el tensor de tensiones [math] \sigma [/math] a través de la fórmula:

• [math]\sigma[/math]=[math]\lambda\nabla \cdot \overrightarrow{u}\;\mathbf{I}+2\mu\epsilon[/math]

Donde [math] \mathbf{I} [/math] es el tensor identidad en el conjunto de vectores libres del espacio [math] R^{3} [/math], y [math] \lambda [/math], [math] \mu [/math] son conocidos como coeficientes de Lamé, que dependen de las propiedades elásticas de cada material. Tomando [math] \lambda [/math][math]\;[/math]=[math]\;[/math][math] \mu [/math][math]\;[/math]=[math]\;[/math]1 se van a calcular y representar las tensiones normales que marcan el eje [math] \overrightarrow{e_{\rho}}[/math] y el eje [math]\frac{1}{\rho}\overrightarrow{e_{\theta}}[/math]
Para obtener dichas tensiones normales hay que realizar, anteriormente, unas operaciones. Se calculará el gradiente del campo vectorial [math]\overrightarrow{u}[/math] y su traspuesto, obteniendo así el tensor identidad:[math]\;[/math][math]\epsilon (\vec u) = \frac{\nabla \vec u + \nabla \vec u ^ t}{2}[/math].
Gradiente del campo [math]\overrightarrow{u}[/math]

• [math]\overrightarrow{u}=\frac{1}{5}(\rho-1)\rho\overrightarrow{e_{\rho}}[/math]

Se va a calcular dicho gradiente de manera matricial: [math]\left[ \nabla\overrightarrow{u}(\rho,\theta,z) \right]=\left( \frac{\partial \overrightarrow{u}}{\partial \rho}\left| \frac{1}{\rho}\frac{\partial \overrightarrow{u}}{\partial \theta} \right| \frac{\partial \overrightarrow{u}}{\partial z}\right)[/math]

• [math]\frac{\partial \overrightarrow{u}}{\partial \rho}=\frac{1}{5}(2\rho-1)\overrightarrow{e_{\rho}}+\frac{1}{5}(\rho-1)\rho\frac{\partial \overrightarrow{e_{\rho}}}{\partial \rho}=\frac{1}{5}(2\rho-1)\overrightarrow{e_{\rho}}[/math]

[math]\frac{\partial \overrightarrow{e_{\rho}}}{\partial \rho}=\Gamma^{k}_{11}\overrightarrow{e_{k}}=\Gamma^{1}_{11}\overrightarrow{e_{\rho}}+\Gamma^{2}_{11}\overrightarrow{e_{\theta}}+\Gamma^{3}_{11}\overrightarrow{e_{z}}=\overrightarrow{0}[/math]

• [math]\frac{\partial \overrightarrow{u}}{\partial \theta}=\frac{1}{5}(\rho-1)\rho\frac{\partial \overrightarrow{e_{\rho}}}{\partial \theta}=\frac{1}{5}(\rho-1)\rho\overrightarrow{e_{\theta}}[/math]

[math]\Gamma^{k}_{12}\overrightarrow{e_{k}}=\Gamma^{1}_{12}\overrightarrow{e_{\rho}}+\Gamma^{2}_{12}\overrightarrow{e_{\theta}}+\Gamma^{3}_{12}\overrightarrow{e_{z}}=1\cdot \overrightarrow{e_{\theta}}[/math]

• [math]\frac{\partial \overrightarrow{u}}{\partial z}=\frac{1}{5}(\rho-1)\rho\frac{\partial \overrightarrow{e_{\rho}}}{\partial z}=\overrightarrow{0}[/math]

[math]\frac{\partial \overrightarrow{e_{\rho}}}{\partial z}=\Gamma^{k}_{13}\overrightarrow{e_{k}}=\Gamma^{1}_{13}\overrightarrow{e_{\rho}}+\Gamma^{2}_{13}\overrightarrow{e_{\theta}}+\Gamma^{3}_{13}\overrightarrow{e_{z}}=\overrightarrow{0}[/math]

Por tanto:
[math]\left[ \nabla\overrightarrow{u}(\rho,\theta,z) \right]=\begin{pmatrix} \frac{1}{5}(2\rho-1)& 0&0 \\ 0 &\frac{1}{\rho}(\frac{1}{5}(\rho-1)\rho)&0 \\ 0 &0&0 \end{pmatrix}[/math][math]\quad[/math]=[math]\quad[/math][math]\frac{1}{5}\begin{pmatrix} 2\rho-1& 0&0 \\ 0 &\rho-1&0 \\ 0 &0&0 \end{pmatrix}[/math]

Cálculo de la matriz gradiente traspuesta

[math]\left[ \nabla\overrightarrow{u}(\rho,\theta,z) \right]^{t}=\frac{1}{5}\begin{pmatrix} 2\rho-1& 0&0 \\ 0 &\rho-1&0 \\ 0 &0&0 \end{pmatrix}[/math]

Como se puede observar[math]\quad[/math][math]\nabla \overrightarrow{u}=\nabla \vec u ^ t[/math]. Por lo tanto, el tensor deformaciones queda definido como:[math]\quad[/math][math]\epsilon (\vec u) = \frac{\nabla \vec u + \nabla \vec u }{2}[/math]
Se va a calcular el tensor deformaciones matricialmente

• [math]\epsilon (\vec u)=\frac{1}{2}\left( \frac{1}{5}\begin{pmatrix} 2\rho-1&0 &0 \\ 0& \rho-1&0 \\ 0& 0&0 \end{pmatrix}+\frac{1}{5}\begin{pmatrix} 2\rho-1 &0 &0 \\ 0& \rho-1&0 \\ 0&0 &0 \end{pmatrix} \right)=\frac{1}{5}\begin{pmatrix} 2\rho-1&0 &0 \\ 0& \rho-1&0 \\ 0& 0&0 \end{pmatrix}[/math]

Conociendo [math]\epsilon (\vec u)[/math], se obtiene [math]\sigma[/math]

• [math]\sigma=1\cdot \nabla \overrightarrow{u}\;\mathbf{I}+2\cdot 1\cdot \epsilon[/math]

• [math]\sigma= 1\cdot\frac{1}{5}\begin{pmatrix} 2\rho-1& 0&0 \\ 0 &\rho-1&0 \\ 0 &0&0 \end{pmatrix}\cdot \begin{pmatrix} 1 &0 &0 \\ 0& 1&0 \\ 0& 0&1 \end{pmatrix}+2\cdot 1\cdot\frac{1}{5} \begin{pmatrix} 2\rho-1& 0&0 \\ 0 &\rho-1&0 \\ 0 &0&0 \end{pmatrix}[/math]=[math]\;[/math][math]\frac{3}{5}\begin{pmatrix} 2\rho-1&0 &0 \\ 0& \rho-1&0 \\ 0& 0&0 \end{pmatrix}[/math]

Finalmente, con esta información, se puede proceder al cálculo de las tensiones normales, que son los valores de la diagonal principal de la matriz tensor de tensiones o [math]\sigma[/math]

• Tensión normal en la dirección del eje[math]\;\overrightarrow{e}_{\rho}:\overrightarrow{e}_{\rho}\cdot \sigma\cdot \overrightarrow{e}_{\rho}=\frac{3}{5}(2\rho-1)[/math]

• Tensión normal en la dirección del eje[math]\;\frac{1}{\rho}\;\overrightarrow{e}_{\theta}:\frac{1}{\rho}\;\overrightarrow{e}_{\theta}\cdot \sigma\cdot \frac{1}{\rho}\overrightarrow{e}_{\theta}=\frac{3}{5}(\rho-1)[/math]

9 Tensiones tangenciales respecto al plano ortogonal a [math]\vec e_\rho[/math]

En este apartado se calcularán las tensiones tangenciales respecto al plano ortogonal a [math]\overrightarrow{e}_{\rho}\;[/math], es decir

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

Haciendo referencia a los datos adquiridos en el apartado anterior:

• [math]\sigma=[/math][math]\frac{3}{5}\begin{pmatrix} 2\rho-1&0 &0 \\ 0& \rho-1&0 \\ 0& 0&0 \end{pmatrix}[/math]

• [math]\left(\vec e_{\rho}\cdot \sigma\cdot \vec e_{\rho}\right)=\frac{3}{5}(2\rho-1)[/math]

Se realizan las operaciones para obtener dichas tensiones tangenciales:

• [math]|σ·\vec e_ρ-(\vec e_ρ·σ·\vec e_ρ)·\vec e_ρ| = \left |\frac{3}{5}\begin{pmatrix} 2\rho-1&0 &0 \\ 0& \rho-1&0 \\ 0& 0&0 \end{pmatrix} ·\begin{pmatrix} 1\\0\\0 \end{pmatrix} - \left(\frac{3}{5}(2\rho-1)\right)\cdot\begin{pmatrix} 1\\0\\0 \end{pmatrix}\right|=0[/math]

Aunque las tensiones tangenciales respecto al plano ortogonal a [math]\overrightarrow{e}_{\rho}\;[/math] sean nulas no implica que no se puedan representar, lo que no será posible es ver una elevación si nos centramos en las tangenciales. Por otra parte las deformaciones producidas en el campo provienen totalmente de parte de las tensiones normales y por ello hay modificaciones en el arco plano inicial.

Representación de tensiones tangenciales

Mayor deformación del campo

% Paso de muestreo
h = 0.1;
% Valor de u (radios)
u = 1:h:2;
% Cálculo de cuántos puntos caben aproximadamente
puntos = round(pi/h) + 1;
% Vector v (ángulos) de 0 a pi
v = linspace(0, pi, puntos);

[rho, theta] = meshgrid(u, v);

%Cambio a coordendas cartesianas
xx = rho .* cos(theta);
yy = rho .* sin(theta);

%Cálculo del vector radial
erx = cos(theta);
ery = sin(theta);
erz = zeros(size(theta));

%sigma
sigma_xx = (3/5) * (2 * ones(size(rho)));
sigma_xy = (3/5) * rho;
sigma_xz = (3/5) * (-1 * ones(size(rho)));

sigma_yx = zeros(size(rho));
sigma_yy = zeros(size(rho));
sigma_yz = zeros(size(rho));

sigma_zx = zeros(size(rho));
sigma_zy = zeros(size(rho));
sigma_zz = zeros(size(rho));

%Vector tensión con la componente normal y la tangente
Tx = sigma_xx .* erx + sigma_xy .* ery + sigma_xz .* erz;
Ty = sigma_yx .* erx + sigma_yy .* ery + sigma_yz .* erz;
Tz = sigma_zx .* erx + sigma_zy .* ery + sigma_zz .* erz;

Tn_scalar = Tx .* erx + Ty .* ery + Tz .* erz;   % = 3/5 (2ρ - 1)

Tnx = Tn_scalar .* erx;
Tny = Tn_scalar .* ery;
Tnz = Tn_scalar .* erz;

Ttx = Tx - Tnx;
Tty = Ty - Tny;
Ttz = Tz - Tnz;

%Dibujo del campo de tensiones tangenciales
figure;
quiver(xx, yy, Ttx, Tty, 'LineWidth', 1.3);
axis equal
title('Tensiones tangenciales')
xlabel('x')
ylabel('y')

%Mayor deformación del campo
figure
surf(xx, yy, Tn_scalar); 
shading interp;
colorbar;
title('Mayor deformación del campo');
xlabel('x'); ylabel('y');

10 . Tensiones tangenciales respecto al plano ortogonal a [math]\frac{1}{\rho}[/math][math]\vec e_\theta[/math]

En este apartado se calcularán las tensiones tangenciales respecto al plano ortogonal a [math]\frac{1}{\rho}[/math][math]\vec e_\theta[/math] , es decir:

• [math]\left |\sigma\cdot\frac{1}{\rho}\vec e_\theta-(\frac{1}{\rho}\vec e_\theta\cdot \sigma \cdot\frac{1}{\rho}\vec e_\theta) \frac{1}{\rho}\vec e_\theta \right |[/math]

Haciendo referencia a los datos adquiridos en el apartado 8:

• [math]\sigma=[/math][math]\frac{3}{5}\begin{pmatrix} 2\rho-1&0 &0 \\ 0& \rho-1&0 \\ 0& 0&0 \end{pmatrix}[/math]

• [math]\left(\frac{1}{\rho}\vec e_{\theta}\cdot \sigma\cdot \frac{1}{\rho}\vec e_{\theta}\right)=\frac{3}{5}(\rho-1)[/math]

Se obtienen, entonces, dichas tensiones tangenciales:

• [math]\left |\sigma\cdot\frac{1}{\rho}\vec e_\theta-(\frac{1}{\rho}\vec e_\theta\cdot \sigma \cdot\frac{1}{\rho}\vec e_\theta) \frac{1}{\rho}\vec e_\theta \right |= \left |\frac{3}{5}\begin{pmatrix} 2\rho-1&0 &0 \\ 0& \rho-1&0 \\ 0& 0&0 \end{pmatrix} ·\begin{pmatrix} 0\\\frac{1}{\rho}\\0 \end{pmatrix} - \left(\frac{3}{5}(\rho-1)\right)\cdot\begin{pmatrix} 0\\\frac{1}{\rho}\\0 \end{pmatrix}\right|= 0 [/math]

Al igual que en el apartado anterior, las tensiones tangenciales son cero por lo que no hay puntos en los que sea mayor. Esto no quita que no sea posible representarlas. Por otra parte, si hay mayor o menor deformación en el campo, aunque depende completamente de las tensiones normales

Representación de tensiones tangenciales

Mayor deformación del campo

% Paso de muestreo
h = 0.1;
% Valor de u (radios)
u = 1:h:2;
% Cálculo de cuántos puntos caben aproximadamente
puntos = round(pi/h) + 1;
% Vector v (ángulos) de 0 a pi
v = linspace(0, pi, puntos);

[rho, theta] = meshgrid(u, v);

%Cambio a coordenadas cartesianas
xx = rho .* cos(theta);
yy = rho .* sin(theta);

%Cálculo del vector radial
erx = cos(theta);
ery = sin(theta);
erz = zeros(size(theta));

%Sigma
sigma_xx = (3/5) * (2*rho - 1);
sigma_xy = zeros(size(rho));
sigma_xz = zeros(size(rho));

sigma_yx = zeros(size(rho));
sigma_yy = (3/5) * rho;
sigma_yz = zeros(size(rho));

sigma_zx = zeros(size(rho));
sigma_zy = zeros(size(rho));
sigma_zz = zeros(size(rho));

%Vector tension con la componente normal y tangencial
Tx = sigma_xx .* erx + sigma_xy .* ery + sigma_xz .* erz;
Ty = sigma_yx .* erx + sigma_yy .* ery + sigma_yz .* erz;
Tz = sigma_zx .* erx + sigma_zy .* ery + sigma_zz .* erz;

Tn_scalar = Tx .* erx + Ty .* ery + Tz .* erz;

Tnx = Tn_scalar .* erx;
Tny = Tn_scalar .* ery;
Tnz = Tn_scalar .* erz;

Ttx = Tx - Tnx;
Tty = Ty - Tny;
Ttz = Tz - Tnz;

%Dibujo del campo de tensiones tangenciales
figure;
quiver(xx, yy, Ttx, Tty, 'LineWidth', 1.3);
axis equal
title('Tensiones tangenciales')
xlabel('x')
ylabel('y')

%Representación de la mayor deformación del campo
figure;
surf(xx, yy, Tn_scalar);
shading interp;
colorbar;
title('Mayor deformación del campo');
xlabel('x'); ylabel('y');

11 Masa aproximando la integral numéricamente

La densidad de la placa viene dada por la función [math] d(\rho,\theta)=1 + e^{\rho^2 \cos\theta} [/math]. A partir de este dato se va a calcular la masa del arco. La coordenada [math] \rho [/math] conservará el radio del arco y [math] \theta [/math] es un ángulo que queda comprendido entre 0 y [math] \pi [/math], por tanto: [math] (\rho,\theta)=\left[ 1;2 \right]\times \left[ 0,\pi \right] [/math]
La integral quedaría expresada de la siguiente manera:
M = [math] \int_{0}^{\pi} \int_{1}^{2} \left(1 + e^{\rho^2 \cos \theta}\right) \rho \, d\rho \, d\theta [/math]
Sin embargo, se va resolver de manera numérica con Matlab.
Masa del arco utilizando métodos numéricos (Método Trapecio):

% MATLAB Code para calcular la masa 
% Densidad: d(rho, theta) = 1 + exp(rho^2 * cos(theta))
% Región: 1 <= rho <= 2, 0 <= theta <= pi
% Parámetros y Límites de integración
rho_min = 1;
rho_max = 2;
theta_min = 0;
theta_max = pi;
% Subintervalos
N_rho = 50;  
N_theta = 100; 
% Tamaño del paso (h) para cada variable
h_rho = (rho_max - rho_min) / N_rho;
h_theta = (theta_max - theta_min) / N_theta;
% Mallado de la superficie 
rho_vec = linspace(rho_min, rho_max, N_rho + 1);
theta_vec = linspace(theta_min, theta_max, N_theta + 1);
[RHO, THETA] = meshgrid(rho_vec, theta_vec);
% Función densidad 
D = RHO .* (1 + exp(RHO.^2 .* cos(THETA)));
% Cálculo de los Pesos de la Regla del Trapecio
% Pesos para la dirección rho
w_rho = ones(N_rho + 1, 1);
w_rho(1) = 0.5;
% Pesos para la dirección theta 
w_theta = ones(N_theta + 1, 1);
w_theta(1) = 0.5;
w_theta(end) = 0.5;
% Fórmula Matricial de la Integral Doble (Método del Trapecio)
M = h_rho * h_theta * (w_theta') * F * w_rho;
% Resultado de la masa 
fprintf('Masa (M) aproximada: %.5f\n', M)