Visualización de campos escalares y vectoriales. Grupo 29.
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Visualización de campos escalares y vectoriales en elasticidad. Grupo 29 |
| Asignatura | Teoría de Campos |
| Curso | 2023-24 |
| Autores |
Oliver Prada Sanchidrián Rafael Garcia Lopez Gonzalo Ramirez Mateo Alejandra Martin Moreno Carlos de Ana de Miguel |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
Contenido
1 Introducción
El propósito de este artículo es analizar el desplazamiento experimentado por una pieza después de aplicar una fuerza.
Para ello, se tienen dos cantidades físicas dependientes de las variables x e y:
- La temperatura T(x,y)
- El campo de desplazamientos [math]\vec u(x,y)[/math], producido por la acción de una fuerza.
Para ello, se utilizará el software de programación y cálculo numérico Matlab/Octave.
2 Presentación de la placa
Para simplificar las operaciones, se va a tomar una sección del objeto a estudiar, esta es la resultante de intersecar a la pieza con un plano ortogonal a su eje axial, lo que resulta una placa rectangular plana, centrada en el origen. A continuación se muestra la sección de la placa.
El código empleado para dibujar la placa es el siguiente:
%MALLADO
%Región de la placa
h=2/10
x=[-1:h:1];
y=[0:h:12];
%Mallado con las matrices Mx e My
[Mx,My]=meshgrid(x,y);
mesh(Mx,My,0*My);
%Ejes
axis([-5,5,-0.5,12.5])
%Region del dibujo
title(‘Mallado de la placa’);
xlabel(‘Eje X’);
ylabel(‘Eje y’);
view(2);
3 Temperatura
La temperatura viene determinada por la siguiente función:
3.1 Curvas de nivel
Esta fórmula nos proporciona la siguiente representación gráfica. La primera imagen muestra por colores la variación en el recinto observable, y la segunda representa las curvas de nivel del campo de temperaturas.
En las gráficas anteriores se aprecia, como la zona mas fría está concentrada en la parte superior derecha. Según se aleja del foco de frío la placa va aumentando la temperatura, el punto en el que alcanza su valor máximo es el (-1,0).
El código empleado para dibujar las gráficas anteriores es el siguiente:
%CURVAS DE NIVEL Y TEMPERATURA DE LA PLACA
x=-1:0.2:1
y=0:0.2:12
[X,Y]=meshgrid(x,y)
%Definición de la temperatura
T=3*log(1+(X-1).^2)+log(1+(Y-8).^2);
%Representación de curvas de nivel
hold on
subplot(1,2,1)
surf(X,Y,T)
view(2)
title(“Temperatura”)
xlabel(“Eje X”)
ylabel(“Eje Y”)
colorbar
subplot(1,2,2)
contour(X,Y,T,0)
colorbar
axis([-1,1,0,12])
title(“Curvas de nivel de la temperatura”)
xlabel(“Eje X”)
ylabel(“Eje Y”)
hold off
3.2 Gradiente de la temperatura
Para realizar el cálculo del gradiente, se emplea su respectiva expresión en cartesianas:
En este caso particular:
A continuación se presenta una representación del gradiente de la temperatura, en el que se puede apreciar como el campo de vectores es ortogonal a las curvas de nivel de la temperatura, ya que marca la dirección de máximo crecimiento de la misma.
El código empleado para dibujar las gráficas anteriores es el siguiente:
%Gradiente de Temperatura
h=0.2;
u=-1:h:1;
v=0:h:12;
[Mx,My]=meshgrid(u,v);
Temp=3*log(1 +(X-1).^2)+log(1+(Y-8).^2);
contour(Mx,My,Temp,30);
%Añadimos al dibujo de las líneas de nivel del gradiente de temperatura
hold on
[Px,Py]=gradient (Temp,0.1,0.1);
quiver(u,v, Px, Py)
axis([-1,1,0,12]);
colorbar
hold off
xlabel('Eje X')
ylabel('Eje Y')
title ('Gradiente de la Temperatura')
4 Ley de Fourier
La ley de Fourier, conocida también como la ley calorífica [math]\vec Q[/math],
afirma que existe una proporcionalidad entre el flujo de la energía y el gradiente de la temperatura, viaja de de acuerdo la fórmula:
Esta ley aplicada a nuestro caso particular, proporciona un campo vectorial. Este campo se asemeja al obtenido en el apartado anterior, pero con la dirección opuesta.
Esta gráfica adjunta corresponde al siguiente código de programación:
%LEY DE FOURIER
h=0.2;
u=-1:h:1;
v=0:h:12;
[Mx,My]=meshgrid(u,v);
Temp=3*log(1 +(X-1).^2)+log(1+(Y-8).^2);
contour(Mx,My,Temp,30);
%Añadimos al dibujo de las líneas de nivel del gradiente de temperatura
hold on
[Px,Py]=gradient (Temp,0.1,0.1);
quiver(u,v, -1.*Px, -1.*Py)
axis([-1,1,0,12]);
colorbar
hold off
xlabel('Eje X')
ylabel('Eje Y')
title ('LEY DE FOURIER')
5 Desplazamiento
5.1 Campo de vectores de desplazamientos
El campo de desplazamiento corresponde al campo de vectores generado en los puntos de mallado del sólido, cuando el tiempo es nulo, es decir, t=0.
El código correspondiente es:
%CAMPO DE DESPLAZAMIENTOS
h=0.2;
u=-1:h:1;
v=0:h:12;
[Mx,My]=meshgrid(u,v);
%Definimos vector desplazamiento
Ux= (Mx/3).*sin((pi*My)/12);
Uy= zeros(size(Mx));
%Representamos campo vectorial
quiver(Mx,My,Ux,Uy);
axis([-1,1,0,12]);
xlabel('Eje X')
ylabel('Eje Y')
view(2)
title ('Campo de desplazamientos')
5.2 Sólido antes y después del desplazamiento
El campo de deformaciones mencionado anteriormente, podemos asumir que se generarán ondas a lo largo de la placa. Por ello, para observar la variación, será necesario la realización de la placa antes y después del desplazamiento.
El código empleado para llevar a cabo la comparación se expone a continuación:
%SOLIDO ANTES Y DESPUÉS DEL DESPLAZAMIENTO
h=0.2;
u=-1:h:1;
v=0:h:12;
[Mx,My]=meshgrid(u,v);
Ux= (Mx/3).*sin((pi*My)/12);
Uy= zeros(size(Mx));
figure(1)
%ANTES
subplot(2,2,1)
mesh(Mx,My,0*Mx);
axis([-1,1,0,12]);
title('Antes del desplazamiento');
view(2)
xlabel('Eje X')
ylabel('Eje Y')
%DESPUÉS
subplot(2,2,2)
mesh(Ux+Mx ,Uy+My,0*Ux);
axis([-1,1,0,12]);
title('Después del desplazamiento')
view(2)
xlabel('Eje X')
ylabel('Eje Y')
%COMPARACIÓN
subplot(2,2,3)
plot3(Mx,My,Mx*0);
hold on
plot3(Ux+Mx,Uy+My,0*Ux);
hold off
view(2)
axis([-1,1,0,12]);
title ('Comparación')
hold off
xlabel('Eje X')
ylabel('Eje Y')
6 Divergencia
La divergencia de un campo vectorial es una medida del cambio de volumen local debido al desplazamiento.
Se deben considerar dos expresiones a la hora del cálculo:
- El campo de de desplazamientos: [math]\vec u (x,y,t) = -\vec a sin(πk( \vec d\cdot\vec r_0 (x,y) - vt)) = \frac{x}{3} ({sen(\frac{πy}{12}}) \vec i [/math]
- La fórmula de la divergencia: [math]\nabla\cdot\vec u= \frac{\partial}{\partial{x}}({\frac{x}{3} ({sen(\frac{πy}{12}})})[/math], es decir, [math]\nabla\cdot\vec u= \frac{1}{3} ({sen(\frac{πy}{12}}) [/math]
La representación gráfica de la divergencia del campo de desplazamiento resulta:
La Divergencia de [math] \vec u \lt\math\gt es máxima en el eje (x,6) y su valor es igual a 0,333. El mínimo punto de Divergencia de \ltmath\gt \vec u \lt\math\gt esta situado en los extremos de la placa, y su valor corresponde a 0
[[Archivo:divergencia2.15.jpeg|600px|thumb|right|2. Divergencia del campo vectorial]]
{{matlab|codigo=
%DIVERGENCIA
h=0.2;
u=-1:h:1;
v=0:h:12;
[Mx,My]=meshgrid(u,v);
figure(2)
%Representamos la divergencia
div= (1/3).*sin((pi*My)/12);
surf(Mx,My, div);
view(2)
xlabel('Eje X')
ylabel('Eje Y')
colorbar
title ('Divergencia')
%Divergencia 3D
h=2/10;
x=-1:h:1;
y=0:h:12;
[Mx,My]=meshgrid(x,y);
%Introducimos la divergencia
div=(1/3).*sin((pi/12).*My);
%Representación
surf(Mx,My,div)
shading interp
colorbar
title('Divergencia')
xlabel('x')
ylabel('y')
zlabel('z')
%Dibujo de la divergencia
figure
pcolor(Mx,My,div)
shading interp
hold on
contour(Mx,My,div,'k')
title('Divergencia')
xlabel('x')
ylabel('y')
colorbar
axis([-1,1,0,12])
view(3)
hold off
maximo=max(max(div))
minimo=min(min(div))
}}
==Rotacional==
El rotacional de un campo vectorial muestra la tendencia de un campo de un camp a inducir rotación alrededor de un punto. Viene determinad por la siguiente expresión:
\ltcenter\gt\ltmath\gt\nabla×\vec u(x,y) = {\begin{vmatrix}
\vec i & \vec j & \vec k \\
\frac{\partial}{\partial x} & \frac{\partial}{\partial y} & \frac{\partial}{\partial z} \\
\vec u_x & \vec u_y & \vec u_z
\end{vmatrix}} [/math]</center>
Si se aplica el campo de desplazamiento:
A continuación se muestra una representación gráfica de este resultado :
Como puede observarse, las zonas en las que el rotacional adopta su valor máximo se encuentran en los polos de la placa. Es decir, en el (-1 , 0) y en el (1 , 12).
El código empleado para dibujar el rotacional del campo de desplazamientos es el siguiente:
%Rotacional
h = 0.2;
u = -1:h:1;
v = 0:h:12;
[Mx, My] = meshgrid(u, v);
% Definimos rotacional
Rot = (-Mx .* pi .* cos((My .* pi) / 12)) / 36; % Corrección en la fórmula
% Representamos rotacional
surf(Mx, My, Rot);
axis([-1, 1, 0, 12]);
view(2)
colorbar
xlabel('Eje X')
ylabel('Eje Y')
title('Rotacional')
7 Tensor de tensiones
El tensor de tensiones viene determinado por la siguiente expresión: [math]σ=λ·\nabla·\vec u I+2με[/math]
Desarrollando:
[math]\sigma = \left(\begin{matrix} (\frac{1}{3}sin(\frac{πy}{12}) & 0 & 0 \\ & \\0 & (\frac{1}{3}sin(\frac{πy}{12}) & 0 \\ & \\0 & 0 & (\frac{1}{3}sin(\frac{πy}{12})\end{matrix} \right) + \left (\begin{matrix} (\frac{2}{3}sin(\frac{πy}{12}) & (\frac{πx}{36})cos(\frac{πy}{12}) & 0 \\&\\ (\frac{πx}{36})cos(\frac{πy}{12}) & 0 & 0 \\ 0 & 0 & 0 \end{matrix} \right) [/math]
Finalmente resulta: [math]σ[/math] = [math] \left (\begin{matrix} (sin(\frac{πy}{12}) & (\frac{πx}{36})cos(\frac{πy}{12}) & 0 \\&\\ (\frac{πx}{36})cos(\frac{πy}{12}) & (\frac{1}{3})sin(\frac{πy}{12}) & 0 \\&\\ 0 & 0 & (\frac{1}{3})sin(\frac{πy}{12})\end{matrix} \right)
[/math]
8.1 Tensiones normales
Los módulos de las tensiones normales aplicadas según las direcciones dadas son:
Las tensiones normales vienen representadas en la siguiente gráfica :
La imagen ha sido obtenida con el siguiente código:
%TENSIONES NORMALES
h = 0.2;
u = -1:h:1;
v = 0:h:12;
[Mx, My] = meshgrid(u, v);
% Coeficientes de Lamé
lambda = 1;
mu = 1;
% Tensor de deformaciones
epsilon = zeros(size(Mx, 1), size(Mx, 2), 2, 2);
epsilon(:, :, 1, 1) = gradient(Mx) / h;
epsilon(:, :, 2, 2) = gradient(My) / h;
epsilon(:, :, 1, 2) = 0.5 * (gradient(My) + gradient(Mx)) / h;
epsilon(:, :, 2, 1) = epsilon(:, :, 1, 2);
% Tensor de tensiones
sigma = zeros(size(Mx, 1), size(Mx, 2), 2, 2);
sigma(:, :, 1, 1) = lambda * sum(epsilon(:, :, 1, 1), 3) + 2 * mu * epsilon(:, :, 1, 1);
sigma(:, :, 2, 2) = lambda * sum(epsilon(:, :, 2, 2), 3) + 2 * mu * epsilon(:, :, 2, 2);
sigma(:, :, 1, 2) = 2 * mu * epsilon(:, :, 1, 2);
sigma(:, :, 2, 1) = sigma(:, :, 1, 2);
% Tensiones normales en la dirección que marca el eje i
subplot(2, 2, 1)
quiver(Mx, My, sigma(:, :, 1, 1), sigma(:, :, 2, 1));
xlabel('Eje X')
ylabel('Eje Y')
axis([-1, 1, 0, 12]);
title('Tensiones normales a la dirección i')
colorbar
% Tensiones normales en la dirección que marca el eje j
subplot(2, 2, 2)
quiver(Mx, My, sigma(:, :, 1, 2), sigma(:, :, 2, 2));
xlabel('Eje X')
ylabel('Eje Y')
axis([-1, 1, 0, 12]);
title('Tensiones normales a la dirección j')
colorbar
% Tensiones normales en la dirección que marca el eje k
subplot(2, 2, 3)
quiver(Mx, My, zeros(size(Mx)), zeros(size(Mx))); % No hay componente en la dirección k
xlabel('Eje X')
ylabel('Eje Y')
axis([-1, 1, 0, 12]);
title('Tensiones normales a la dirección k')
colorbar
8.2 Tensiones tangenciales
Para calcular las tensiones tangenciales que sufre la placa en la dirección del plano ortogonal a [math]\vec{i}[/math], es decir [math]|σ ·\vec{i} − (\vec{i} · σ ·\vec{i})\vec{i}|[/math], en t = 0.
[math] \left| \sigma \cdot \vec{i} - (\vec{i} \cdot \sigma \cdot \vec{i}) \vec{i} \right|=\left| (\sigma) \cdot \begin{pmatrix}1\\0 \\ 0\end{pmatrix} - sen(\frac{\Pi y}{12})\cdot \begin{pmatrix}1\\0 \\ 0\end{pmatrix} \right |= \left| ( \frac{πx}{36}cos(\frac{πy}{12}) \vec{j} ) \right| [/math]
%Tensores tangenciales
h= 2/10;
u=-1:h:1;
v= 0:h:12;
%Creación de matriz x e y
[Mx,My]=meshgrid(u,v);
%Tensión tangencial en cada punto
tn=(Mx*pi/36).*cos((pi/12).*My);
%Representación gráfica
quiver(Mx,My,tn,tn.*0);
axis([-1,1,0,12]);
title('Tension tangencial');
xlabel('X');
ylabel('Y');8 La tensión de Von Mises
9 El campo de fuerzas [math]\vec F[/math]
El campo de fuerzas que actúa sobre la placa y causa el desplazamiento se aproxima usando:
∇·σ es el campo vectorial que se obtiene al hacer la divergencia de los vectores cuyas componentes son las filas de la matriz σ (previamente calculada)
el tensor deformaciones es finalmente igual a:
Teniendo en cuenta a matriz anterior, el tensor de tensiones resulta:
Con nuestro vector [math]\vec{u}=\vec{a}·sin(pi\cdot k(d·\vec{r0}(x, y)−vt))[/math] y los desplazamientos:
[math]\vec{a}=x/3 \vec{i} [/math];
[math]\vec{d}=1/12 \vec{j}[/math];
[math]k=1[/math];
Sustituimos y nos da el vector [math]\vec u[/math]:
[math]\vec{u}=\frac{x}{3}sin (\frac{\pi \cdot y}{12}-vt)\vec{i}.[/math]
Debemos derivarlo dos veces en función de t:
[math]\frac{∂\vec{u}}{∂t}= vx/3 ·cos (\frac{\pi \cdot y}{12}-vt)\vec{i}.[/math]
[math]\frac{∂^2\vec{u}}{∂t^2}= -v^2x/3 (sin(\frac{\pi \cdot y}{12}-vt))\vec{i}.[/math]
Sustituyendo los dos términos en la ecuación de la elasticidad lineal de [math]\vec F[/math] y suponiendo que [math]\vec F[/math]=0;
[math]\vec{F}=\frac{∂^2\vec{u}}{∂t^2}-∇· σ = \frac{-v^2x}{3} (sin(\frac{Πy}{12}-vt))\vec{i}. - \begin{pmatrix} \ \frac{-xΠ^2}{432}sin(\frac{Πy}{12}-vt)\vec{i} & \frac{Π}{18}cos(\frac{Πy}{12})\vec{j} & 0\vec{k} \end{pmatrix} =0 [/math]
Para que se pueda cumplir
[math]\frac{Π}{18}cos(\frac{Πy}{12})\vec{j}.[/math] =0
Y para ello [math]cos(\frac{Πy}{12}) [/math]= 0
Despejando se llega a la igualdad de [math]V= \frac{Π}{12}[/math]
Si la onda fuera longitudinal, es decir, tomando [math]\vec{a}=1/3 \vec{i} [/math]
Nuestro nuevo vector sustituído [math]\vec u[/math] nos da:
[math]\vec{u}=\frac{1}{3}sin (\frac{\pi \cdot y}{12}-vt)\vec{i}.[/math]
Debemos derivarlo dos veces en función de t:
[math]\frac{∂\vec{u}}{∂t}= v/3 ·cos (\frac{\pi \cdot y}{12}-vt)\vec{i}.[/math]
[math]\frac{∂^2\vec{u}}{∂t^2}= -v^2/3 (sin(\frac{\pi \cdot y}{12}-vt))\vec{i}.[/math]
El nuevo tensor de tensiones es igual a [math]σ=λ∇·\vec{u}1 + 2μ∈[/math] [math] \mathbb = \; \begin{pmatrix} {sin(\frac{Πy}{12}-vt)} & {\frac{Π}{36}cos(\frac{Πy}{12}-vt)} & {0}\\ {\frac{Π}{36}cos(\frac{Πy}{12}-vt)} & {\frac{1}{3}sin(\frac{Πy}{12}-vt)} & {0}\\ {0} & {0} & {\frac{1}{3}sin(\frac{Πy}{12}-vt)}\\ \end{pmatrix} [/math]
Siendo [math] ∇· σ= \begin{pmatrix} \ \frac{-Π^2}{432}sin(\frac{Πy}{12}-vt)\vec{i} & \frac{Π}{18}*cos(\frac{Π*y}{12})\vec{j} & 0\vec{k} \end{pmatrix} [/math]
Sustituyendo los dos términos en la ecuación de la elasticidad lineal de [math]\vec F[/math] y suponiendo que [math]\vec F[/math]=0;
[math]\vec{F}=\frac{∂^2\vec{u}}{∂t^2}-∇· σ = \frac{-v^2}{3} (sin(\frac{Πy}{12}-vt))\vec{i}. - \begin{pmatrix} \ \frac{-Π^2}{432}sin(\frac{Πy}{12}-vt)\vec{i} & \frac{Π}{18}cos(\frac{Πy}{12})\vec{j} & 0\vec{k} \end{pmatrix} =0 [/math]
Despejando se llega a la misma solución que antes de [math]V= \frac{Π}{12}[/math]
