Visualización de Campos Escalares y Vectoriales en Elasticidad (grupo 14)

Trabajo realizado por estudiantes
Título	Campos escalares y vectoriales en elasticidad. Grupo 24
Asignatura	Teoría de Campos
Curso	2014-15
Autores	Paula de Santos Muñoz Ciro Rodriguez Matamoros Joaquín Sánchez Molina Íñigo Uraga Palacio Jorge Martín Sebastián
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

Para este análisis y representación de campos escalares en elasticidad consideramos una placa plana (en 2 dimensiones) que ocupa la región comprendida entre las parábolas :

• P1: 18y -81x²-1=0
• P2: 2y +x²-1=0

Para representarla se utilizará el siguiente sistema de coordenadas curvilíneas adaptado a la geometría que nos dan:

• x = uv
• y = 1/2(u²-v²)

Considerando que el dominio en el que estarán comprendidas u y v será:

centre

En ella vamos a suponer que tenemos definidas dos cantidades físicas. Por un lado la temperatura T(u,v), dependiente de las dos coordenadas curvilíneas (u,v), y por otro lado los desplazamientos [math]vec{u}[/math](x,y) producidos por la acción de una fuerza determinada. De esta forma, si definimos r₀(u,v) el vector de posición de los puntos de la placa antes de la deformación, la posición de cada punto (u,v) de la placa después de la deformación viene dada por:

vector posición

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

centre

donde [math]vec{a}[/math] y [math]vec{b}[/math] son vectores dados.

Además, en este trabajo supondremos lo siguiente:

centre

1 Situación inicial de la placa

Primero haremos un mallado para representar los puntos interiores del sólido, utilizando, un paso de muestreo h = 1/20 para las coordenadas u y v. El intervalo en el que representaremos comprende:

[math](x,y) \in [-1,1]*[-1,1][/math]

h=1/20;                             % Muestreo.
u=1/3:h:1;                          % Intervalo [1/3,1].
v=-1:h:1;                           % Intervalo [-1,1].
[uu,vv]=meshgrid(u,v);              % Malla.
xx=uu.*vv;                          % Parametrización X.
yy=(1/2).*((uu.^2)-(vv.^2));        % Parametrización Y.
plot(xx,yy);                        % Dibujo.
mesh(xx,yy,0*xx)                   
axis([-1,1,-1,1])                   
view(2)

Representación del sólido.

2 Líneas Coordenadas y Base Natural

2.1 Líneas Coordenadas

Las líneas coordenadas sirve para entender mejor la transformación a coordenadas curvilíneas, y se obtienen variando una de las coordenadas de la transformación u o v y manteniendo fija la restante. Hemos representado varias líneas coordenadas a base de dar un valor concreto a u o a v, dentro de sus respectivos intervalos:

xx11=uu.*0.5 ;xx12=uu.*-0.5;xx13=uu.*1;xx14=uu.*-1;
xx15=uu.*0.75;xx16=uu.*-0.75;xx17=uu.*0;                                % Parametrización X fijando v (0.5/-0.5/1/-1/0.75/-0.75/0) y cambiando u.
yy11=(1/2).*((uu.^2)-(0.5.^2));yy12=(1/2).*((uu.^2)-((-0.5).^2));
yy13=(1/2).*((uu.^2)-(1.^2));yy14=(1/2).*((uu.^2)-((-1).^2));
yy15=(1/2).*((uu.^2)-(0.75.^2));yy16=(1/2).*((uu.^2)-((-0.75).^2));     
yy17=(1/2).*((uu.^2)-(0.^2));                                           % Parametrización Y fijando v (0.5/-0.5/1/-1/0.75/-0.75/0) y cambiando u.
xx21=vv.*0.5;xx22=vv.*0.4;xx23=vv.*1;xx24=vv.*0.9;xx25=vv.*0.75;
xx26=vv.*0.65;xx27=vv.*(1/3);                                           % Parametrización X fijando u (0.5/0.4/1/0.9/0.75/0.65/0.333) y cambiando v.
yy21=(1/2).*((0.5.^2)-(vv.^2));yy22=(1/2).*((0.4.^2)-(vv.^2));
yy23=(1/2).*((1.^2)-(vv.^2));yy24=(1/2).*((0.9.^2)-(vv.^2));
yy25=(1/2).*((0.75.^2)-(vv.^2));yy26=(1/2).*((0.65.^2)-(vv.^2));
yy27=(1/2).*(((1/3).^2)-(vv.^2));                                       % Parametrización X fijando u (0.5/0.4/1/0.9/0.75/0.65/0.333) y cambiando v.
subplot(1,2,1);                                                         % Dibujaremos las líneas coordenadas en dos gráficas (u y v) en la misma imagen.
hold on                                                                 
mesh(xx11,yy11,0*xx);mesh(xx12,yy12,0*xx);mesh(xx13,yy13,0*xx);
mesh(xx14,yy14,0*xx);mesh(xx15,yy15,0*xx);mesh(xx16,yy16,0*xx);
mesh(xx17,yy17,0*xx);                                                   
axis([-1,1,-1,1])                                                       
view(2)                                                                
hold off                                                               
subplot(1,2,2);                                                         
hold on                                                                
mesh(xx21,yy21,0*xx);mesh(xx22,yy22,0*xx);mesh(xx23,yy23,0*xx);
mesh(xx24,yy24,0*xx);mesh(xx25,yy25,0*xx);mesh(xx26,yy26,0*xx);
mesh(xx27,yy27,0*xx);                                                  
axis([-1,1,-1,1])                                                      
view(2)                                                                 
hold off

A continuación se muestra las gráficas resultantes:

Líneas coordenadas fijando la variable v (izquierda) y u (derecha), respectivamente.

2.2 Base Natural

Al realizar una transformación a coordenadas curvilíneas, en nuestro caso de x e y a u y v, el vector de posición [math] \vec{r_o}[/math] se expresará así:

centre

La base natural [math]\vec{g_u}, \vec{g_v}[/math] tiene como vectores la derivada del vector posición [math] \vec{r_o}[/math] según las nuevas coordenadas u y v. Al tratarse de una placa plana (2 dimensiones solamente), sólo se requieren los vectores, ([math]\vec{g_u}, \vec{g_v}[/math]). Aún así, más adelante en el trabajo tendremos que considerar una tercera coordenada, por eso también incluiremos en nuestra base natural el vector [math]\vec{g_w}[/math].

• [math] \vec{g_u}=v\hat{i} +u \hat{j}[/math]
• [math] \vec{g_v}=u\hat{i} -v \hat{j}[/math]

plot(xx,yy);               % Dibujo.
hold on                    
mesh(xx,yy,0*xx)           
quiver(xx,yy,vv,uu);       % Representación del primer vector de la base natural en cada punto.
quiver(xx,yy,uu,-vv);      % Representación del segundo vector de la base natural en cada punto. 
axis([-1,1,-1,1])          % Región que dibujamos.
view(2)                    
hold off

Vectores de la base natural [math]\vec{g_u}[/math] y [math]\vec{g_v}[/math]

Sabemos que está bién porque los vectores que acabamos de representar son tangentes a las líneas coordenadas (Recordamos que [math]\vec{g_u}[/math] y [math]\vec{g_v}[/math] son las derivadas del vector de posición [math] \vec{r_0}[/math] respecto u y v).

3 Curvas de nivel por influencia de un foco de calor

La temperatura el producto de un foco de calor dado por el campo escalar

centre

T=exp(-(xx-yy).^2);                  % Función Temperatura.
subplot(1,2,1);                      
contour(xx,yy,T,20);                 % 20 líneas de nivel.
axis([-1,1,-1,1])                    
view(2)                             
subplot(1,2,2);                      
surf(xx,yy,T); colorbar;             
axis([-1,1,-1,1])                   
view(2)                             
max(max(T))                          % Valor máximo de la temperatura.

Gráficas de la variación de la temperatura en la placa.

Por la forma de la función temperatura, el valor máximo de temperatura en la placa va a ser 1, y se da cuando x es igual a y.

4 Gradiente de T

El gradiente de una función escalar es la dirección en la cual el campo crece más rápido. Nuestras curvas de nivel representan los puntos que tienen la misma temperatura. Así que nuestro vector gradiente de temperatura será siempre perpendicular a estas líneas de nivel, y lo obtendremos derivando la función temperatura respecto x e y respectivamente.

Función temperatura y sus respectivas derevidas parciales respecto x (izquierda) e y (derecha).

h=1/20;                                               % Muestreo.
u=1/3:h:1;                                            % Intervalo [1/3,1].
v=-1:h:1;                                             % Intervalo [-1,1].
[uu,vv]=meshgrid(u,v);                                
xx=uu.*vv;                                            
yy=(1/2).*((uu.^2)-(vv.^2));                          
f=exp(-(xx-yy).^2);                                   % Función Temperatura.
fx=(-2.*xx+2.*yy).*exp(-(xx-yy).^2);                  % Derivada con respecto a x de la función Temperatura.
fy=(2.*xx-2.*yy).*exp(-(xx-yy).^2);                   % Derivada con respecto a y de la función Temperatura.
hold on                                               
quiver(xx,yy,fx,fy)                                   % Representación de los vectores gradiente. 
contour(xx,yy,f,20);colorbar;                         
view(2)                                               
hold off

Vectores del gradiente superpuestos sobre las líneas de temperatura.

5 Campo de desplazamientos

Ahora aplicamos el campo de desplazamientos descrito en la introducción del artículo sobre la placa. Para ver como queda el sólido después de ésta aplicación habrá que hacer una serie de cálculos:

centre

centre

centre

centre

Todo esto llevado a Matlab para obtener las gráficas que nos piden se traduce en esto:

Ux=-2.*uu.*vv;                          % Componente x del campo de desplazamientos u.
Uy=-2.*uu.^2;                           % Componente y del campo de desplazamientos u.
subplot(1,3,1)                          % Muestra varias imágenes. 1ª Imagen.
quiver(xx,yy,Ux,Uy);                    % Representación del campo vectorial de desplazamientos u.
axis([-1,1,-1,1]);                      % Selecciona la región a dibujar.
xd=xx+Ux;                               % Componente x final del sólido deformado.
yd=yy+Uy;                               % Componente y final del sólido deformado.
subplot(1,3,2), mesh(xx,yy,0*xx);       % Muestra varias imágenes (2ª Imagen) y mallado completo.
axis([-1,1,-1,1])                       % Selecciona la región a dibujar.
view(2)                                 % Ver imagen desde arriba.
subplot(1,3,3), mesh(xd,yd,0*xd);       % Muestra varias imágenes (3ª Imagen) y mallado completo.
axis([-1,1,-2,1]);                      % Selecciona la región a dibujar.
view(2)                                 % Ver imagen desde arriba.

Las imágenes obtenidas fueron las siguientes:

Gráficas que muestran, de izquierda a derecha, el campo de vectores de desplazamiento, el mallado original de la placa, y el resultado final de la placa tras la aplicación de la fuerza.

6 Divergencia de un campo

La divergencia de un campo vectorial controla la diferencia entre el flujo saliente y el flujo entrante de un campo vectorial (sobre la superficie que rodea a un volumen de control), con lo cual, si la divergencia es positiva, el campo tiene "fuentes" , y si la divergencia es negativa se dice que tiene "sumideros" .

Pasamos a calcularla apoyándonos en cálculos previos.

centre

Como u^u es igual a 0, el cálculo se nos simplifica y finalmente queda:

centre

Divergencia=(-2.*(3*uu.^2+vv.^2)./(uu.^2+vv.^2));   % Divergencia.
surf(xx,yy,Divergencia);colorbar;                    
axis([-1,1,-1,1])

Imagen de la divergencia [math]\nabla\cdot\vec{u}[/math].

Podemos corroborar con la imagen de la placa deformada, que el mayor cambio de área se produce en la zona superior de la misma, la cual se deforma en sentido descendente; mientras que la menor variación se produce en los picos inferiores, los cuales se trasladan al lado opuesto, "doblándose" la placa sobre el eje de ordenadas.

7 Cálculo del rotacional de un campo vectorial

El rotacional o rotor es un operador vectorial que muestra la tendencia de un campo vectorial a inducir rotación alrededor de un punto.

Como dijimos previamente en el apartado de las bases naturales, para calcular el rotacional vamos a necesitar definir una tercera componenete [math]\vec{g_w}[/math]. Los cálculos serán los siguientes:

centre

centre

centre

8 Tensor de tensiones

En un medio elástico, lineal, isótropo y homogéneo, los desplazamientos permiten escribir el tensor de tensiones [math]\sigma ^ij=λ\nabla·\vec{u}1+2μ\epsilon[/math] siendo:

• [math]\epsilon (\vec{u})[/math] la parte simétrica del tensor gradiente de [math]\vec{u}[/math], [math]\nabla\vec{u}[/math]::

[math]\epsilon (\vec{u})=\frac{\nabla\vec{u}+ \nabla\vec{u}^t}{2}[/math]

• [math]λ[/math] y [math]μ[/math] los coeficientes de Lamé que dependen de las propiedades elásticas de cada material.

Se utilizarán estas expresiones para dibujar las tensiones normales en la dirección que marca [math]\vec{g_u}[/math] y la dirección que marca [math]\vec{g_v}[/math], es decir:

[math]\frac{\vec{g}_u}{|\vec{g_u}|}·σ·\frac{\vec{g}_u}{|\vec{g_u}|}\qquad\frac{\vec{g}_v}{|\vec{g_v}|}·σ·\frac{\vec{g}_v}{|\vec{g_v}|}[/math]

Bajo estas instrucciones se empieza a definir [math]\sigma ^ij[/math]. Para ello nos serviremos de la fórmula de las derivadas parciales covariantes [math]u^{i}, _j = \frac{\partial u^i}{\partial x^j}+\Gamma^i_{jk}u^k[/math]

1. Tomaremos [math]λ=μ=1[/math].

2. Se calcula el tensor gradiente de [math]\vec{u}[/math], [math]\nabla\vec{u}[/math]:

[math]\nabla\vec{u}= \begin{bmatrix} \frac{\partial u^u}{\partial u} + \Gamma ^1_{11} u^u & \frac{\partial u^u}{\partial v} + \Gamma ^1_{12} u^u & \frac{\partial u^u}{\partial w} + \Gamma ^1_{13} u^u \\ \frac{\partial u^v }{\partial u} + \Gamma ^{2}_{11}u^u & \frac{\partial u^v }{\partial v} + \Gamma ^{2}_{12}u^u & \frac{\partial u^v }{\partial w} + \Gamma ^{2}_{13}u^u \\ \frac{\partial u^w }{\partial u} + \Gamma ^{3}_{11}u^u &\frac{\partial u^w }{\partial v} + \Gamma ^{3}_{12}u^u&\frac{\partial u^w }{\partial w} + \Gamma ^{3}_{13}u^u\end{bmatrix}=\begin{bmatrix}-2- \frac{2u^2}{u^2+v^2} & \frac{-2uv}{u^2+v^2} & 0 \\\frac{2uv}{u^2+v^2} & \frac{-2u^2}{u^2+v^2} & 0\\0&0&0 \end{bmatrix}\qquad\nabla \vec{u}^t=\begin{bmatrix}-2- \frac{2u^2}{u^2+v^2} & \frac{2uv}{u^2+v^2} & 0 \\\frac{-2uv}{u^2+v^2} & \frac{-2u^2}{u^2+v^2} & 0\\0&0&0 \end{bmatrix}[/math]:

3. Parte simétrica del tensor gradiente [math]\nabla\vec{u}[/math], [math]\epsilon (\vec{u})[/math] :

Por lo tanto [math]\epsilon (\vec{u})=\begin{bmatrix}\frac{-4u^2-2v^2}{(u^2+v^2)} & 0 & 0 \\\ 0 & \frac{-2u^2}{(u^2+v^2)} & 0\\0&0&0 \end{bmatrix} [/math]

4. Cálculo de [math]σ[/math]: Se recuerda que la divergencia de [math]\vec{u}[/math] es [math]\nabla\cdot\vec u = -2\frac{3u^2+v^2}{u^2+v^2}[/math]

[math]\sigma ^ij=λ\nabla·\vec{u}1+2μЄ=\begin{bmatrix}\frac{-14u^2-6v^2}{u^2+v^2} & 0 & 0 \\\ 0 & \frac{-10u^2-2v^2}{u^2+v^2} & 0\\0 &0& \frac{-6u^2-2v^2}{u^2+v^2} \end{bmatrix}[/math]

Por último, para la representación de las tensiones normales en la dirección de [math]\vec{g_u}[/math] y la dirección de [math]\vec{g_v}[/math], se necesita definir:[math]\frac{\vec{g}_u}{|\vec{g_u}|}= \frac{v\vec {e_1}+ u\vec {e_2}}{\sqrt{v^2+u^2}}\qquad\frac{\vec{g}_v}{|\vec{g_v}|}=\frac{u\vec {e_1}-v\vec {e_2}}{\sqrt{v^2+u^2}}[/math]

Sabiendo esto, se puede proceder en Matlab al cálculo de las matrices que permitirán posteriormente la representación de las tensiones normales en la dirección que marca [math]\vec{g_u}[/math] y la dirección de [math]\vec{g_v}[/math]:

• [math]\frac{\vec{g}_u}{|\vec{g_u}|}·σ·\frac{\vec{g}_u}{|\vec{g_u}|}[/math]
• [math]\frac{\vec{g}_v}{|\vec{g_v}|}·σ·\frac{\vec{g}_v}{|\vec{g_v}|}[/math]

9 Tensión de Von Mises

La tensión de Von Miss es una magnitud escalar que se emplea como indicador de cuando un material inicia un comportamiento plástico (y no elástico puro). La tensión de Von Mises se define por la siguiente fórmula:

Tensión de Von Mises.

En la cual σ1, σ2 y σ3 son autovalores de σ.

sigma1=(-14*uu.^2-6*vv.^2)./(vv.^2+uu.^2);                                        % Autovalor 1.
sigma2=(-10*uu.^2-2*vv.^2)./(vv.^2+uu.^2);                                        % Autovalor 2.
sigma3=(-6*uu.^2-2*vv.^2)./(vv.^2+uu.^2);                                         % Autovalor 3.
Mises=sqrt(((sigma1-sigma2).^2+(sigma2-sigma3).^2+(sigma3-sigma1).^2)./2);        % Fórmula de Von Misses.
surf(xx,yy,Mises);                                                                % Visualización de superficie en 3D.
axis([-1,1,-1,1])                                                                 % Selecciona la región a dibujar.
max(max(Mises))                                                                   % Valor máximo de la tensión de Von Mises en toda la región.

Este fue el resultado:

Visualización en 3D de la placa tras la tensión de Von Mises.

El valor máximo de la tensión de Von Mises es 6.9282 y se alcanza en los puntos de la placa pertenecientes al eje de ordenadas.

10 Masa de la placa

Por último, hallaremos la masa de la placa utilizando la función densidad que nos dan:

Función densidad.

Para ello utilizaremos el Método del Trapecio en Matlab. Este método se basa en aplicar la función en cada punto y obtener una matriz de valores de la densidad en cada punto de la malla. Después esa matriz hay que multiplicarla por un vector fila y columna y sumando cada uno de los elementos de esa matriz obtendremos finalmente la masa total del sólido.

Como x e y pueden ser valores negativos, la función densidad podría ser negativa en algunos puntos de la placa. Para evitarlo, los resultados finales se obtendrán convirtiendo cada valor de la matriz de densidades en su valor absoluto, multiplicándolos más tarde por los pasos y sumándolos todos entre sí.

N1=200; N2=200;                          % 200 puntos.
a=1/3; b=1; c=-1; d=1;                   % Extremos de los intervalos.
h1=(b-a)/N1; h2=(d-c)/N2;                % Pasos.
u=a:h1:b; v=c:h2:d;                      % Intervalos.
[uu,vv]=meshgrid(u,v);                   % Malla.
xx=uu.*vv;                               % Parametrización X.
yy=(1/2).*((uu.^2)-(vv.^2));             % Parametrización Y.
d=(xx.^2+yy.^2)*log(1.+1./(xx.^2+1));    % Función Densidad.
D=abs(d);                                % Valor absoluto.
w1=ones(N1+1,1);                         
w1(1)=1/2; w1(N1+1)=1/2;                  
w2=ones(N2+1,1);                         
w2(1)=1/2; w2(N1+1)=1/2;                 
result=h1*h2*w2'*D*w1                    % Resultado.

El valor final de la masa obtenido es 34.1843.