Elasticidad en Campos Escalares y Vectoriales (Coordenadas Cilíndricas)
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Elasticidad en Campos Escalares y Vectoriales (Coordenadas Cilíndricas) - Grupo 4-C |
| Asignatura | Teoría de Campos |
| Curso | 2023-24 |
| Autores | Juan Casteres Cortiñas Jorge Martínez Rodríguez-Malo Manuel Martín-Oar Faustmann Sandra Fuzhen Rodríguez Ibáñez |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
En el siguiente trabajo se expondrán los cálculos relacionados con el estudio de diferentes campos escalares y vectoriales, así como su influencia en un sólido curvado. Para ello se hará uso de coordenadas cartesianas y cilíndricas.
Se considera una placa rectangular plana (en dimensión 2) que ocupa la región [math](x, y) ∈ [−3, 3] × [−1, 3][/math]. También se suponen definidas las cantidades físicas de: la temperatura [math]T(x,y)=sin(x\cdot 2+(y−3)\cdot2)[/math] y el campo de desplazamientos: [math] u (ρ, θ) = log (3 − ρ) 2 cos(2θ)\vec{e_ρ} [/math] (Proporcionadas en coordenadas cartesianas y cilíndricas respectivamente).
Con estos datos, se procede a responder a las siguientes cuestiones:
Contenido
- 1 Dibujo Lámina
- 2 Cálculo de las curvas de nivel del campo de temperaturas
- 3 Ley de Fourier aplicada a la energía calorífica [math]\vec{Q}[/math]
- 4 Dibujo de un campo vectorial sobre el mallado, en t=0
- 5 Desplazamiento provocado por el campo vectorial [math]\vec{u}[/math]
- 6 Dibujo de [math]\nabla \cdot \vec{u}[/math] en t=0
- 7 Cálculo y dibujo de [math]|\nabla \times \vec{u}|[/math] en t=0
- 8 Tensor Tensiones
- 9 Cálculo de las tensiones tangenciales respecto al plano ortogonal a [math]\vec{e_\rho}[/math]
- 10 Tensión de Von Mises
1 Dibujo Lámina
Para el dibujo del mallado se tomará la región dada: una parte de un anillo. Como paso de muestreo para las variables [math]x[/math] e [math]y[/math], tomaremos [math]\frac{1}{10}[/math], contrariamente a los [math]\frac{2}{10}[/math] que propone el enunciado. Esto, es debido a que un paso de muestreo menor produce una malla con más puntos, y por lo tanto, con más detalle. Es por ello, que la elección de [math]\frac{1}{10}[/math] proporciona un buen balance entre detalle y claridad para todos los problemas planteados a lo largo de este artículo.
Para llevar a cabo el dibujo de la lámina (sólido), se ha empleado el programa informático MATLAB, debido a sus buenas prestaciones y sus múltiples facilidades a la hora de graficar problemas matemáticos.
Al tratarse de una geometría circular, se emplearán coordenadas cilíndricas. Sin embargo, MATLAB únicamente trabaja en coordenadas cartesianas, por lo que todos los resultados se parametrizarán para expresarlos de tal forma.
%Dibujo Lámina (sólido)
%% Datos y Variables
h=1/10; %Coordenadas r y t
r=1:h:2;
t=pi/4:h:(3*pi)/4;
%% 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);
%% Dibujo de la malla
mesh(X,Y,0*X);
axis([-3,3,-1,3]); % definimos ejes
title('Lámina');
xlabel('X');
ylabel('Y');
grid on;
view(2);
2 Cálculo de las curvas de nivel del campo de temperaturas
La temperatura del sólido viene dada por el siguiente campo escalar:
Si introducimos la expresión del campo en MATLAB, podemos representarlo gráficamente con el comando "surf", obteniendo la imagen que se presenta a continuación.
%% Cálculo función de Temperatura
T=sin(X.^2+(Y-3).^2);
%% Gráfico Mapa Térmico
figure(1);
surf(X,Y,T);
title('Temperatura');
view(2);
axis equal;
axis([-3,3,-1,3]);
xlabel('X');
ylabel('Y');
2.1 Representación curvas de nivel de T
Las curvas de nivel son líneas que unen los puntos de la misma altitud. Esto nos muestra de una forma muy visual de la variación de la función temperatura. Es posible representar estas curvas en MATLAB mediante el comando "contour".
Las zonas de la gráfica que presentan colores más azulados son aquellas en las que la temperatura es menor, mientras que en las más amarillentas la temperatura es mayor.
Las partes en las que aparecen zonas blancas son aquellas en las que la temperatura es nula.
%% Gráfico Curvas de Nivel
figure(2);
contour(X,Y,T,50);
title('Curvas de Nivel');
axis equal
axis([-3,3,-1,3]);
xlabel('X');
ylabel('Y');
colorbar
2.2 Punto máximo de T
Para la función [math]T(x,y)[/math] cabe recalcar que los valores máximos y mínimos se alcanzan en el intervalo [math]x ∈ [−1,1][/math] lo que nos da a entender que se trata de una función trigonométrica seno.
En los dos gráficos obtenidos, tanto el del campo escalar temperatura como el de las curvas de nivel de dicho campo, podemos hacer zoom y observar que en ambos casos que el punto de la gráfica en el que la temperatura es mayor es aproximadamente el (0,2).
INTRODUCIR IMAGENES AMPLIADAS
2.3 Cálculo de [math]\vec{\nabla T}[/math]
Podríamos calcular el gradiente del campo escalar T(x,y) analíticamente mediante las fórmulas vistas en clase. Sin embargo, introduciendo el campo escalar "F" y el paso "h" utilizado en apartados anteriores, el comando "gradient" nos lo calcula directamente.
El comando "contour" nos permite dibujar las curvas de nivel (apartado 2.1) y el comando "quiver" nos permite representar el campo vectorial gradiente. Por tanto, si representamos ambos a la vez mediante el comando "hold on", obtenemos el gráfico de [math]\vec{\nabla T}[/math] sobre las curvas de nivel de [math]T(x,y)[/math].
%% Gradiente de la Temperatura
[FX,FY]=gradient(T,h); %Cálculo del gradiente discretizando h
%Derivada parcial respecto de X=Fx=(2.*X)/(X.^2+1);
%Derivada parcial respecto de Y=Fy=((2.*Y)-8)/((Y-4).^2+1);
figure(4);
contour(X,Y,T,50);
title('\nablaT');
xlabel('X');
ylabel('Y');
view(2);
axis equal
axis([-3,3,-1,3]);
colorbar
hold on
quiver(X,Y,FX,FY)
axis equal
axis([-3,3,-1,3]);
view(2)
2.4 Ortogonalidad de [math]\nabla T[/math] respecto a las curvas de nivel
Representando el campo vectorial [math]\nabla T[/math] sobre las curvas de nivel, se puede apreciar que los vectores del gradiente son ortogonales a las curvas de nivel.
%% Gradiente de la Temperatura
[FX,FY]=gradient(T,h); %Cálculo del gradiente discretizando h
%Derivadas parciales
%Fx=(2.*X)/(X.^2+1);
%Fy=((2.*Y)-8)/((Y-4).^2+1);
figure(4);
contour(X,Y,T,50);
title('\nablaT');
xlabel('X');
ylabel('Y');
view(2);
axis equal
axis([-3,3,-1,3]);
colorbar
hold on
quiver(X,Y,FX,FY)
axis equal
axis([-3,3,-1,3]);
view(2)
3 Ley de Fourier aplicada a la energía calorífica [math]\vec{Q}[/math]
La energía calorífica es un término utilizado en el campo de la termodinámica, para referirse a aquella energía asociada al calor que recibe, desprende o transmite un cuerpo. Un ejemplo físico muy común de la aplicación de esta fórmula es la transmisión de calor, y por tanto, variación de temperatura a lo largo de un sólido.
Según la Ley de Fourier, la energía calorífica [math]\vec{Q}[/math] viaja de acuerdo a la fórmula [math]\vec{Q} = −\kappa\nabla T[/math], donde [math]\kappa[/math] es la constante de conductividad térmica de la placa que supondremos κ = 1.
Sabiendo que el gradiente de la temperatura representa la dirección de más velocidad de cambio de la función. El cual viene dado por:
Por tanto, queda:
Una propiedad característica del gradiente es la ortogonalidad a las curvas de nivel que hemos representado en el apartado anterior.
%% Ley de Fourier
figure(5);
title('Q=-k\nablaT');
xlabel('X');
ylabel('Y');
view(2);
hold on
%Representación campo vectorial Fourier
quiver(X,Y,-FX,-FY)
axis equal
axis([-3,3,-1,3])
view(2)
4 Dibujo de un campo vectorial sobre el mallado, en t=0
En primer lugar, tomamos el campo vectorial [math]\vec{u}(ρ, θ) = log (3 − ρ) 2 cos(2θ)\vec{e_ρ} [/math] y lo parametrizamos como se muestra en las líneas 3 y 5 del código. Esto lo hacemos para tener [math]\vec{u}[/math] expresado en coordenadas cartesianas, ya que MATLAB no puede representarlo si está expresado en cilíndricas.
Posteriormente, representamos dicho campo vectorial sobre nuestro sólido mediante la función "quiver", obteniendo el gráfico que se muestra a continuación.
%%Campo en t=0
%Componente x campo vectorial (Nueva parametrización)
ui=((log(3-U))./2).*cos(2.*V).*cos(V);
%Componente y campo vectorial
uj=((log(3-U))./2).*cos(2.*V).*sin(V);
figure(6);
title('Campo de vectores u en t=0');
xlabel('X');
ylabel('Y');
view(2);
hold on
%Representación campo vectorial u
quiver(X,Y,ui,uj)
axis equal
axis([-3,3,-1,3])
view(2)
5 Desplazamiento provocado por el campo vectorial [math]\vec{u}[/math]
Dado el campo: [math] u (ρ, θ) = log (3 − ρ) 2 cos(2θ)\vec{e_ρ} [/math] queremos representar el desplazamiento que experimenta nuestro sólido a causa de él.
5.1 Dibujo antes del desplazamiento
Previo al desplazamiento, la posición del sólido es idéntica a la representada en el apartado 1.
%Inicio
subplot(2,2,1);
j=mesh(X,Y,X*0);
set(j,'EdgeColor','g')
title('Situación antes del desplazamiento');
xlabel('X');
ylabel('Y');
axis equal;
axis([-3,3,-1,3])
view(2);
5.2 Dibujo después del desplazamiento
Tras la actuación del campo vectorial [math]\vec{u}[/math], la posición del sólido es la siguiente:
%Después del desplazamiento
subplot(2,2,2);
mesh(X+ui,Y+uj,X*0);
title('Situación después del desplazamiento');
xlabel('X');
ylabel('Y');
axis equal;
axis([-3,3,-1,3])
view(2);
5.3 Comparación pre y post deplazamiento
%Comparación entre ambas situaciones
subplot(2,2,3.5);
mesh(X+ui,Y+uj,X*0);
hold on
j=mesh(X,Y,X*0);
set(j,'EdgeColor','g');
title('Comaparación');
xlabel('X');
ylabel('Y');
axis equal;
axis([-3,3,-1,3])
view(2);
6 Dibujo de [math]\nabla \cdot \vec{u}[/math] en t=0
En coordenadas cilíndricas, el operador diferencial de la divergencia de campos vectoriales viene dado por:
Tomando nuestro campo vectorial: [math] u (\rho, \theta) = log (3 − \rho) 2 cos(2\theta)\vec{e_\rho} [/math]
(Al introducir el campo en MATLAB, parametrizamos poniéndolo en función de U y V)
%% Divergencia
figure(8)
div=((-cos(2.*V))./(2.*(3-U)))+((log(3-U).*cos(2.*V))./2.*U);
surf(X,Y,div)
axis equal
axis([-3,3,-1,3]);
xlabel('X');
ylabel('Y');
title('divergencia')
colorbar
view(2)6.1 Análisis de [math]\nabla \cdot \vec{u}[/math]
Se pide determinar de forma analítica los puntos en los que la divergencia de [math]\vec{u}[/math] es máxima, mínima y nula. Por lo que habiendo calculado el valor de [math]\nabla\cdot\vec{u}[/math] en el apartado anterior. Se procede a sacar sus derivadas:
6.2 Representación del cambio de volumen
7 Cálculo y dibujo de [math]|\nabla \times \vec{u}|[/math] en t=0
El rotacional de un campo vectorial muestra la tendencia de un campo a inducir rotación alrededor de un punto. En coordenadas cilíndricas, el operador diferencial del rotacional de campos vectoriales viene dado por:
Sustituyendo para los valores de nuestro campo y operando, obtenemos:
Tras calcular [math]\nabla \times \vec{u}[/math], sabemos que es un campo vectorial. Sin embargo, lo pedido es [math]|\nabla \times \vec{u}|[/math], que es un campo escalar. Por tanto, se hace se procede a hacer el módulo del rotacional calculado anteriormente.
Al introducir el campo en MATLAB, parametrizamos poniéndolo en función de U y V:
Si representamos gráficamente el campo escalar calculado, obtenemos las figuras que se muestran debajo.
figure(9)
ROT=sqrt(((1./U).*(-sin(2*V)).*log(3.-U)).^2); %Módulo del campo escalar.
%%ROT=abs((1./U).*(-sin(2*V)).*log(3.-U)); %Valor absoluto del campo.
subplot(1,2,1)
%Rotacional en 2D
pcolor(X,Y,ROT);
colorbar;
title('Rotacional en 2D');
xlabel('Eje X');
ylabel('Eje Y');
axis on;
axis([-3,3,-1,3]);
hold on
plot(0.7,0.7,'ro') %Máximo 1 en gráfico 2D
hold on
plot(-0.65,0.75,'ro') %Máximo 2 en gráfico 2D
%Rotacional en 3D
subplot(1,2,2)
surf(X,Y,ROT);
view(3)
colorbar;
title('Rotacional en 3D');
xlabel('Eje X');
ylabel('Eje Y');
zlabel('Eje Z');
axis on
axis([-3,3,-1,3]);
hold on
plot3(0.7,0.7,0.693,'ro') %Máximo 1 en gráfico 3D
hold on
plot3(-0.65,0.75,0.693,'ro') %Máximo 2 en gráfico 3D
7.1 Cálculo punto máximo del rotacional
Como puntos máximos de rotacional, como se puede observar en el gráfico, nos quedan los siguientes puntos:
Los cuales están señalados en ambos gráficos con un círculo rojo.
8 Tensor Tensiones
En un medio elástico lineal, isótropo y homogéneo se puede describir el tensor de tensiones mediante la expresión:
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].
Partiendo de la expresión matricial del gradiente:
Podemos calcular fácilmente [math]\epsilon(\vec{u})[/math]. Además, se tiene que la matriz identidad multiplicada por la divergencia del campo [math]\vec{u}[/math] da como resultado la matriz:
Por lo tanto, operando con las matrices anteriores se obtiene la matriz del tensor tensiones:
8.1 Tensiones normales en la dirección de [math]\vec{e_\rho}[/math]
Vienen definidas por [math]\vec{e_\rho}·\sigma·\vec{e_\rho}[/math]. Por tanto, multiplicando la matriz [math]\sigma[/math] a ambos lados por el vector unitario [math]\vec{i}[/math], obtenemos las tensiones normales en la dirección de [math]\vec{i}[/math].
8.2 Tensiones normales en la dirección de [math]\vec{e_\theta}[/math]
Vienen definidas por [math]\vec{e_\theta}·\sigma·\vec{e_\theta}[/math]. Por tanto, multiplicando la matriz [math]\sigma[/math] a ambos lados por el vector unitario [math]\vec{j}[/math], obtenemos las tensiones normales en la dirección de [math]\vec{j}[/math].
8.3 Tensiones normales en la dirección de [math]\vec{e_z}[/math]
Vienen definidas por [math]\vec{e_z}·\sigma·\vec{e_z}[/math]. Por tanto, multiplicando la matriz [math]\sigma[/math] a ambos lados por el vector unitario [math]\vec{k}[/math], obtenemos las tensiones normales en la dirección de [math]\vec{k}[/math].
8.4 Código Matlab
%% Apartado 8
i_sig_i=(-3./2).*(cos(2.*V)./(3-U))+(log(3-U).*cos(2.*V))./(2.*U);
j_sig_j=(-cos(2.*V)./(2.*(3-U)))+(3./2).*((log(3-U).*cos(2.*V))./U);
k_sig_k=(-cos(2.*V)./(2.*(3-U)))+((log(3-U).*cos(2.*V))./(2.*U));
figure(10)
subplot(2,2,1)
surf(X,Y,i_sig_i)
view(2)
axis equal
axis([-3,3,-1,3]);
colorbar
title('i·\sigma·i')
subplot(2,2,2)
surf(X,Y,j_sig_j)
view(2)
axis equal
axis([-3,3,-1,3]);
colorbar
title('j·\sigma·j')
subplot(2,2,3.5)
surf(X,Y,k_sig_k)
view(2)
axis equal
axis([-3,3,-1,3]);
colorbar
title('k·\sigma·k')
Los tres campos escalares resultantes representan las tensiones normales en la dirección que marca el eje [math]\vec{e_\rho}[/math],las tensiones normales en la dirección que marca el eje [math]\vec{e_\theta}[/math], y las correspondientes al eje [math]\vec{e_z}[/math].
9 Cálculo de las tensiones tangenciales respecto al plano ortogonal a [math]\vec{e_\rho}[/math]
Dichas tensiones vienen definidas por la expresión:
Puesto que [math](\vec{e_\rho}\cdot\sigma\cdot\vec{e_\rho})[/math] lo tenemos calculado del apartado 8.2, bastará con multiplicar el valor hallado por [math]\vec{e_\rho}[/math].
[math]\sigma\cdot\vec{e_\rho}[/math] se calcula de la siguiente forma:
Por lo tanto, el vector solución resultante de la resta de los anteriores es:
Dado que en Matlab se trabaja en coordenadas cartesianas, Es necesario pasar de la base física cilíndrica a la cartesiana.
Por lo tanto el módulo queda tal que:
tti=((log(3-U).*sin(2.*V).*sin(V)./U));
ttj=((-log(3-U).*sin(2.*V).*cos(V)./U));
mod=sqrt((tti.^2)+(ttj.^2));
figure(11)
hold on
mesh(X,Y,X.*0)
surf(X,Y,mod)
axis equal
axis([-3,3,-1,3]);
title('módulo de tensiones tangenciales respecto al plano ortogonal a i');
10 Tensión de Von Mises
La tensión de Von Mises viene definida por la fórmula:
donde [math]\sigma_1 , \sigma_2 , \sigma_3 [/math] son los autovalores de σ (también conocidos como tensiones principales). Se trata de una magnitud escalar que se suele usar como indicador para saber cuando un material inicia un comportamiento plástico (y no elástico puro).
Para sacar los autovalores en MATLAB se utiliza el comando eig.m
%% Matriz Tensión de Von Misses
sig=[];
M_VonMises=zeros(length(t),length(r));
a=@(UU,VV)(-3./2).*(cos(2.*VV)./(3-UU))+(log(3-UU).*cos(2.*VV))./(2.*UU);
b=@(UU,VV)(-cos(2.*VV)./(2.*(3-UU)))+(3./2).*((log(3-UU).*cos(2.*VV))./UU);
c=@(UU,VV)(-cos(2.*VV)./(2.*(3-UU)))+((log(3-UU).*cos(2.*VV))./(2.*UU));
d=@(UU,VV)((-log(3-UU).*sin(2.*VV))./UU);
%fila1=[a,d,0];
%fila2=[d,b,0];
%fila3=[0,0,c];
%sigma=[fila1;fila2;fila3]
% Se ha de asignar a la matriz "M_VonMises" los valores
% de la tensión de Von Mises en cada punto de la placa.
for i=1:length(t)
for j=1:length(r)
sig(1,1)=a(U(i,j),V(i,j));
sig(1,2)=d(U(i,j),V(i,j));
sig(1,3)=0;
sig(2,1)=d(U(i,j),V(i,j));
sig(2,2)=b(U(i,j),V(i,j));
sig(2,3)=0;
sig(3,1)=0;
sig(3,2)=0;
sig(3,3)=c(U(i,j),V(i,j));
[v,u]=eig(sig);%Autovalores
M_VonMises(i,j)=sqrt(((u(1,1)-u(2,2)).^2+((u(2,2)-u(3,3)).^2)+((u(3,3)-u(1,1)).^2))/2);
end
end
% Von Mises.10.1 Representación de la Tensión de Von Mises
%Representación Von Mises
figure(12)
surf(X,Y,M_VonMises);
% Se pone título a la gráfica.
title('Von Mises');
% Se da nombre a los ejes.
xlabel('Eje de las X');
ylabel('Eje de las Y');
axis equal
axis([-3,3,-1,3]);
% Se da equidistancia a los ejes.
% Se le aplica una barra con colores en función
% de la tension de Von Mises en cada punto.
colorbar;
view(2)
figure(13)
surf(X,Y,M_VonMises);
% Se pone título a la gráfica.
title('Von Mises');
% Se da nombre a los ejes.
xlabel('Eje de las X');
ylabel('Eje de las Y');
axis equal
axis([-3,3,-1,3]);
% Se da equidistancia a los ejes.
% Se le aplica una barra con colores en función
% de la tension de Von Mises en cada punto.
colorbar;
view(3)
