Visualización de campos escalares y vectoriales en elasticidad. (Grupo 25C)

Trabajo realizado por estudiantes
Título	Visualización de campos vectoriales y escalares en elasticidad. Grupo 25C
Asignatura	Teoría de Campos
Curso	2014-15
Autores	María Moreno, Pablo Goenechea, Andrés Vila, Victor Molina, Jose Ramos
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Presentación

La siguiente práctica tiene como objetivo el estudio de los campos sobre un sólido elástico tras ser sometido a tensión, en función de los desplazamientos y de la temperatura de cada punto del sólido.

Un campo físico representa la distribución de una magnitud física que varía en el espacio, en otras palabras, es una aplicación que asocia a cada punto del espacio un tensor. Hablamos de campos escalares y vectoriales.

El sólido que vamos a estudiar es una placa rectangular que ocupa la región plana producto cartesiano de los intervalos [math][-1/2,1/2]\times[0,4][/math].

Previo al estudio, debemos definir la placa en nuestro programa informático. Para ello se configura un mallado sobre dicha placa que determine los puntos de análisis, de forma que cada intervalo es dividido en subintervalos de valor constante [math]H = \frac{1}{10}[/math].

2 Representación de la placa y mallado

% Representación tridimensional y vista ortogonal de la placa

x = -0.5:0.1:0.5;      % Intervalo de puntos en el eje "x"
y = 0:0.1:4;           % Intervalo de puntos en el eje "y"
[X,Y] = meshgrid(x,y); % Creación de las matrices que representan 
                       % los puntos de la malla
subplot(1,2,1);
mesh(X,Y,0*X);         % Representación espacial
title('Representación espacial');
xlabel('Eje "x"');
ylabel('Eje "y"');
zlabel('Eje "z"');
subplot(1,2,2);
mesh(X,Y,0*X); 
axis([-1.5 1.5 -1 5]);
title('Vista Ortogonal');
xlabel('Eje "x"');
ylabel('Eje "y"');
view(2);               % Vista ortogonal a la placa

La malla nos proporciona una idea de la distribución espacial del sólido

3 Campo escalar Temperatura

Considerando que sobre el sólido actúa un campo escalar que representa la temperatura, independiente del tiempo y que está definido para cada punto de la malla por la función;

[math]T(x,y) = (8-y^2+2y)e^{-(x+1)^2}[/math]

Vamos a representar el campo.

% Campo escalar Temperatura de los puntos de la malla

x = -0.5:0.1:0.5;      % Intervalo de puntos en el eje "x"
y = 0:0.1:4;           % Intervalo de puntos en el eje "y"
[X,Y] = meshgrid(x,y); % Creación de las matrices que representan 
                       % los puntos de la malla

ft = inline('(8-y.^2+2*y).*exp(-(x+1).^2)','x','y'); % Se crea la función 
% temperatura por medio del comando "inline"
T = ft(X,Y); % Cálculo del campo de temperaturas

subplot(1,2,1);
mesh(X,Y,0*X);
hold on;
surfc(X,Y,T);
title('Función de temperatura');
xlabel('Eje "x"');
ylabel('Eje "y"');
zlabel('Eje "z"');
colorbar;

subplot(1,2,2);
[C,h] = contour(X,Y,T,20);
set(h,'ShowText','on','TextStep',get(h,'LevelStep')*2);
title('Vista ortogonal de las curvas de nivel');
xlabel('Eje "x"');
ylabel('Eje "y"');
colormap hot;
colorbar;

La representación gráfica de las curvas de nivel, junto con sus respectivos valores, de la función Temperatura, permite observar las zonas de crecimiento, pudiendo aproximar a simple vista la zona más caliente de la placa. Una buena aproximación del punto de máxima temperatura puede ser (x , y) = (-0.5 , 1) dónde la función vale T(-0.5 , 1)= 7.01

Sin embargo para conocer el valor exacto recurrimos al cálculo vectorial del campo gradiente de un campo escalar (siendo este diferenciable), en nuestro caso la temperatura, lo que nos permite conocer la dirección de máximo crecimiento de este campo escalar, desde cualquier punto del mismo. Las diferentes direcciones deben ser ortogonales a las curvas de nivel del campo escalar.

4 Campo vectorial Gradiente del campo Temperatura

El gradiente es un operador que, aplicado a cualquier campo escalar, da como resultado un campo vectorial, de vectores, cuyas direcciones y sentido son las de máximo crecimiento de la función que es campo escalar, en nuestro caso la temperatura. Luego la expresión del gradiente para el campo escalar temperatura es:

[math]∇.T= \frac{∂T}{∂x} i ⃗+ \frac{∂T}{∂y}j ⃗+ \frac{∂T}{∂z} \vec k[/math]

que en nuestro caso, para una determinada función de temperatura, hemos obtenido:

[math]∇T= \frac{∂T}{∂x} i ⃗+ \frac{∂T}{∂y}j ⃗+ \frac{∂T}{∂z} \vec k = (8-y^2+2y).e^{-(x+1)^2}.(-2(x+1)).(\vec î) + e^{-(x+1)^2}.(-2y+2).(\vec ĵ)[/math]

De esta manera podemos representar el campo vectorial por medio del siguiente programa;

%Cálculo y representación del campo vectorial gradiente de temperatura

x = -0.5:0.1:0.5;      % Intervalo de puntos en el eje "x"
y = 0:0.1:4;           % Intervalo de puntos en el eje "y"
[X,Y] = meshgrid(x,y); % Creación de las matrices que representan 
                       % los puntos de la malla
                       
ft = inline('(8-y.^2+2*y).*exp(-(x+1).^2)','x','y');
T = ft(X,Y);

% Cálculo del campo gradiente
gx = inline('(8-y.^2+2*y).*exp(-(x+1).^2).*(-2*(x+1))','x','y'); % Función que define las coordenadas en x del gradiente
gy = inline('exp(-(x+1).^2).*(-2*y+2)','x','y'); % Función de las coordenadas en y del gradiente

Gx = gx(X,Y); % Cálculo de las coordenadas en x del gradiente para los valores de la malla
Gy = gy(X,Y); % Coordenadas en y

% Representación del gradiente y curvas de nivel de la temperatura sobre la
% placa

quiver3(X,Y,0*X,Gx,Gy,0*Gx);
hold on;
contour(X,Y,T,10);
hold on;
mesh(X,Y,0*X);
view(2);
title('Campo Gradiente de temperatura');
xlabel('Eje "x"');
ylabel('Eje "y"');

Como podemos observar los vectores que conforman el campo son ortogonales a las curvas de nivel de la función Temperatura. Este hecho se debe a que la manera más rápida de recorrer la superficie es en la dirección de máxima pendiente, que corresponde a la mínima distancia entre curvas de nivel, es decir, perpendicular a éstas, y que está definida por el Gradiente.

5 Campo vectorial Desplazamientos y Representación de la placa antes y después de la transformación

Además de la temperatura, también se ha definido una deformación sobre la placa, debida a la acción de una fuerza, que podemos representar como un campo vectorial al que llamamos Campo de Desplazamientos [math]U(x,y)[/math].

De esta forma si definimos [math]R_0(x,y)[/math] el vector de posición de los puntos de la placa antes de la deformación, la posición de cada punto [math](x,y)[/math] de la placa después de la deformación viene dada por [math]R(x,y) = R_0(x,y) + U(x,y)[/math].

Vamos a suponer que la fuerza [math]F[/math] aplicada sobre la placa, ha provocado un desplazamiento de los puntos de la misma dado por el vector de desplazamientos [math]U(x,y) = \vec a.(\vec b.R_0)^2[/math] Dónde [math]\vec a[/math] y [math]\vec b[/math] son dos vectores dados tal que [math]\vec a = î/20 , \vec b = ĵ[/math].

% Campo vectorial de desplazamientos y forma de la placa tras la
% deformación

x = -0.5:0.1:0.5;      % Intervalo de puntos en el eje "x"
y = 0:0.1:4;           % Intervalo de puntos en el eje "y"
[X,Y] = meshgrid(x,y); % Creación de las matrices que representan 
                       % los puntos de la malla
                       
% El vector de posición para cada punto inicial de la placa queda definido
% por la matriz de coordenadas
% Los desplazamientos se producen en la dirección del vector unitario i,
% esto es en la dirección del eje "x"

U = (Y.^2)./20; % Es la matriz de los desplazamientos para el eje "x"
R = X + U; % Es la matriz con las coordenadas en x, después del desplazamiento

figure(1);

subplot(1,2,1);
mesh(X,Y,0*X);
view(2);
axis([-1.5 1.5 -1 5]);
title('Puntos de la placa sin deformación');
xlabel('Eje "x"');
ylabel('Eje "y"');

subplot(1,2,2);
mesh(R,Y,0*R);
view(2);
axis([-1.5 1.5 -1 5]);
title('Puntos después de la transformación');
xlabel('Eje "x"');
ylabel('Eje "y"');

figure(2);
mesh(X,Y,0*X);
hold on;
mesh(R,Y,0*R);
view(2);
title('Comparamos la placa antes y depués de la transformación');
xlabel('Eje "x"');
ylabel('Eje "y"');

figure(3);
mesh(X,Y,0*X);
hold on;
quiver3(X,Y,0*X,U,0*Y,0*U,5);
view(2);
title('Campo vectorial de desplazamientos sobre la placa');
xlabel('Eje "x"');
ylabel('Eje "y"');

6 Divergencia del campo Desplazamientos

La Divergencia de un campo vectorial, en nuestro caso del campo de desplazamientos, está definida por la siguiente expresión:

[math]∇.\vec U= \frac{∂U_1}{∂x} i ⃗+ \frac{∂U_2}{∂y}j ⃗+ \frac{∂U_3}{∂z} \vec k[/math]

Geométrica-mente la divergencia de un campo vectorial es una medida del cambio de volumen local debido al desplazamiento, en este caso, el desplazamiento de la placa rectangular en el plano XY produce un cambio de volumen nulo, luego la divergencia ha de ser cero.

Además, podemos comprobar en la definición de divergencia que para nuestro caso, efectivamente toma valor nulo, ya que el desplazamiento se produce únicamente en la dirección de “x” y además, éste sólo depende de los valores de “y”, luego la derivada es cero.

7 Rotacional del campo Desplazamientos

El rotacional de un campo vectorial está definido por la expresión:

[math]\nabla \times \vec U = \begin {bmatrix} \vec i & \vec j & \vec k \\ \frac{\delta}{\delta x} & \frac{\delta}{\delta y} & \frac{\delta}{\delta z} \\ \ U_1 & U_2 & U_3 \end{bmatrix} = \begin {bmatrix} \vec i & \vec j & \vec k \\ \frac{\delta}{\delta x} & \frac{\delta}{\delta y} & \frac{\delta}{\delta z} \\ \ y^2/20 & 0 & 0 \end{bmatrix} = -\frac{\delta \frac{y^2}{20}}{\delta y} \vec k = \frac{-y}{10} \vec k[/math]

Para representar el rotacional se ha de hallar su módulo, en este caso el módulo del campo vectorial U sólo tiene componente para el eje "z", luego se aplica el valor absoluto:

[math]|\nabla \times \vec U| = |\frac{-y}{10}| = \frac{y}{10}[/math]

% Rotacional del campo vectorial de desplazamientos

x = -0.5:0.1:0.5;      % Intervalo de puntos en el eje "x"
y = 0:0.1:4;           % Intervalo de puntos en el eje "y"
[X,Y] = meshgrid(x,y); % Creación de las matrices que representan 
                       % los puntos de la malla
                       
rz = inline('-y/10','x','y');
Zr = abs(rz(X,Y));

subplot(1,2,1);
surf(X,Y,Zr);
hold on;
mesh(X,Y,0*X);
title('Representación del rotacional');
ylabel('Eje "y"');
xlabel('Eje "x"');

subplot(1,2,2);
surf(X,Y,Zr);
view(2);
axis([-1.5 1.5 -1 5]);
title('Vista ortogonal del rotacional sobre la placa');
ylabel('Eje "y"');
xlabel('Eje "x"');

Los puntos con que sufren mayor rotacional son los [math](x,y) = (x,4)[/math]. Cuyo valor de rotacional es 0.4

8 Tensiones debidas al desplazamiento

8.1 Tensiones normales

Sabiendo que el medio que tratamos, es un medio elástico lineal, isótropo y homogéneo, entonces el tensor de tensiones sobre el sólido elástico depende de los coeficientes de Lamé (\(\lambda\) y \(\mu\)):

\( \sigma_{ij}=\lambda \nabla \cdot \vec u \delta_{ij} + 2\mu \epsilon_{ij} \)

Para el cálculo de las tensiones tomamos como valor de éstos coeficientes 1.

Las tensiones normales se calculan con la fórmula anterior.Como la divergencia del campo vectorial es cero, y como los coeficientes de Lamé son iguales a 1, la tensión normal es \( \sigma_{ij}=2*\epsilon_{ij} \), donde \(\epsilon_{ij} \) es la parte simétrica del tensor gradiente de \(\vec u\).

Así, tenemos que la tensión normal al eje \(\vec i\) nos queda \(\vec i\)*\(\sigma\)*\(\vec i\), que es la proyección de

\(\sigma\)*\(\vec i\) (que depende solo de \(\vec j\)) sobre \(\vec i\), dando como resultado cero.

Para la tensión normal al eje \(\vec j\), ocurre lo mismo, siendo \(\vec j\)*\(\sigma\)*\(\vec j\) la proyección de

\(\sigma\)*\(\vec j\) (que depende solo de \(\vec i\)) sobre \(\vec j\), dando también cero.

8.2 Tensiones tangenciales respecto al plano ortogonal a î

Las tensiones tangenciales están definidas por la ecuación \(|\sigma \cdot \vec i-(\vec i \cdot \sigma \cdot \vec i) \vec i|\) respecto al plano ortogonal a \(\vec i\). El segundo sumando es cero. Entonces queda \(|\sigma \cdot \vec i|\) = [math][0,\frac{y}{10},0][/math]

% Tensiones tangenciales respecto al plano ortogonal a î

x = -0.5:0.1:0.5;      % Intervalo de puntos en el eje "x"
y = 0:0.1:4;           % Intervalo de puntos en el eje "y"
[X,Y] = meshgrid(x,y); % Creación de las matrices que representan 
                       % los puntos de la malla
                       
Fx = 0*X; % Componentes x del campo vectorial
Fy = Y./10; % Componentes y

subplot(1,2,1);
quiver(X,Y,Fx,Fy,2.5);
axis([-1.5 1.5 -1 5]);
view(2);
title('Tensiones tangenciales plano ortogonal a i');
xlabel('Eje "x"');
ylabel('Eje "y"');

subplot(1,2,2);
quiver(X,Y,Fx,Fy,2.5);
hold on;
mesh(X,Y,0*X);
axis([-1.5 1.5 -1 5]);
view(2);
title('Visualización de las tensiones en la placa')
xlabel('Eje "x"');
ylabel('Eje "y"');

8.3 Tensiones tangenciales respecto al plano ortogonal a ĵ

Las tensiones están definidas por la ecuación \(|\sigma \cdot \vec j-(\vec j \cdot \sigma \cdot \vec j) \vec j|\) respecto al plano ortogonal a \(\vec j\). El segundo sumando es cero por ello queda \(|\sigma \cdot \vec j|\) = [math][\frac{y}{10},0,0][/math]

% Tensiones tangenciales respecto al plano ortogonal a j

x = -0.5:0.1:0.5;      % Intervalo de puntos en el eje "x"
y = 0:0.1:4;           % Intervalo de puntos en el eje "y"
[X,Y] = meshgrid(x,y); % Creación de las matrices que representan 
                       % los puntos de la malla
                       
Fx = Y./10; % Componentes x del campo vectorial
Fy = 0*Y; % Componentes y

subplot(1,2,1);
quiver(X,Y,Fx,Fy,2.5);
axis([-1.5 1.5 -1 5]);
view(2);
title('Tensiones tangenciales plano ortogonal a j');
xlabel('Eje "x"');
ylabel('Eje "y"');

subplot(1,2,2);
quiver(X,Y,Fx,Fy,2.5);
hold on;
mesh(X,Y,0*X);
axis([-1.5 1.5 -1 5]);
view(2);
title('Visualización de las tensiones en la placa')
xlabel('Eje "x"');
ylabel('Eje "y"');

8.4 Tensión de Von Mises

La tensión de Von Mises se define por la fórmula

[math] \sigma _{VM}=\sqrt{\frac{(\sigma_1-\sigma_2)^{2}+(\sigma_2-\sigma_3)^{2}+(\sigma_3-\sigma_1)^{2}}{2}}[/math]

[math]\sigma_1, \sigma_2, \sigma_3[/math] Son los autovalores de la matriz tensor de tensiones σ (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.

Los valores obtenidos son: [math]\sigma_1 = 0 , \sigma_2 = \frac{y}{10} , \sigma_3 = \frac{-y}{10}[/math]

Una vez obtenidos los autovalores, la tensión de Von Mises sobre la placa queda representada así:

% Intervalos de trabajo
x=[-1/2:0.1:1/2];
y=[0:0.1:4];
% Mallado
[X,Y]=meshgrid(x,y);
[m,n]=size(X);
% Definimos la matriz asociada al tensor de tensiones
sigma=[0,1/10,0;1/10,0,0;0,0,0];

% Autovalores 

auto = eig(sigma);

t1 = Y.*auto(1);
t2 = Y.*auto(2);
t3 = Y.*auto(3);

% Tensión de Von Mises
TVM=sqrt(((t1-t2).^2+(t2-t3).^2+(t3-t1).^2)/2);


subplot(2,1,1)
mesh(X,Y,TVM)
subplot(2,1,2)
surf(X,Y,TVM)

Tensión de Von Mises máxima en y=4

9 Masa total de la placa y función de densidad

Sabemos que la función que define la densidad de la placa está dada por la expresión:

[math]f_d(x,y) = x.y.e^{-x^-2}[/math]

Vamos a representar dicha función para ver como se distribuye en la placa la densidad.

% Representación de la función densidad

x = -0.5:0.1:0.5;      % Intervalo de puntos en el eje "x"
y = 0:0.1:4;           % Intervalo de puntos en el eje "y"
[X,Y] = meshgrid(x,y); % Creación de las matrices que representan 
                       % los puntos de la malla
                       
Zd = (X.*Y.*exp(-X.^(-2)));
mesh(X,Y,Zd);
title('Función de densidad de la placa');
xlabel('Eje "x"');
ylabel('Eje "y"');

Al no tener sentido hablar de densidad negativa, se toman como absolutos todos los valores de la función de densidad.

Así la masa total de la placa es:

[math]M = \int_{S}f_d dS[/math]

La integral define el volumen bajo la superficie, cuyo valor corresponde a la masa total de la placa.

Para el cálculo del volumen vamos a proceder por el método de aproximación numérica que considera los volúmenes de los paralelepípedos, cuyas bases tienen área [math]H^2 ;H=0.1[/math] y volúmenes las imágenes de la función de densidad. La suma de todos los volúmenes nos da el valor deseado.

% Cálculo de la masa total de la placa a partir de la función de densidad

x = -0.5:0.1:0.5;      % Intervalo de puntos en el eje "x"
y = 0:0.1:4;           % Intervalo de puntos en el eje "y"
[X,Y] = meshgrid(x,y); % Creación de las matrices que representan 
                       % los puntos de la malla
                       
Zd = (X.*Y.*exp(-X.^(-2))); % Función de densidad
V = 0.1^2*Zd; % Matriz cuyas componentes son los volúmenes de los 
% paralelepípedos en cuestión
V = abs(V); % Toma valores absolutos de las componentes de la matriz
M = sum(sum(V)); % La masa total es la suma de todos los volúmenes
disp(M); % Muestra en pantalla la masa total de la placa

La masa total de la placa tiene un valor de [math]M = 0.0163[/math]