Visualización de campos escalares y vectoriales en elasticidad (Grupo A12)

Trabajo realizado por estudiantes
Título	Visualización de campos escalares y vectoriales en elasticidad. Grupo 12-A
Asignatura	Teoría de Campos
Curso	2014-15
Autores	Adrián Pascual Blázquez, Alejandro Giménez Alves, Miguel Aparicio Martín-Romo, Juan José Olivas Caballero
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Enunciado

Visualización de campos escalares y vectoriales en elasticidad. Consideramos una placa rectangular plana (de dimensión 2) que ocupa la región [-1/2, 1/2]x[0,4]. En ella vamos a suponer que tenemos definidas dos cantidades físicas: la temperatura T(x,y), que depende de las dos variables espaciales (x,y), y los desplazamientos [math]\vec u(x,y,t)[/math] producidos por la acción de una fuerza determinada.

De esta forma, si definimos [math]\vec 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 (x,y) de la placa después de la deformación viene dada por:

[math]\vec r(x,y) = \vec r_{0}(x,y) + \vec u(x,y)[/math]

Vamos a suponer que la fuerza aplicada sobre la placa ha provocado un desplazamiento de los puntos de la misma dado por el vector de desplazamientos:

[math]\vec u(x,y)= \vec a(\vec b · \vec r_{0})^2 [/math] donde [math]\vec a [/math] y [math]\vec b [/math] son vectores dados.

Por último supondremos:

[math] \vec a = \frac{\vec i}{20} , \vec b = \vec j [/math]

En este caso, [math]\vec u(x,y)=\frac{y^2}{20} \vec i[/math].

2 Representación del mallado del sólido

Para la representación del mallado utilizaremos el programa Matlab. En el introduciremos el comando Linspace para representar los intervalos, dónde tomaremos como paso de muestreo [math] h = \frac{1}{10} . [/math] La malla la creamos con el comando Meshgrid. Como la placa es plana, deberemos tomar una matriz Z de ceros. Por último, para su representación utilizaremos el comando surf.

x=linspace(-1/2,1/2,11);  %Límites de la placa en el eje X
y=linspace(0,4,41);       %Límites de la placa en el eje Y
[X,Y]=meshgrid(x,y);      %Mallado de la placa
Z=zeros(41,11);           %Altura de cada punto de la placa
surf(X,Y,Z)               %Representación de la placa
axis([-1,1,-1,5])         %Límites de representación
view(2)                   %Punto de vista vertical

Placa horizontal [(-0.5,0.5) x (0,4)]

3 Representación de un campo escalar

Se pide la representación de la temperatura (campo escalar) sobre una placa. La temperatura viene definida por la función:

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

Utilizamos la función inline. Aplicando esta función a las matrices de coordenadas de la placa obtenemos Z1, que es la matriz imagen de la función T. Por último, representamos la superficie con el comando surf y calculamos el valor máximo de la temperatura sobre la placa.

x=linspace(-1/2,1/2,11);                                    %Límites de la placa en el eje X
y=linspace(0,4,41);                                         %Límites de la placa en el eje Y
[X,Y]=meshgrid(x,y);                                        %Mallado de la placa
f=inline('(8-y.^2+2.*y).*exp(-(x+1).^2) ','x','y');         %Función de temperatura
Z1=f(X,Y);                                                  %Matriz de la temperatura de los puntos de la placa
surf(X,Y,Z1)                                                %Representación de la temperatura en la placa
axis([ -1,1 , -1,5 ])                                       %Límites de representación
view(2)                                                     %Punto de vista vertical
max(max(Z1))                                                %Valor máximo de la temperatura

El valor máximo de la temperatura resulta 7,0092.

Temperatura con foco cercano a la esquina izquierda

4 Cálculo y representación del gradiente y curvas de nivel de una función escalar.

Como la temperatura de la placa esta definida por la función:

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

Calculamos el gradiente con:

[math] \nabla T(x,y,z)= \frac{\partial T}{\partial x} \vec i + \frac{\partial T}{\partial y} \vec j + \frac{\partial T}{\partial z} \vec k = -2(x+1)(8-y^2+2y)*e^{-(x+1)^2} \vec i + (-2y+2)*e^{-(x+1)^2} \vec j [/math]

Una curva de nivel es aquella línea que en una figura, une todos los puntos que tienen igualdad de altura, para representarlas, se usará el comando “contour”. Procedemos pues al cálculo y representación del gradiente, para ello se han de calcular las derivadas parciales con respecto a X e Y de la función T, en este caso, serán las funciones dx, dy, las cuales son traducidas a matlab con “inline”. Aplicando las matrices del mallado a estas funciones, se obtienen las matrices de las coordenadas de el campo vectorial gradiente de T, las cuales pueden ser representadas en el recinto de la placa, (limitado por las matrices del mallado) con el comando “quiver”.

x=linspace(-1/2,1/2,11);                                              %Intervalos de trabajo para las curvas de nivel
y=linspace(0,4,41);
f=inline('(8-y.^2+2.*y).*exp(-(x+1).^2)','x','y');                    %Formula de temperatura y sus derivadas parciales
fx=inline('-2.*(x+1).*(8-y.^2+2.*y).*exp(-(x+1).^2)','x','y');
fy=inline('(-2.*y+2).*exp(-(x+1).^2)','x','y');
[X,Y]=meshgrid(x,y);                                                  %Mallado para las curvas de nivel
T=f(X,Y);
xx=linspace(-1/2,1/2,11);                                             %Intervalos para el gradiente
yy=linspace(0,4,41);
[XX,YY]=meshgrid(xx,yy);                                              %Mallado para las curvas de nivel
U=fx(XX,YY);                                                          %Componentes del gradiente
V=fy(XX,YY);
%%%%%%%                                                               %Dibujo del gradiente con las curvas de nivel
contour(X,Y,T);                                                       %Curvas de nivel
hold on
quiver(XX,YY,U,V)                                                     %Gradiente

Lineas de nivel y campo gradiente de la temperatura

5 Representación de campos vectoriales

Ahora se va a suponer que sobre la placa se ha aplicado una fuerza que ha provocado el desplazamiento de los puntos de la misma. El desplazamiento viene dado en cada punto por el vector:

[math]\vec u(x,y)=\frac{y^2}{20} \vec i[/math]

En este caso, el campo vectorial [math]\vec u [/math] solo está definido en la dirección [math]\vec i[/math]. Esto quiere decir que únicamente hay desplazamiento de los puntos de la placa paralelamente al eje X. Como no hay desplazamiento de los puntos de la placa en la dirección del eje Y, el vector que definiría dichos desplazamientos vendría definido por [math]\vec v(x,y)=0[/math].

x=linspace(-1/2,1/2,11);       
y=linspace(0,4,41);
[X,Y]=meshgrid(x,y);       %Mallado de la placa   
u=(Y.^2)/20;               %Desplazamientos en la dirección del eje x
v=0*X;                     %Desplazamientos en la dirección del eje y
quiver(x,y,u,v)            %Representación del campo vectorial de desplazamientos

Vector de desplazamientos

6 Representación de cada punto del sólido antes y despúes del desplazamiento determinado por el campo vectorial [math]\vec u[/math].

Definimos el vector de posición de los puntos de la placa, después de aplicar la fuerza, según la expresión [math]\vec r (x,y)= \vec r_{0}(x,y)+\vec u(x,y)[/math], donde [math]vec r_0(x,y)[/math] es el vector de posición de los puntos de la placa en reposo y [math]vec u(x,y)[/math] el vector desplazamiento definido antes.

xc=linspace(-1/2,1/2,11);            
yc=linspace(0,4,41);                 
[XC,YC]=meshgrid(xc,yc);             %Mallado de la placa
X=XC+(YC.^2)/20;                     %Ecuaciones de los puntos desplazados
Y=YC;
figure(1)                            %Dibujo de las placas
subplot(1,2,1)                       %Placa sin desplazamiento
mesh(XC,YC,0*XC)
axis([-1.5,1.5,-1,5])
subplot(1,2,2)                       %Placa con desplazamiento
mesh(X,Y,0*X)                        %Suponemos que el campo de desplazamientos no tiene componente en z
axis([-1.5,1.5,-1,5])
view(2)

Placa horizontal [(-0.5,0.5) x (0,4)]

7 Calculo y representación de la divergencia en todos los puntos del sólido

Tenemos [math]\vec u(x,y) = \vec a (\vec b · \vec r_0)^2 [/math] y [math] \vec r_0 (x,y)= x \vec i + y \vec j [/math]

Luego [math]\vec u(x,y) = \frac{\vec i}{20} · (\vec j ·(x \vec i + y \vec j))^2 = \frac{\vec i}{20} · y^2 = \frac{y^2}{20} \vec i [/math]

En coordenadas cartesianas, podemos definir la divergencia como:

[math]\nabla\cdot\vec u= \frac{ \partial \vec u_x }{ \partial x} + \frac{ \partial \vec u_y }{ \partial y}+ \frac{ \partial \vec u_z }{ \partial z}=0 [/math]

La divergencia de un campo vectorial es una medida del cambio de volumen local debido al desplazamiento, en este caso, la divergencia es 0, por tanto la representación es 0 en todos los puntos. Lo comprobamos visualmente en el gráfico obtenido con matlab:

x=linspace(-1/2,1/2,11);           %Intervalos de trabajo
y=linspace(0,4,41);
[X,Y]=meshgrid(x,y);               %Mallado
u=(Y.^2)/20;                       %Componentes del campo vectorial
v=0.*Y;
div=divergence(X,Y,u,v);           %Divergencia
mesh(X,Y,div)
axis([-1,1,-0.5,4.5])
view(2)

Representación de la divergencia en todos los puntos de la placa

8 Cálculo del rotacional

En coordenadas cartesianas, se puede definir el rotacional como:

[math]\nabla\times\vec{u}= \begin{pmatrix} \vec{i} & \vec{j} & \vec{k} \\ \frac{ \partial}{\partial x} & \frac{\partial }{\partial y} & \frac{\partial }{\partial z}\\ u_x & u_y & u_z \end{pmatrix}= \begin{pmatrix} \vec{i} & \vec{j} & \vec{k} \\ \frac{ \partial}{\partial x} & \frac{\partial }{\partial y} & \frac{\partial }{\partial z}\\ \frac{y^2}{20}& 0 & 0 \end{pmatrix}= - \frac{y}{10} \vec k [/math]

Para representar el rotacional hay que hallar su módulo. En este caso, el módulo del campo vectorial U sólo tiene componente en z, luego se aplica el valor absoluto:

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

Introducimos con el comando “inline” la componente Z del rotacional, y la aplicamos a los puntos del sólido.

x=linspace(-1/2,1/2,11);                   %Límites de la placa en el eje X
y=linspace(0,4,41);                        %Límites de la placa en el eje Y
[X,Y]=meshgrid(x,y);                       %Mallado de la placa
Z=0*X;                                     %Definición del eje Z
hz=inline('-y/10','x','y','z');            %Función del rotacional en cada punto de la placa
Z2=abs(hz(X,Y,Z));                         %Módulo del rotacional
surf(X,Y,Z2);                              %Representación del módulo del rotacional en cada punto de la placa
axis([-1,1,-1,5])                          %Límites de representación
view(2)

Representación del modulo del rotacional del campo U en cada punto de la placa

Se puede observar que los puntos que sufren un mayor rotacional (tendencia de un campo al giro) son los de mayor y (y=4). [math]|\nabla\times\vec{u}|= \frac{y}{10} = 0.4 [/math]

9 Tensiones

Definimos la función épsilon como la parte simétrica del tensor gradiente de [math]\vec u[/math] como [math] \epsilon (\vec u) = \frac{(\nabla \vec u + \nabla \vec u ^t)}{2} [/math] siendo [math] \nabla \vec u = \frac{\delta u_i}{\delta x^j} \vec e_i \otimes \vec e_j = \frac{y}{10} \vec i \otimes \vec j [/math] y su transpuesto [math] \nabla \vec u^t = \frac{y}{10} \vec j \otimes \vec i [/math].

Por lo tanto [math] \epsilon (\vec u) = \frac{y}{20} (\vec i \otimes \vec j + \vec j \otimes \vec i [/math]) = \begin{pmatrix} 0 & \frac{y}{20} & 0 \\ \frac {y}{20} & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}

Se define el tensor de tensiones a través de la fórmula [math] \sigma = \lambda \nabla \cdot \vec u \mathbb{1} + 2 \mu \epsilon [/math] , la cual tomando los valores [math] \lambda = \mu \ = 1 [/math] y como [math] \nabla \cdot \vec u = 0 [/math]:

[math] \sigma = 2 \epsilon = \frac{y}{10} (\vec i \otimes \vec j + \vec j \otimes \vec i) [/math]

9.1 Tensiones normales

9.1.1 Normales a [math] \vec i[/math]

Las tensiones normales con respecto a [math] \vec i [/math] se definen como, [math] \vec i \cdot \sigma \cdot \vec i [/math], operando con el tensor [math] \sigma = 2 \epsilon = \frac{y}{10} (\vec i \otimes \vec j + \vec j \otimes \vec i) [/math]: [math] \vec i \cdot \sigma \cdot \vec i = 0 [/math]

9.1.2 Normales a [math] \vec j[/math]

Las tensiones normales con respecto a [math] \vec j [/math] se definen como [math] \vec j \cdot \sigma \cdot \vec j [/math], operando igualmente con el tensor [math] \sigma = 2 \epsilon = \frac{y}{10} (\vec i \otimes \vec j + \vec j \otimes \vec i) [/math]: [math] \vec j \cdot \sigma \cdot \vec j = 0 [/math]

9.2 Tensiones tangenciales

9.2.1 Tensiones tangenciales respecto al plano ortogonal a [math]\vec i[/math].

Para hallar las tensiones tangenciales al plano ortogonal a [math]\vec i[/math], usamos la fórmula [math]|\sigma \cdot \vec i-(\vec i \cdot \sigma \cdot \vec i) \vec i|[/math] para calcularlas, obteniendo una tensión dada por un vector en la dirección [math]\vec j[/math]. Construimos las matrices Kx, Ky y Kz con las componentes de las tensiones, y después las dibujamos usando el comando quiver3 para representar campos vectoriales en 3 dimensiones.

x=linspace(-1/2,1/2,11);       %Límites de la placa en el eje X
y=linspace(0,4,41);            %Límites de la placa en el eje Y
[X,Y]=meshgrid(x,y);           %Mallado de la placa
Z=zeros(41,11); 
Ky=Y/10;                       %Componente Y de las tensiones tangenciales respecto al plano ortogonal a i
Kx=zeros(size(Ky));            %Componente X de las tensiones tangenciales respecto al plano ortogonal a i             
Kz=zeros(size(Ky));            %Componente Z de las tensiones tangenciales respecto al plano ortogonal a i
quiver3(X,Y,Z,Kx,Ky,Kz);       %Representación del campo de las tensiones tangenciales respecto al plano ortogonal a i
axis([-1,1,-1,5])              %Límites de representación
view(2)                        %Punto de vista vertical

Representación de las tensiones tangenciales provocadas por el campo U en el eje x

9.2.2 Tensiones tangenciales respecto al plano ortogonal a [math]\vec j[/math].

Para obtener y dibujar las tensiones tangenciales respecto al plano ortogonal a [math]\vec j[/math], el método es análogo al de las tensiones tangenciales respeto al plano ortognal a [math]\vec i[/math], con la correspondiente fórmula [math]|\sigma \cdot \vec j-(\vec j \cdot \sigma \cdot \vec j) \vec j|[/math]. En este caso las tensiones obtenidas vienen representadas por un vector en la dirección [math]\vec i[/math]. Para representarlas lo haremos igual también, construyendo matrices con las componentes de las tensiones y representándolas sobre los puntos de coordenadas X, Y, Z con el comando quiver3.

x=linspace(-1/2,1/2,11);       %Límites de la placa en el eje X
y=linspace(0,4,41);            %Límites de la placa en el eje Y
[X,Y]=meshgrid(x,y);           %Mallado de la placa
Z=zeros(41,11); 
Kx=Y/10;                       %Componente X de las tensiones tangenciales respecto al plano ortogonal a j                  
Ky=zeros(size(Kx));            %Componente Y de las tensiones tangenciales respecto al plano ortogonal a j
Kz=zeros(size(Kx));            %Componente Z de las tensiones tangenciales respecto al plano ortogonal a j
quiver3(X,Y,Z,Kx,Ky,Kz);       %Representación del campo de las tensiones tangenciales respecto al plano ortogonal a j
axis([-1,1,-1,5])              %Límites de representación
view(2)                        %Punto de vista vertical

Representación de las tensiones tangenciales provocadas por el campo U en el eje y

10 Representación de la tensión de Von Mises

La tensión de Von Mises viene definida por la fórmula:

[math] σ_{VM}=\sqrt[]{\frac{(σ_1-σ_2)^2+(σ_2-σ_3)^2+(σ_3-σ_1)^2}{2}}[/math]

donde [math] σ_1 [/math], [math] σ_2 [/math] y [math] σ_3 [/math], también conocidos como tensiones principales, son los autovalores del tensor de tensiones [math] σ [/math]. Se trata de una magnitud escalar que indica cuándo un material inicia comportamiento plástico (y no elástico puro).

sigmatotal=[];                                                          %Vector vacío para la primera interacción.
for y=linspace(0,4,41);                                                 %Bucle para crear el vector sigmatotal
    sigma=[0 y/10 0;y/10 0 0;0 0 0];                                    %Matriz sigma
    E=eig(sigma);                                                       %Función para sacar los autovalores de la matriz sigma
    sigmavm= sqrt(((E(1)-E(2))^2+(E(2)-E(3))^2+(E(3)-E(1))^2)/2);       %Ecuación de la tensión de Von Mises
    sigmatotal=[sigmatotal,sigmavm]                                     %Creación del vector sigmatotal a partir de las tensiones de las interacciones anteriores
end
plot(sigmatotal)                                                        %Dibujo de la tensión

Representación de la tensión de Von Mises

Se puede observar que los puntos que alcanzan un valor mayor son los de mayor y (y=4), cuyo valor es [math] σ_{VMmax}={\frac{\sqrt{3}*y}{10}}=0.693 [/math] que se sacaría añadiendo en matlab el comando max(max(sigmatotal))

11 Masa de la placa

11.1 Método del trapecio para integrales dobles

Consideramos [math] [a,b]\times [c,d] [/math] un rectángulo y [math] f:[a,b]\times[c,d]\to \mathbb{R} [/math] una función real. Para aproximar la integral [math] \int_a^b\int_c^df(u) \; dv \; du [/math], aplicamos el método del trapecio. En primer lugar, consideramos la partición del intervalo [math] [a,b] [/math] en [math] N_1 [/math] subintervalos de longitud [math] h_1=\frac{b-a}{N_1}[/math]. Definimos [math] u_n=a+nh, [/math] dónde [math] n=0,1,...,N_1 [/math]. El método del trapecio nos da:

[math] \int_a^b \int_c^d f(u,v) \; du \; dv \sim h_1\sum_{n=0}^{N}w_{n}\int_c^df(u_n,v)dv. [/math]

Ahora, usamos otra vez el método del trapecio para resolver la integral restante: [math] \int_a^b \int_c^d f(u,v) \; du \; dv \sim h_1h_2\sum_{n=0}^{N}\sum_{m=0}^{N}w_{n}\hat w_mf(u_n,v_m)dv =h_1h_2w^t\cdot f\cdot \hat w. [/math]

donde [math] f [/math] es la [math] N_1\times N_2 [/math] matriz de componentes [math] f_{ij}=f(u_i,v_j)[/math] y [math] \hat w [/math] es similar a [math] w [/math] pero con [math] N_2+1 [/math] filas.

11.1.1 Cálculo de la masa de la placa

En Matlab/Octave aproximamos la integral con la función de densidad [math] f(u,v)= u·v·e^{-1/u^2} [/math] en el intervalo [math] [-0.5,0.5]\times [0,4] [/math]

N1=200; N2=100;                       %Número de puntos
a=-0.5; b=0.5; c=0; d=4;              %Extremos del intervalo/Valores de integración
h1=(b-a)/N1; h2=(d-c)/N2;
u=a:h1:b; v=c:h2:d;                   %Coordenadas de la partición
[uu,vv]=meshgrid(u,v);                %Coordenadas del rectángulo
f=abs((uu.*vv).*exp(-1./(uu.^2)));    %Función en valor absoluto debido a que toma valores negativos en algunos puntos
w1=ones(N1+1,1);                      %Vector
w1(1)=1/2; w1(N1+1)=1/2;
w2=ones(N2+1,1);                      %Vector
w2(1)=1/2; w2(N2+1)=1/2;
masa=h1*h2*w2'*f*w1                   %Valor de la masa

Por lo tanto, el resultado que nos da Matlab para la masa de la placa es:

[math] \int_a^b \int_c^d f(u,v) \; du \; dv = \int_a^b \int_c^d u*v*e^{-1/u^2} \; du \; dv = 0.0064 [/math]