Campos de temperaturas y deformaciones en superficie plana semicircular (Grupo 23)
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Campos de temperaturas y deformaciones en superficie plana semicircular (Grupo 23) |
| Asignatura | Teoría de Campos |
| Curso | 2025-26 |
| Autores |
|
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
Un Campo de Temperaturas es un tipo de Campo Escalar definido sobre un objeto físico para describir la temperatura existente en cada uno de sus puntos.
La visualización de un campo escalar puede realizarse con herramientas como Octave y MATLAB, que emplean operaciones entre escalares y matrices para llevar a cabo cálculos de mayor complejidad.
Enlace de descarga al póster del trabajo en formato PDF: Póster Grupo 23
Contenido
- 1 Creación de la superficie 2D
- 2 Representación de Campos Escalares
- 3 Calculo del ∇T y su representación
- 4 Representación del campo de vectores en los puntos del mallado del sólido
- 5 Representación del desplazamiento del campo ū
- 6 Calculo de la ∇*u y su representación
- 7 Cálculo del Rotacional del campo u
- 8 Cálculo de tensiones normales
- 9 Cálculo de tensiones tangenciales
- 10 Cálculo de las tensiones tangenciales respecto al plano ortogonal a [math]\vec{e_\theta}[/math]
- 11 Cálculo de la masa de un campo en coordenadas polares cilíndricas
- 12 Aplicaciones reales en ingeniería
- 13 Póster del trabajo
- 14 Bibliografía
1 Creación de la superficie 2D
Para poder mostrar el campo de temperaturas de un sólido en MATLAB, lo primero es construir una malla que represente la geometría del objeto. Esto se hace creando una red de puntos mediante la instrucción meshgrid, propio de MATLAB.
Pensemos ahora en una discretización que modele los puntos internos de un sólido rectangular bidimensional, situado en el plano cartesiano XY (en este caso en coordenadas polares), utilizando un paso de muestreo [math](h) = 1/10[/math].
En este caso, el arco esta comprendido entre los radios para: [math]p=[1,2][/math], y el angulo en XY, [math]𝜃=[0,π][/math].
El código MATALB será el siguiente:
%Primero definimos las variables que vamos a utilizar (ρ,θ) y h (paso de
%muestreo)
h = 0.1;
r = (1:h:2);
angulo = (0:h:pi);
%Definimos los ejes, ajustandolos ambos a un tamaño equitativo para que sea
%proporcional la figura y establecemos con el comando view que se vea en 2D
axis equal;
axis([-3,3,-1.5,3]);
view(2);
hold on;
%Para dibujar la gráfico establecemos RR y AA como variables con meshgrid
[RR, AA] = meshgrid(r, angulo);
x = RR.*cos(AA);
y = RR.*sin(AA);
%Para evitar que aparezca el eje z reflejado ponemos 0.*x
mesh(x,y,0.*x);
%Añadimos finalmente información a los ejes y reflejamos la cuadrícula
xlabel('Eje X');
ylabel('Eje Y');
grid on;
hold off;
2 Representación de Campos Escalares
El campo escalar de temperaturas de un objeto físico en 2D se obtiene designando un valor por medio una función a cada punto del mallado que previamente se ha creado para representar el objeto. En este caso se ha generado una funcion dada en coordeandas polares [math]p[/math],[math]𝜃[/math]
%Mantemeos parte del apartado 1 solo que en vez de definir la malla del dominio nos centramos en representar la temperatura.
h = 0.1;
r = (1:h:2); angulo = (0:h:pi); %Parametrización
[RR, AA] = meshgrid(r, angulo); %Creacion del mallado
x = RR .* cos(AA); y = RR .* sin(AA);
axis equal; axis([-3,3,-1.5,3]);
view(2); hold on;
%Definimos la función Temperatura
T = (x - y).^2;
%Dibujo de la función T aplicada al dominio definido en el apartado 1.
contourf(x, y, T, 40);
xlabel('Eje X'); ylabel('Eje Y');
grid on; hold off3 Calculo del ∇T y su representación
Para el cálculo del gradiennte de un campo escalar dado en coordenadas polares (cilíndicas) sea
Al no haber componente [math]ez[/math] esta se anula, con el gradiente ya calculado se dibuja el campo vectorial. Este se realiza mediante la función quiver. Que utiliza cuatro variables. Las dos primeras (X;Y) serán desde donde va a partir el campo vectorial y las dos siguientes. (Gradientx, Gradiente) hacia donde va el vector.
%Mantenemos parte de la información del apartado 1 y 2
h = 0.1; r = (1:h:2);
angulo = (0:h:pi);
[RR, AA] = meshgrid(r, angulo);
x = RR .* cos(AA); y = RR .* sin(AA);
T = (x - y).^2;
axis equal; axis([-3,3,-1.5,3]);
hold on;
contourf(x, y, T, 40);
%Gradiente de T
[derivadaradio, derivadaangulo] = gradient(T, h, h);
%Como estamos en polares se tiene que dividir la derivada del angulo entre el radio
derivadaangulo = derivadaangulo ./ RR;
%Se convierten las coordenadas de polares a cartesianas para usar quiver
ux = derivadaradio .* cos(AA) - derivadaangulo .* sin(AA);
uy = derivadaradio .* sin(AA) + derivadaangulo .* cos(AA);
quiver(x, y, ux, uy, 'r');
xlabel('Eje X'); ylabel('Eje Y');
grid on; hold off4 Representación del campo de vectores en los puntos del mallado del sólido
El campo de vectores del mallado en este caso dado el vector
h = 0.1; r = (1:h:2); angulo=(0:h:pi); %Parametrización
[RR, AA]=meshgrid(r, angulo);
x=RR.*cos(M);
y=RR.*sin(M);
hold on; axis equal; axis([-3,3,-1.5,3]);
u=(1/5) (RR-1).* RR;
%Se calcula tanto ux como uy(1,3) de u tras multiplicarse por ep;
ux=u.*cos(AA);
uy=u.* sin(AA);
quiver(x, y, ux, uy, 'black');%Se dibuje la gráfica
xlabel("Eje= X"); ylabel("Eje Y");
grid off; hold off;
5 Representación del desplazamiento del campo ū
En este caso se nos pide comparar la representación del campo escalar ū, tras sufrir un desplazamineto, a través de la función subplot en MATLAB Se define la nueva posición de las coordenadas XX e YY sumándole dicho desplazamiento
h=0.1; r=(1:h:2);
angulo=(0:h:pi);
[RR, AA]=meshgrid(r, angulo); %Creación del mallado
x =RR.*cos(A); y=RA.*sin(AA);
u=(1/5)*(RR-1). *RR; %Campo escalar
ux=u.*cos(AA); uy=u.*sin(AA);
%Definimos la nueva posición de las x mediante la suma de la original más el desplazamiento que sufre este.
xnuevo=x+ux; ynuevo=y+uy;
subplot(1,2,1) %Dibujenos priser el sólido original
mesh(xnuevo,ynuevo,0.*x);
axis equal; axis (-1,3,-1.5,3]);
xlabel('Eje X'); ylabel('Eje Y');
view(2)
6 Calculo de la ∇*u y su representación
Dado el campo dado [math]u(p,𝜃)=1/5*(p − 1)p𝑒p[/math] , dado que se encuentre en coordenadas cilíndricas para el cálculo de la divergencia, se utilizara la siguiente fórmula:
En este caso la componente ɞz se anula al no haber en el campo ū
A través de cálculos matemáticas sea la [math]∇*u= 1/5*(3p-2)[/math] Sea creciente con p, por tanto los puntos donde alcance una mayor divergencia sea en el borde 'exterior' con p=2
6.1 La representación de la divergencia del campo:
h=0.1; r=(1:h:2);
angulo=(0:h:pi);
[RR, AA]=meshgrid(r, angulo); %Creación del mallado
x =RR.*cos(A); y=RA.*sin(AA);
u=(1/5)*(RR-1). *RR; %Campo escalar
%Se calcula la derivada de u respecto al radio. No hace falta calcularla
respecto al ángulo ya que como no esta incluida en la u original va a ser
[derivadau,valor_inecesario] = gradient(u, h, h);
%Calculamos la divergencia
divergencia=derivadau+u./ RR;
%Dibujanos la divergencia una vez que esta ya está calculada
contourf(x, y, divergencia, 20);
axis equal; axis([-3,3,-1.5,3]);
xlabel('Eje X'), ylabel('Eje Y'); grid on;
7 Cálculo del Rotacional del campo u
Sea ∇×u el rotacional de un campo de desplazamientos expresado en un sistema de coordenadas en vectores unitarios ortogonales , dado nuestro campo en coordenadas cilíndricas se verifica que el rotacional del campo 'u se calcule partir de la fórmula: [math]\ \bigtriangledown \times \vec{u} = 1/p*\begin{vmatrix}\vec{ep} & \vec{peθ} &\vec{ez} \\ d/dp & d/dθ & d/dz\\ Fp & pFθ & Fz \end{vmatrix}[/math]
En este caso el cálculo del ∇×u=0 por tanto su valor absoluto también será cero. Debido a que en nuestro campo u no hay tanto componente eθ y ez, se anula las derivadas parciales.
7.1 Representación del valor absoluto del rotacional
%Mantenemos parte de la información de apartados anteriores
h = 0.1;
r = (1:h:2);
angulo = (0:h:pi);
[RR, AA] = meshgrid(r, angulo);
x = RR .* cos(AA);
y = RR .* sin(AA);
%Apartado 7
%Dividimos u respecto de r y el ángulo. Teniendo en cuenta que del ángulo
%es inexistente
urespectoradio = (1/5) * (RR - 1) .* RR;
urespectoangulo = 0;
%Calculamos las derivadas parciales
[derivadauradio, derivadauangulo] = gradient(urespectoradio, h, h);
%Calculamos el rotacional
rotacional = (1 ./ RR) .* derivadauangulo.* urespectoangulo - (1 ./ RR) .* derivadauangulo .* urespectoangulo;
%Como el enunciado pide el valor absoluto del rotacional y no el rotacional
rotacionalabs = abs(rotacional);
%Graficamos el rotacional
contourf(x, y, rotacionalabs, 20);
axis equal;
axis([-3,3,-1.5,3]);
xlabel('Eje X'); ylabel('Eje Y');
title('Rotacional en Valor Absoluto');
grid on;
colorbar;8 Cálculo de tensiones normales
El tensor de tensiones asociado a un medio elástico, isótropo y homogéneo, donde [math]e=\frac{\bigtriangledown\vec{u}+\bigtriangledown\vec{u}^t}{2}[/math] representa la parte simétrica del gradiente del campo vectorial [math]\vec u[/math]. Tomando los coeficientes de Lamé [math]λ,µ =1[/math] y considerando que [math]\bigtriangledown \cdot \vec{u}[/math] es la divergencia del campo [math]\vec u[/math], es posible obtener las deformaciones proyectadas sobre los vectores de la base ortonormal orientada positivamente [math]{\vec i,\vec j,\vec k}[/math].
8.1 Cálculo del gradiente del campo [math]\vec{u}[/math]
En este caso dado el campo, eθ y ez sean nulos. En este caso para poder calcular el gradiente de un campo vectorial en este caso, en coordenadas cilíndricas sea a partir de esta fórmula.
∇u =\begin{pmatrix}dup/dp & 1/p*dup/dθ-uθ/p\\duθ/dp & 1/p*duθ/dθ+up/p\end{pmatrix}
De tal manera que ∇u =\begin{pmatrix}1/5*2p-1 & 0\\0 & 1/5*p-1\end{pmatrix} Sea la traspuesta de ∇u igual a ∇u , por tanto sea el tensor de deformación simétrico sea:
Considerando el tensor de deformación , siendo I, la matriz identidad. [math]sigma=λ\bigtriangledown \cdot \vec{u}1 + 2µe[/math]
8.2 Cálculo del Tensor de deformaciones [math]\sigma[/math]
Sea en el apartado 6, la divergencia sea [math]∇*u= 1/5*(3p-2)[/math], por tanto considerando que [math]λ,µ =1[/math], y sustituyendo la parte simétrica del gradiente sea
8.3 Representación de las tensiones normales
Dibujar las tensiones normales en la dirección de ep y 1/p*eθ De tal manero que sea:
8.3.1 Tensor normal radial
[math]ơ radial=ep*ơ*ep[/math], sustituyendo los valores de ơ siendo ep=(1,0) esto se debe que al trabajar en coordenadas cilíndricas, la base física deba ser ortogonal y unitaria, considerando la base de referencia [math]{\vec i,\vec j,\vec k}[/math].
8.3.2 Tensor normal tangencial
[math]ơtangencial=1/peθ*ơ*1/peθ[/math],sustituyendo los valores de ơ siendo eθ=(0,1)
Mediante la representación de MATLAB:
%Mantenemos parte de la información de apartados anteriores
h = 0.1; r = (1:h:2);
angulo = (0:h:pi);
[RR, AA] = meshgrid(r, angulo);
x = RR .* cos(AA); y = RR .* sin(AA);
%Apartado 8
%Calculamos las dos tensiones normales
tensionnormal1 = (7*RR - 4)/5;
tensionnormal2 = (5*RR - 4) ./ (5 * RR.^2);
%Graficamos ambas funciones usando subplot
subplot(1,2,1);
contourf(x, y, tensionnormal1, 20);
axis equal;
xlabel('Eje X'); ylabel('Eje Y');
title('Tensión normal radial');
colorbar; grid on;
%Ahora vamos con la segunda
subplot(1,2,2); contourf(x, y, tensionnormal2, 20);
axis equal;9 Cálculo de tensiones tangenciales
Una vez calculadas las tensiones normales, nos interesa también calcular la tensión en su componente tangencial, es decir, aquella cuya dirección es paralela a la superficie del plano. Este dato es muy importante para saber qué partes de la superficie sufren más tensiones y así poder anticiparse a posibles fallos en las estructuras debidas a sobreesfuerzos de tensión.
9.1 Respecto al plano ortogonal a [math]\vec{e_\rho}[/math]
Seguimos trabajando sobre la misma superficie plana:
[math]\tau_{\rho} = \left|\, \sigma \cdot \vec{e}_{\rho} \;-\; \left( \vec{e}_{\rho} \cdot \sigma \cdot \vec{e}_{\rho} \right)\vec{e}_{\rho} \,\right|[/math]
Además, el tensor de deformaciones [math]\sigma[/math] ya lo conocemos pues fue calculado en el apartado anterior:
Para simplificar la operación, la ejecutamos en varias partes. Primero calculamos el producto [math]\sigma·\vec{e_\rho}[/math], donde conocemos el vector [math]\vec{e_\rho}=\begin {pmatrix}1\\0\end{pmatrix}[/math]. Por tanto:
[math]\sigma \cdot \vec{e}_{\rho}
=
\begin{pmatrix}
\sigma_{\rho\rho} & 0 \\
0 & \sigma_{\theta\theta}
\end{pmatrix}
\begin{pmatrix}
1 \\ 0
\end{pmatrix}
=\begin{pmatrix}\dfrac{7\rho - 4}{5} \\[6pt]0\end{pmatrix}[/math]
El producto escalar [math]( \vec{e}_{\rho} \cdot \sigma \cdot \vec{e}_{\rho})[/math] se corresponde con la tensión normal anteriormente calculada, que al multiplicar por el vector[math]\vec{e_\rho}=\begin{pmatrix}1\\0\end{pmatrix}[/math] se obtiene como resultado:[math]\begin{pmatrix}\dfrac{7\rho - 4}{5} \\[6pt]0\end{pmatrix}[/math]
Se aprecia entonces que ambos productos son iguales y al restarlos para calcular la tensión tangencial se anulan. Se concluye que:
[math]\sigma \cdot \vec{e}_{p} - (\vec{e}_{p} \cdot \sigma \cdot \vec{e}_{p}) \, \vec{e}_{p} =
\begin{pmatrix}\frac{7\rho - 4}{5} \\[6pt]0\end{pmatrix}-\begin{pmatrix}\frac{7\rho - 4}{5} \\[6pt]0\end{pmatrix}=\begin{pmatrix}0 \\[6pt]0\end{pmatrix}[/math]
Que la tensión tangencial resultante sea nula y unifome para toda la superficie, significa que las deformaciones que pueda sufrir la placa serán debidas exclusivamente a las tensiones normales, que serán las únicas que aplican esfuerzos sobre el anillo. En una situación real, eta información nos permitiría aplicar medidas para reducir los esfuerzos únicamente en la componente normal, e ignorar los tangenciales puesto que en ningún punto causarán deformación.
Se representa en la siguiente gráfica la tensión tangencial. Naturalmente, el color sobre la superficie es uniforme puesto que aplica la misma tensión en toda su superficie.
%Mantenemos parte de la información de apartados anteriores
h = 0.1;
r = (1:h:2);
angulo = (0:h:pi);
[RR, AA] = meshgrid(r, angulo);
x = RR .* cos(AA);
y = RR .* sin(AA);
%Apartado 8
%Calculamos las dos tensiones normales
tensionnormal1 = (7*RR - 4)/5;
tensionnormal2 = (5*RR - 4) ./ (5 * RR.^2);
%Graficamos ambas funciones usando subplot
subplot(1,2,1);
contourf(x, y, tensionnormal1, 20);
axis equal;
xlabel('Eje X'); ylabel('Eje Y');
title('Tensión normal radial');
colorbar;grid on;
%Ahora vamos con la segunda
subplot(1,2,2);
contourf(x, y, tensionnormal2, 20);
axis equal;
10 Cálculo de las tensiones tangenciales respecto al plano ortogonal a [math]\vec{e_\theta}[/math]
Al igual que en el apartado 9 seguimos trabajando sobre la misma superficie plana:
[math]\tau_{\rho} = \left|\, \sigma \cdot \vec{e}_{θ} \;-\; \left( \vec{e}_{θ} \cdot \sigma \cdot \vec{e}_{θ} \right)\vec{e}_{θ} \,\right|[/math]
Además, el tensor de deformaciones [math]\sigma[/math] ya lo conocemos pues fue calculado en el apartado anterior:
[math]\sigma = \begin{pmatrix}\dfrac{7\rho - 4}{5} & 0 \\0 & \dfrac{5\rho - 4}{5}\end{pmatrix}[/math]
Donde en este caso sea [math] eθ =1/p*(0,1)[/math].
A partir de los cálculos sea <center>[math]\sigma*eθ=\begin{pmatrix} 0\\(5\rho - 4)5\rho^3\end{pmatrix}[/math]Posteriormente se calcula [math]eθ*\sigma*eθ=1/p*(0,1)*\sigma*1/p*(0,1)[/math]. Sustituyendo en todos los valores sea el tensor tangencial al plano ortogonal a1/p*eθ:
Su código en MATLAB y su representación sea la siguiente:
% Mantenemos parte de la información de apartados anteriores
h = 0.1;
r = (1:h:2);
angulo = (0:h:pi);
[RR, AA] = meshgrid(r, angulo);
x = RR .* cos(AA);
y = RR .* sin(AA);
% Apartado 10
% Tensión tangencial respecto al plano ortogonal a e(theta)
tensiontangencialangulo = ((5*RR - 4) .* (RR.^2 - 1)) ./ (5 * RR.^5);
% Graficamos
contourf(x, y, tensiontangencialangulo, 20);
axis equal;
xlabel('Eje X'); ylabel('Eje Y');
title('Tensión tangencial respecto al plano ortogonal a e(theta)');
colorbar;
grid on;
11 Cálculo de la masa de un campo en coordenadas polares cilíndricas
La masa se calcula como la integral de una función de densidad sobre el dominio de la superficie. En este caso, se aplica una integral doble puesto que la densidad viene definida en dos variables, [math]\rho[/math] y [math]\theta[/math].
Dado el campo escalar [math]d(\rho, \theta) = 1 + e^{\rho^{2} \cos \theta}[/math], se tiene en cuenta que los límites de integración para cada variable vendrán dados por la definición de la superficie plana, hallada en apartados anteriores. Dichos límites serán: [math]\rho=[1, 2][/math] y [math]\theta=[0,\pi][/math]
Ahora ya podemos definir la integral doble que debemos calcular para hallar la masa:
Para simplificar, dividimos la integral en dos términos de la siguiente manera:
El segundo término se calcula como una aproximación por el método del trapecio. Para ello utilizamos el siguiente código en Matlab, que calcula el área de una superficie muy próxima al valor del área bajo la curva:
% Definimos la función a integrar como anónima y los límites de integración
f = @(rho, theta) rho .* exp(rho.^2 .* cos(theta));
a = 1; b = 2;
c = 0; d = pi;
N1 = 100;
N2 = 100;
% Se subdivide el espacio creando una cuadrícula y se evalúa la función en
% cada punto
rho_vec = linspace(a, b, N1);
theta_vec = linspace(c, d, N2);
[RHO, THETA] = meshgrid(rho_vec, theta_vec);
Z = f(RHO, THETA);
%Se aplica la integral en dos partes (integral doble)
integral_rho = trapz(rho_vec, Z, 2);
result = trapz(theta_vec, integral_rho, 1);
Tras ejecutar el código, el programa nos devuelve como resultado de la aproximación el valor de 19,9337 unidades. Por tanto [math]I_2=19,9337[/math]
Por tanto:
12 Aplicaciones reales en ingeniería
Estudiar el comportamiento de ondas, sólidos y fluidos es de especial importancia en ingeniería civil para garantizar la seguridad de las estructuras. Es por ello que el uso de superficies planas como la que hemos estudiado es muy habitual para monitorizar el estado en timpo real de estructuras. Además, son datos valiosos para anticipàrse a los posibles problemas que puedan surgir y de esta forma construir obras más resilientes y seguras. Aunque el campo de aplicaciones es muy amplio, nos vamos a centrar en dos de las más útiles e interesantes:
12.1 Análisis de vibraciones
Se trata de una aplicación muy utilizada en obras subterráneas como túneles de ferrocarril, donde los trenes producen importantes vibraciones a su paso y es muy necesario estudiar su magnitud y dirección para evitar movimientos en la superficie (con la correspondiente afección a edificios) o en la propia estructura del túnel. Es un entorno especialmente delicado por su tipología subterránea, donde se debe monitorizar cualquier movimiento.
Se suelen colocar sensores que introducen placas similares a la estudiada en este trabajo. El paso de los trenes genera un vector desplazamiento y las afecciones al túnel se calculan con la divergencia y especialmente el rotacional. De esta manera, se monitoriza el túnel en distintos puntos y así se realiza un mantenimiento preventivo de la vía y se pueden establecer limitaciones de velocidad en tramos donde se superen los valores de seguridad (especialmente en curvas, donde las vibraciones aumentan por mayor rozamiento entre rueda y carril).
12.2 Diseño de puentes y voladizos
Los tableros de puentes o las estructuras voladizas deben ser muy resistentes porque se encuentran suspendidos en el aire y deben soportar tensiones de tracción y compresión garantizando su estabilidad. Por ello, es muy común utilizar sensores que miden las deformaciones producidas en el tablero de ujn purnte o voladizo cada vez que actúa una carga sobre él (por ejemplo, un vehículo). Según la magnitud y dirección del desplazamiento (magnitud vectorial) y las tensiones que aparecen en la superficie plana, es posible determinar qué zonas del tablero se encuentran más expuestas a las cargas y anticiparse a los problemas mediante refuerzos en la estructura o mantenimiento en zonas deterioradas.
12.3 Estudio de fatiga en piezas de maquinaria
En grandes máquinas existen multitud de piezas expuestas a fatiga cuyo estado debe controlarse en todo momento para evitar defectos. Es útil estudiar esas piezas como si se tratasen de placas planas y analizar cómo se distribuyen las tensiones y deformaciones por ellas. De esta manera se puede predecir la vida útil del elemennto y llevar a cabo el reemplazo antes de que muestre signos de fatiga. Se emplea frecuentemente en alas de avión, soportes de superficies o cuchillas de turbinas. centro
13 Póster del trabajo
14 Bibliografía
- Apuntes Teoría de Campos, Tema 2 (Campos tensoriales y operadores diferenciales), Curso Académico 2025/26, ETSICCP
- Google Gemini, Asistente IA
- Matlab documentation, MathWorks
