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

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.

Mallado de nuestra placa rectangular plana

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:

[math] T(x,y)=3log(1+(x-1)^2)+log(1+(y-8)^2)[/math]

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.

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:

[math]\nabla T(x,y,z) =\frac{d∂}{dx} +\frac{d∂}{dy} + \frac{d∂}{dz}[/math]

En este caso particular:

[math]\nabla T(x,y) =\frac{d∂}{dx} +\frac{d∂}{dy}=\gt3 \cdot \frac{2(x-1)}{1+(x-1)^2} \vec i + \frac{2(y-8)}{1+(y-8)^2} \vec j [/math]

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 gradiente de temperatura

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:

[math]\vec Q = -k *\nabla T[/math]
(siendo k una constante con valor 1)

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.

Ley de Fourier aplicada

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.

Campo de vectores de desplazamientos

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.

Sólido 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:

Divergencia del campo vectorial

En la imagen anterior se aprecia como el eje (x,6), es la zona de máxima divergencia. El código empleado se muestra a continuación:

%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')

7 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:

[math]\nabla×\vec u(x,y) = |\begin{matrix} \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{matrix}|[/math]

Si se aplica el campo de desplazamiento:

[math] \nabla×\vec u(ρ,θ) = -\frac{cos(2θ)(ln(ρ)+1)}{2ρ}[/math]

A continuación se muestra una representación gráfica de este resultado :

Rotacional del campo escalar

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')

8 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:

[math] \vec i · \sigma · \vec i = \sin (\frac{y·\Pi}{12}) [/math]

[math] \vec j · \sigma · \vec j = (\frac{1}{3})· \sin (\frac{y·\Pi}{12}) [/math]

[math] \vec k · \sigma · \vec k = (\frac{1}{3})· \sin (\frac{y·\Pi}{12}) [/math]

Las tensiones normales vienen representadas en la siguiente gráfica :

Tensiones normales

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]

Tensiones tangenciales

%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');

11 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:

[math]\vec{F}=\frac{∂^2\vec{u}}{∂t^2}-∇· σ[/math]

El vector [math]\vec{u}[/math] es [math]\vec{u}=\vec{a}·sin(pi*k(d·\vec{r0}(x, y)−vt))[/math] [math]\vec{a}=x/3 \vec{i} [/math] [math]\vec{d}=1/12 \vec{j}[/math] [math]k=1[/math] Sustituimos y nos da: [math]\vec{u}=\frac{x}{3}sin (\frac{\pi*y}{12}-vt)\vec{i}.[/math]