Visualización de campos vectoriales y escalares en elasticidad. grupo 3 A

Trabajo realizado por estudiantes
Título	Visualización de campos escalares y vectoriales en elasticidad.Grupo A3
Asignatura	Teoría de Campos
Curso	2014-15
Autores	Miguel Domínguez Trufero (955) Jorge Bautista Mendo (991) Enrique Crespo Ferrer (1015) Jaime Chueca Rincón (1019) Miguel Fandiño Álvarez (1113) Jorge García Saenz (1126)
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

Consideramos una placa circular comprendida entre las circunferencias de radio 1 y radio 2, centradas en el origen. En ella vamos a suponer que tenemos definidas dos cantidades físicas: la temperatura \(T(\rho,\theta,t)\) y el tiempo \(t\), y los desplazamientos \(\vec u(\rho,\theta,t)\).De esta forma, si definimos \(r_0(\rho,\theta)\) el vector de posición de los puntos de la placa en reposo, la posición de cada punto \((\rho,\theta)\) de la placa en un instante de tiempo \(t\) viene dada por\[ \vec r (\rho,\theta,t)= \vec r_{0}(\rho,\theta)+\vec u(\rho,\theta,t). \] Vamos a suponer que sobre la placa se ha aplicado una fuerza que ha provocado una vibración de manera que los desplazamientos en un instante \(t_0\) vienen dados por: [math]\vec u(\rho,\theta)=(1-\rho)^2 \vec g_{\theta}. [/math]

1 Mallado del sólido

Mallado de la corona circular.

El mallado del sólido representa el dominio en el cual se encuentran los puntos que componen el sólido. Por lo tanto, para realizar el mallado de deben crear dos matrices cuyos componentes serán distintos valores que pueden tomar cada una de las variables (ro, theta). Después por interpolación hallaremos el cojunto total de los puntos que componen el sólido,tomando como paso de muestreo h=1/10 y procedemos a dibujar el mallado:

h=0.1                                % Establecemos paso de muestra
ro=1:h:2;                            % Valor de ro según el muestreo del intervalo
teta=0:h:2*pi+h;                     % Valor de teta según el muestreo del intervalo
[roro,tetateta]=meshgrid(ro,teta);   % Matrices de la placa en coordenadas polares
figure(1)
xx=roro.*cos(tetateta);              % Parametrización
yy=roro.*sin(tetateta);
mesh(xx,yy,0*xx)                     % Dibujar la malla
axis([-3,3,-3,3])                    % Seleccionar región 
view(2)                              % Visualización

2 Temperatura del sólido

Distribución de la temperatura en la placa

La temperatura del sólido varía a lo largo de este según la función [math]T(x,y)=e^{-y}[/math]. Esta depende solo de la variable Y, por lo que la X se mantendrá constante.

h=0.1                              % Establecemos paso de muestra
ro=1:h:2;                          % Valor de ro según el muestreo del intervalo
teta=0:h:2*pi+h;                   % Valor de teta según el muestreo del intervalo
[roro,tetateta]=meshgrid(ro,teta); % Matrices de la placa en coordenadas polares
figure(1)
xx=roro.*cos(tetateta);            % Parametrización
yy=roro.*sin(tetateta);
f=exp(-yy);                        % Campo escalar
surf(xx,yy,f)                      % Dibuja la malla
axis([-2,2,-2,2])                  % Selecciona los limites de los ejes x e y en representación de la placa
colorbar                           % Ponemos en la figura la leyenda de los colores
view(2)                            % Visualización

Como podemos observar en la gráfica, la temperatura desciende de forma negativa en el eje de ordenadas,que es donde se produce la máxima variación,y ya que la X se mantiene igual podemos comprobar que la temperatura es máxima en el punto (x,y)=(0,-2). También podemos observar que las curvas de nivel referentes a la superficie de la temperatura no son equidistantes porque la superficie no tiene una pendiente lineal,de forma que cuanto mayor es el gradiente mas cercanas se encontraran.

3 Gradiente del campo de temperaturas

El gradiente es el campo vectorial cuyas funciones coordenadas son las derivadas parciales del campo escalar dado, que representa la temperatura en cada punto de nuestro sólido. En este caso, el gradiente en cada uno de los puntos apunta en la dirección en la cual la temperatura aumenta más rápido (es decir, la dirección en que la derivada direccional es máxima). El módulo del gradiente indica la rapidez con que aumenta la temperatura en esa dirección. En este caso el gradiente del campo es [math]\nabla T = (0,-e^{-y})[/math] Como sabemos por las propiedades del gradiente, el vector gradiente será perpendicular a las líneas de contorno (líneas "equiescalares", definidas por ϕ =cte) del sólido. Podemos observarlo fácilmente en la figura en 3D del sólido que se muestra más abajo.

Campo vectorial del gradiente de T

h=0.1                         % Establecemos paso de muestra
ro=1:h:2;                     % Valor de ro según el muestreo del intervalo
teta=0:h:2*pi+h;              % Valor de teta según el muestreo del intervalo
[roro,tetateta]=meshgrid(ro,teta); % Matrices de la placa en coordenadas polares
figure(1)
x=roro.*cos(tetateta);        % Parametrization
yy=roro.*sin(tetateta);
hold on                       % Para dibujar el campo escalar con sus lineas equiescalares y su gradiente en la misma figura
f=exp(-yy);                   % El campo escalar
contour(xx,yy,f,30)           % Dibuja la malla
axis([-2,2,-2,2])             % Selecciona la región para dibujar(Dimensiones de los ejes)
% El gradiente de T viene dado por el campo -exp(-y) según el versor j
fx=xx*0;                      % x-Introducimos la función de la derivada parcial respecto de x 
fy=-exp(-yy);                 % y-Introducimos la función de la derivada parcial respecto de y 
quiver(xx,yy,fx,fy)           % Dibujamos el campo vectorial del gradiente de temperaturas
axis([-2,2,-2,2])             % Selecciona la región para dibujar(Dimensiones de los ejes)
view(2)                       % Visualización
colorbar                      % Ponemos en la figura la leyenda de los colores
hold off                      % Para dibujar el campo escalar con sus lineas equiescalares y su gradiente en la misma figura

La temperatura es máxima en el punto (x,y)=(0,-2)

4 Campo de desplazamientos u

Ahora trabajaremos con un campo de desplazamientos [math]\vec u(\rho,\theta)=(1-\rho)^2 \vec g_{\theta}. [/math]. Dibujaremos su correspondiente campo de vectores, y lo introducimos en los puntos del mallado del sólido del apartado 1.

Campo vectorial del desplazamiento [math]\vec u(\rho,\theta)=(1-\rho)^2 \vec g_{\theta}. [/math]

h=0.1;                              % Establece paso de muestra
ro=1:h:2;                           % Establece intervalo de ro
teta=0:h:2*pi+h;                    % Establece intervalo de teta 
[roro,tetateta]=meshgrid(ro,teta);  % Matriz en coordenadas polares
figure(1)
xx=roro.*cos(tetateta);             % Parametrización
yy=roro.*sin(tetateta); 
fx=((1-roro).^2).*(-yy);            % Componentes del campo vecorial x
fy=((1-roro).^2).*(xx);             % Componentes del campo vecorial y
quiver(xx,yy,fx,fy)                 % Dibuja el campo vectorial
axis([-2,2,-2,2])                   % Selecciona región del dibujo
view(2)                             % Muestra la figura

Podemos observar en el dibujo que los desplazamiento de las partículas del sólido es mayor cuanto mas cerca del centro del sólido que en el exterior, ya que ahí los vectores del campo tienen un módulo mayor, siendo los puntos del radio interior donde no hay desplazamiento.Aparte podemos observar como todos giran en el mismo sentido antihorario.

5 Sólido antes y después del desplazamiento

En este apartado vamos a representar nuestro sólido antes y después del desplazamiento dado por el campo de vectores anterior. Representamos ambas gráficas juntas para comparar el sólido anterior del desplazado.

Sólido antes y después desplazamineto

h=0.1;                              % Establecemos paso de muestra
ro=1:h:2;                           % Valor de ro según el muestreo del intervalo
teta=0:h:2*pi+h;                    % Valor de teta según el muestreo del intervalo
[roro,tetateta]=meshgrid(ro,teta);  % Matrices de la placa en coordenadas polares
figure(1)
xx=roro.*cos(tetateta);             % Parametrización
yy=roro.*sin(tetateta);             % Parametrización
subplot(1,3,1)                      % Figura solido momento inicial
mesh(xx,yy,0*xx)                  % Mallado
axis([-4,4,-4,4]);                  % Definimos la dimesión de los ejes 
view(2);                            % Visualización
Fx=xx+((1-roro).^2).*(-yy);         % Movimiento parametrizado 
Fy=yy+((1-roro).^2).*(xx);          % Movimiento parametrizado
subplot(1,3,2)                      % Figura solido en movimiento
mesh(Fx,Fy,yy.*0)                 % Mallado
axis([-4,4,-4,4]);                  % Definimos los ejes (selección de la región para dibujar)
view(2);                            % Visualización
subplot(1,3,3);                     % Figura solido en movimiento y campo vectorial 
hold on                             % Función que dibuja el campo vectorial y el solido en movimiento en una misma figura
mesh(Fx,Fy,yy.*0)                 % Mallado 
axis([-4,4,-4,4]);                  % Definimos los ejes
view(2);                            % Visualización
quiver(xx,yy,fx,fy)                 % Dibujamos el campo vectorial del gradiente de temperaturas
axis([-4,4,-4,4])                   % Definimos los ejes (selección de la región para dibujar)
view(2)                             % Visualización
hold off                            % Función que dibuja el campo vectorial y el sólido en movimiento en una misma figura

Podemos observar en las figuras como el solido se cambia de forma aumentando su area estirándose hacia el exterior, a la vez que la rotación lo deforma.

6 Divergencia

Aquí calculamos la divergencia del campo de vectores u en todos los puntos del sólido y la representamos. A través de esta gráfica veremos qué puntos tienen mayor divergencia. En la gráfica se debe observar que la divergencia es una medida del cambio de volumen local debido al desplazamiento.

                      [math]\bigtriangledown{}\cdot{}\overrightarrow{u}
=\frac{1}{\rho}[(\frac{\partial}{\partial\theta})(\rho\cdot{}(1-\rho)^2)] =0 [/math]

Este al salirnos nulo [math]\bigtriangledown{}\cdot{}\overrightarrow{u} =0[/math] para todo punto del dominio entonces el campo se le llama incompresible.

Divergencia del campo del desplazamineto

h=(0.1);                            % Paso de muestra
ro=1:h:2;                           % Valor de rho en el intervalo [1,2]
teta=0:h:2*pi+h;                    % Valor de theta en el intervalo [0,2*pi]
[roro,tetateta]=meshgrid(ro,teta);  % Matriz de la plaza en coordenadas polares
figure(1)
xx=roro.*cos(tetateta);             % Parametrizacion
yy=roro.*sin(tetateta);
f=0*yy;                             % Parametrizacion de la divergencia
surf(xx,yy,f)                       % Dibujo
axis([-2,2,-2,2])                   % Region en la que se situa el dibujo
colorbar
title('Divergencia del campo del desplazamiento')
view(2)                             % Visualizacion

7 Rotacional

En este apartado calculamos el rotacional del campo de desplazamientos u en todos los puntos del sólido y lo representamos.Este nos mostrara la tendencia del campo vectorial a inducir rotación alrededor de un punto. En nuestro caso quedara el rotacional:

               [math]\overrightarrow{\bigtriangledown{}}\times{}\overrightarrow{u}
=\frac{1}{\rho}\begin{vmatrix} \overrightarrow{g_\rho} & \overrightarrow{g_\theta} & \overrightarrow{g_z} \\ \frac{\partial}{\partial\rho} & \frac{\partial}{\partial\theta} & \frac{\partial}{\partial z} \\ 0 & \rho^2\cdot{}(1-\rho)^2 & 0 \end{vmatrix}=(2+4\cdot{}\rho^2-6\cdot{}\rho)\overrightarrow{g_z}[/math]

y su modulo quedando igual a [math]|\overrightarrow{\bigtriangledown{}}\times{}\overrightarrow{u}| =(2+4\cdot{}\rho^2-6\cdot{}\rho)[/math]

A través de la gráfica se observa cuáles son los puntos que sufren un mayor rotacional.

Rotacional del campo del desplazamineto

h=(0.1);                            % Paso de muestra
ro=1:h:2;                           % Valor de rho en el intervalo [1,2]
teta=0:h:2*pi+h;                    % Valor de theta en el intervalo [0,2*pi]
[roro,tetateta]=meshgrid(ro,teta);  % Malla de la placa en polares
figure(1)
xx=roro.*cos(tetateta);             % Parametrización
yy=roro.*sin(tetateta);
f=(2+4.*roro.^2-6.*roro);           % Función de la Rotación
surf(xx,yy,f)                       % Dibuja la malla
axis([-2,2,-2,2])                   % Región de del dibujo
colorbar
title('Rotacional del campo del desplazamiento')
view(2)                             % Visualización

En la imagen podemos observar que los puntos que tienen mayor rotacional son los más cercanos al radio mayor de la placa, siendo nulo en los puntos del radio menor.

8 Tensor de tensiones

Definimos la parte simétrica del tensor gradiente de u y calculamos el tensor de tensiones σ_ij a través de la fórmula: σ_ij = λ [math]\bigtriangledown{}\cdot{}\overrightarrow{u}[/math]δ_ij + 2 μ ε_ij. Tomamos λ=μ=1 y los vectores [math]\bar{i}[/math] = [math]\bar{g}[/math]_ρ y [math]\bar{j}[/math] = [math]\bar{g}[/math]_θ/ρ. Se nos pide dibujar las tensiones normales en la dirección que marca el eje [math]\bar{i}[/math], es decir [math]\bar{g}[/math]_ρ · σ · [math]\bar{g}[/math]_ρ y las tensiones normales en la dirección que marca el eje [math]\bar{j}[/math], es decir [math]\bar{g}[/math]_θ/ρ · σ · [math]\bar{g}[/math]_θ/ρ. Comparamos esta gráfica con las del módulo de la divergencia y el rotacional.

Comparacion Modulo Tensiones, Divergencia y Rotacional

h=0.1;                                   % Establecemos paso de muestra
ro=1:h:2;                                % Valor de ro según el muestreo del intervalo
teta=0:h:2*pi+h;                         % Valor de teta según el muestreo del intervalo
[roro,tetateta]=meshgrid(ro,teta);       % Matrices de la placa en coordenadas polares
subplot(1,3,1)
xx=roro.*cos(tetateta);                  % Parametrización
yy=roro.*sin(tetateta);  
f=0.*roro;                               % Modulo de las tensiones tangenciales
surf(xx,yy,f)                            % Dibuja la malla
axis                                     % Selecciona los limites de los ejes x e y en representación de la placa
colorbar                                 % Ponemos en la figura la leyenda de los colores
subplot(1,3,2) 
g=0.*roro;                               % Modulo de la divergencia
surf(xx,yy,g);          
subplot(1,3,3)
h=(2+4.*roro.^2-6.*roro);                % Modulo del rotacional
surf(xx,yy,h)

9 Tensiones Tangenciales respecto al plano ortogonal a gρ

Aquí dibujamos las tensiones tangenciales respecto al plano ortogonal gρ a través de la fórmula |σ · [math]\bar{g}[/math]_ρ - ([math]\bar{g}[/math]_ρ · σ · [math]\bar{g}[/math]_ρ)· [math]\bar{g}[/math]_ρ|. Indicamos donde son mayores y comparamos esta gráfica también con las del módulo de la divergencia y el rotacional.

Comparacion de las tensiones tangenciales al plano ortogonal gρ,Divergencia y rotacional

h=0.1;                                   % Establecemos paso de muestra
ro=1:h:2;                                % Valor de ro según el muestreo del intervalo
teta=0:h:2*pi+h;                         % Valor de teta según el muestreo del intervalo
[roro,tetateta]=meshgrid(ro,teta);       % Matrices de la placa en coordenadas polares
subplot(1,3,1)
xx=roro.*cos(tetateta);                  % Parametrización
yy=roro.*sin(tetateta);
f=roro.*(1-roro).*(-2+roro.*(1-roro).^2);% Modulo de las tensiones tangenciales
surf(xx,yy,f)                            % Dibuja la malla
axis                                     % Selecciona los limites de los ejes x e y en representación de la placa
colorbar                                 % Ponemos en la figura la leyenda de los colores
subplot(1,3,2)
g=0.*roro;                               % Modulo de la divergencia
surf(xx,yy,g);          
subplot(1,3,3)
h=(2+4.*roro.^2-6.*roro);                % Modulo del rotacional
surf(xx,yy,h)

10 Tensiones Tangenciales por el plano ortogonal a gθ/ρ

Representamos las tensiones tangenciales por el plano ortogonal a gθ/ρ a través de la fórmula |σ · [math]\bar{g}[/math]_ρ - ([math]\bar{g}[/math]_ρ · σ · [math]\bar{g}[/math]_ρ)· [math]\bar{g}[/math]_ρ|. Igual que en el apartado anterior, diremos dónde es mayor y las comparamos con el módulo de divergencia y rotacional.

Comparacion de las tensiones tangenciales al plano ortogonal gθ/ρ,Divergencia y rotacional

h=0.1;                                       % Establecemos paso de muestra
ro=1:h:2;                                    % Valor de ro según el muestreo del intervalo
teta=0:h:2*pi+h;                             % Valor de teta según el muestreo del intervalo
[roro,tetateta]=meshgrid(ro,teta);           % Matrices de la placa en coordenadas polares
subplot(1,3,1)
xx=roro.*cos(tetateta);                      % Parametrización
yy=roro.*sin(tetateta);
f=1./roro.*(1-roro).*(-2+roro.*(1-roro).^2); % Módulo de las tensiones tangenciales
surf(xx,yy,f)                                % Dibuja la malla
colorbar                                     % Ponemos en la figura la leyenda de los colores
subplot(1,3,2)
g=0.*roro;                                   % Módulo de la divergencia
surf(xx,yy,g);          
subplot(1,3,3)
h=(2+4.*roro.^2-6.*roro);                    % Módulo del rotacional
surf(xx,yy,h)

11 Tensión de Von Mises

Calculamos la tensión de Von Mises a partir de 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]. Esta es una magnitud escalar que pintaremos y analizándola indicamos en qué punto se alcanza el mayor valor.

Tension de Von Mises.

h=0.1;                                        % Establecemos paso de muestra
ro=1:h:2;                                     % Valor de ro según el muestreo del intervalo
teta=0:h:2*pi+h;                              % Valor de teta según el muestreo del intervalo
[roro,tetateta]=meshgrid(ro,teta);            % Matrices de la placa en coordenadas polares
subplot(1,1,1)
xx=roro.*cos(tetateta);                       % Parametrización
yy=roro.*sin(tetateta);
f=sqrt(3).*(1-roro).*(-2+roro.*(1-roro).^2);  % Modulo de las tensiones Von Mises
surf(xx,yy,f)                                 % Dibuja la malla
axis                                          % Selecciona los limites de los ejes x e y en representación de la placa
colorbar                                      % Ponemos en la figura la leyenda de los colores

Es mayor en los puntos de una circunferencia de radio rho=1.6, de valor máximo 1,5.

12 Masa total

En este apartado debemos hallar la masa. Para ello nos indican la función densidad [math]d(x,y,z)= xylog(x+2)[/math] de nuestra placa a partir de la cual podremos hallar la masa y representarla.

h=1/1000;                            % Establecemos paso de muestra
ro=1:h:2;                           
teta=0:h:2*pi+h;                     
[roro,tetateta]=meshgrid(ro,teta);   % Matrices de la placa en coordenadas polares
xx=roro.*cos(tetateta);              % coordenadas cartesianas
yy=roro.*sin(tetateta);         
d=xx.*yy.*log(xx+2);                 % Funcion de densidad d(x,y)
a=h^2*d;                             % Nos da la masa en cada punto
mas=sum(sum(a))                      % Este comando nos suma los elementos de la matriz a, la masa total

La masa total es 2.2097e-006

Explicación: Al crear la matriz de la placa en coordenadas cartesianas, hemos creado la matriz de la placa circular con una serie de puntos que separados por una distancia h (=Paso de muestra). Tomamos como dA el espacio comprendido entre 4 puntos que forman el cuadrado de menor dimensión posible, siendo h el lado de ese cuadrado. Por tanto, para calcular la masa en cada dA habrá que multiplicar ese área (=h^2) por la densidad. Posteriormente se suman todos los diferenciales de masa para obtener la masa total. Se puede observar que cuanto menor sea el paso de muestra, más exacta sera la aproximación de la masa ya que el dA será menor.