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

De MateWiki
Saltar a: navegación, buscar
Trabajo realizado por estudiantes
Título Visualización de campos escalares y vectoriales en elasticidad
Asignatura Teoría de Campos
Curso 2014-15
Autores

Kevin Rincón Crespo

Alejandro Sans Jiménez

María Victoria Sesto Muñoz

José Manuel Vallejo Asín

Álvaro Villarino Redondo

Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

La finalidad del trabajo es el estudio a partir del uso de operadores diferenciales como la divergencia y el rotacional, cambios producidos por la temperatura [math]T(x,y,z)[/math], asi como los desplazamientos [math]\vec u(\rho,\theta)[/math], sobre una placa plana (en dimensión 2) que ocupa el anillo comprendido entre las circunferencias centradas en el origen de radios 1 y 2 .


1 Situación inicial

Se toma una placa plana que ocupa el anillo circular centrado en el origen y comprendido entre los radios 1 y 2. Para definir gráficamente dicho sólido, se diseñó el siguiente programa, haciendo uso de las coordenadas polares:

Mallado del anillo
h=0.1;                      %Paso del muestreo
u=1:h:2;                    %u [1,2]
v=0:h:2*pi+h;               %v [0,2*pi]
[RO,THETA]=meshgrid(u,v);   %Se generan las matrices de u y v
figure(1)
X=RO.*cos(THETA);           %Parametrizamos el anillo
Y=RO.*sin(THETA);
mesh(X,Y,0*X)               %Dibujamos la malla
axis([-3,3,-3,3])           %Definimos los ejes con su dominio correspondiente
axis equal
view(2)                     %Visualizamos el mallado


2 Comportamiento de la temperatura

La temperatura del sólido viene dada por la siguiente función: [math]T(x, y, x)=(8−y^2+2y)e^{−x^2+2x−1}[/math] El resultado obtenido es el gráfico que se puede observar a continuación: 


h=0.1;                                 %Paso del muestreo
u=1:h:2;                               %u [1,2]
v=0:h:2*pi+h;                          %v [0,2*pi]
[RO,THETA]=meshgrid(u,v);              %Se generan las matrices de u y v
figure(1)
X=RO.*cos(THETA);                     %Parametrizamos el anillo
Y=RO.*sin(THETA);
T=((8-Y.^2+2.*Y).*exp(-X.^2+2.*X-1)); %Definimos el campo escalar
subplot(1,2,1)
surf(X,Y,T)
axis([-3,3,-3,3]) ;                   %Definimos los ejes con su dominio correspondiente
axis equal
view(2)
subplot(1,2,2)
surf(X,Y,T)   
axis([-3,3,-3,3])
axis equal
view(3)
maximo=max(max(T))                            %Visualizamos el mallado


Temperatura máxima

El programa nos devuelve como máximo valor de la temperatura 8.9945 y este se sitúa en torno al punto (1,1). Sin embargo, este aspecto se puede tratar con más precisión estudiando el gradiente de temperaturas, que trataremos a continuación.

2.1 Variación de la temperatura

El gradiente de temperaturas indica la dirección de máxima variación de la misma en cada punto, es decir, en que dirección (y sentido) aumenta en mayor proporción. La expresión del gradiente queda como:

[math]∇T= \frac{∂T}{∂x} i ⃗+ \frac{∂T}{∂y}j ⃗+ \frac{∂T}{∂z} \vec k=(8-y^2+ 2y)(-2x+2) e^{-x^2+2x-1}\vec i+ (-2y+2)e^{-x^2+2x-1}j ⃗ [/math] Con la ayuda de MatLAB, podemos visualizar gráficamente el gradiente:

Gradiente de la temperatura
h=0.1;                                  % Paso del muestreo
u=1:h:2;                                % u [1,2]
v=0:h:2*pi+h;                           % v [0,2*pi]
[RO,THETA]=meshgrid(u,v);               % Se generan las matrices de u y v
figure(1)
X=RO.*cos(THETA);                       % Parametrizamos el anillo
Y=RO.*sin(THETA);
T=((8-Y.^2+2.*Y).*exp(-X.^2+2.*X-1));   % Definimos el campo escalar
contour(X,Y,T)                          % Dibujamos las líneas de nivel
hold on
Tx=(8-Y.^2+2.*Y).*(-2.*X+2).*exp(-X.^2+2.*X-1);     % componente x del campo vectorial 
Ty=(-2.*Y+2).*exp(-X.^2+2.*X-1);        % componente y del campo vectorial
quiver(X,Y,Tx,Ty)                       % Dibujamos el campo vectorial
axis([-3,3,-3,3])                       % Definimos los ejes con su dominio correspondiente
view(2)                                 % Vista en planta
colorbar


Tal y como vemos en la gráfica, los vectores gradiente asociados a cada punto del sólido tienen una dirección ortogonal a las líneas equipotenciales, es decir, aquellas que unen puntos con igual temperatura. Por otro lado, en las proximidades de la zona más caliente del anillo convergen los vectores, ya que el gradiente se anula. Esto se debe a que en el máximo, el valor de la temperatura no crece.

3 Deformaciones

3.1 Campo de desplazamiento [math]u ⃗ [/math]

[math]u ⃗= (1- ρ)^2 g ⃗_ρ[/math] [math]Recordemos que: g ⃗_ρ=cos⁡(θ)•i ⃗+sin⁡(θ)•j ⃗+0•k ⃗ [/math]

Campo de desplazamiento[math]u ⃗ [/math]
h=0.1;                  % paso de muestreo
u=1:h:2;               % u [1,2]
v=0:h:2*pi+h;            % v [0,2*pi]
[RO,THETA]=meshgrid(u,v); % matrices de las coorenadas RO y THETA
figure(1)
X=RO.*cos(THETA);        % parametrización
Y=RO.*sin(THETA);    
hold on
mesh(X,Y,0*X)               %Dibujamos la malla
axis([-3,3,-3,3]) %Definimos los ejes con su dominio correspondiente
axis equal
view(2)           % Vista en planta
fx=((1-RO).^2.*cos(THETA));                 
fy=((1-RO).^2.*sin(THETA));  
quiver(X,Y,fx,fy)     % Dibujamos el campo vectorial
axis([-3,3,-3,3])      % Seleccionamos el dominio del dibujo
axis equal
view(2)
hold off

El campo [math]u ⃗ [/math] actúa en dirección radial, ya que sólo tiene componentes según [math]g ⃗_ρ[/math]. Este comportamiento se estudiará posteriormente con la divergencia y el rotacional. Si utilizamos el comando subplot, se puede comparar la figura antes y después de aplicarle el campo[math]u ⃗ [/math]:

Efectos del campo [math]u ⃗ [/math] sobre el mallado
h=0.1;                      %Paso del muestreo
u=1:h:2;                    %u [1,2]
v=0:h:2*pi+h;               %v [0,2*pi]
[RO,THETA]=meshgrid(u,v);   %Se generan las matrices de u y v
figure(1)
X=RO.*cos(THETA);           %Parametrizamos el anillo
Y=RO.*sin(THETA);
subplot(1,3,1);
mesh(X,Y,0*X)               %Dibujamos la malla del anillo sin deformar
axis([-3,3,-3,3])           %Definimos los ejes con su dominio correspondiente
axis equal
view(2) 
subplot(1,3,2);             %se abre el subplot
fx=(RO.^2-RO+1).*cos(THETA);%tras calcular a mano el r desplazado, se introduce
fy=(RO.^2-RO+1).*sin(THETA);
mesh(fx,fy,0*fx);           % Dibujo de malla desplazada
axis([-3,3,-3,3])           % Definimos los ejes con su dominio correspondiente
axis equal
view(2);
subplot(1,3,3);
hold on                     %abrimos hold on para superponer las dos graficos siguientes
mesh(fx,fy,0*fx)            %Dibujamos la malla deformada
axis([-3,3,-3,3])           %Definimos los ejes con su dominio correspondiente
axis equal
view(2)
quiver(X,Y,fx,fy)           % Dibujamos el campo vectorial
axis([-3,3,-3,3])           % Definimos los ejes con su dominio correspondiente
axis equal
view(2)
hold off


3.2 Divergencia

La divergencia indica la variación de volumen local en el sólido debido al desplazamiento. Ésta viene dada por: [math]∇• u ⃗= \frac{1}{√g} \frac{∂}{∂x^i} (√g u^i )=\frac{1}{p}\frac{∂}{∂ρ}[ρ(1-ρ)^2]= \frac{1}{ρ}+3ρ-4[/math]

Divergencia
h=0.1;                       %Paso del muestreo
u=1:h:2;                    %u [1,2]
v=0:h:2*pi+h;               %v [0,2*pi]
[RO,THETA]=meshgrid(u,v);   %Se generan las matrices de u y v
figure(1)
X=RO.*cos(THETA);           %Parametrizamos el anillo
Y=RO.*sin(THETA);
mesh(X,Y,0*X)               %Dibujamos la malla
f=(1./RO)-4+3.*RO;% Módulo de las tensiones en la dirección radial
subplot(1,2,1)
surf(X,Y,f)
axis equal %Definimos los ejes con su dominio correspondiente
view(3)
subplot(1,2,2)
surf(X,Y,f)
axis equal %Definimos los ejes con su dominio correspondiente
view(2)

Cuando el valor de ρ es igual a 1, la divergencia de [math]u ⃗ [/math]se anula (y también el campo), tal y como se puede ver en ambos gráficos. Por otro lado, cuando ρ es igual a 2, la divergencia adquiere su valor máximo. Esto indica que el campo de deformación tiene más incidencia en la zona exterior del anillo que no en la interior.

3.3 Rotacional

El rotacional indica la cantidad de giro sobre una superficie y su expresión es la siguiente:

[math] \nabla \times \vec u=\begin {bmatrix} \vec g_{\rho} & \vec g_{\theta} & \vec g_{z} \\ \frac{∂}{∂ρ} & \frac{∂}{∂\theta} & \frac{∂}{∂z} \\ \ (1- ρ)^2 & 0 & 0 \end{bmatrix}=\vec {0}\ [/math]

Esto es, el campo [math]u ⃗ [/math] es irrotacional. La interpretación que podemos dar es que no se produce giro entorno al vector normal de la superficie. Y es que es lógico, ya que, como indicamos anteriormente, el campo [math]u ⃗ [/math] viene expresado según [math]g ⃗_p[/math].

4 Tensiones

4.1 Tensor de tensiones σ

[math]σ= λ∇•u ⃗ 1+2μϵ[/math]

Siendo λ=μ=1 y;


[math]ϵ(u ⃗ )= \frac{∇u ⃗+ ∇u ⃗^t }{2} [/math]

[math]ϵ(u ⃗ )=(2ρ-2) g ⃗_p⨂g ⃗_p+ (\frac{1}{ρ}- \frac{2}{ρ^2} + \frac{1}{ρ^3} ) g ⃗_ϑ⨂g ⃗_ϑ[/math]

[math]{σ_{ij}} =(\frac{1}{ρ}+ 3ρ-4)''1''+(4ρ-4) g ⃗_p⨂g ⃗_p+ (\frac{2}{ρ}- \frac{4}{ρ^2} + \frac{2}{ρ^3} ) g ⃗_ϑ⨂g ⃗_ϑ [/math]

Pasamos a una base ortonormal tal que [math]{e ⃗_1,e ⃗_2,e ⃗_3} = {g ⃗_p,g ⃗_ϑ/p,g ⃗_z}[/math]; operamos poniendo que[math] ''1'' =e ⃗_1⨂e ⃗_1,+e ⃗_2⨂e ⃗_2+e ⃗_3⨂e ⃗_3[/math] y obtenemos la siguiente matriz [math]σ_ij[/math]:

[math]{σ_ij}=\begin{bmatrix} 7ρ-8+ \frac{1}{ρ} & 0 & 0 \\ 0 & 5ρ-8+ \frac{3}{ρ} & 0 \\ 0 & \ 0 & 3ρ-4+ \frac{1}{ρ} \end{bmatrix} [/math]:

4.2 Tensiones normales en la dirección de [math]g ⃗_ρ[/math]

Tensiones normales en la dirección de [math]g ⃗_p→ g ⃗_p• σ• g ⃗_p= 7ρ-8+ \frac{1}{ρ} [/math]

Tensiones normales en la dirección de [math]g ⃗_p→ g ⃗_p• σ• g ⃗_p[/math]
h=0.1;                      %Paso del muestreo
u=1:h:2;                    %u [1,2]
v=0:h:2*pi+h;               %v [0,2*pi]
[RO,THETA]=meshgrid(u,v);   %Se generan las matrices de u y v
figure(1)
X=RO.*cos(THETA);           %Parametrizamos el anillo
Y=RO.*sin(THETA);
hold on
mesh(X,Y,0*X)               %Dibujamos la malla
axis([-3,3,-3,3])           %Definimos los ejes con su dominio correspondiente
axis equal
f=(1./RO)+7*RO-8;           % Módulo de las tensiones en la dirección radial
surf(X,Y,f)
axis([-3,3,-3,3])
view(3)
hold off


Las tensiones normales, como ya habíamos visto en la divergencia, son mayores a medida que nos acercamos a la zona exterior del anillo.

4.3 Tensiones normales en la dirección de [math]g ⃗_ϑ/ρ[/math]

Tensiones normales en la dirección de:: [math]\frac{g ⃗_ϑ}{ρ}→ \frac{g ⃗_ϑ}{ρ}•σ• \frac{ g ⃗_ϑ}{ρ}= 5ρ-8+ \frac{3}{ρ} [/math]

Tensiones normales en la dirección de [math]\frac{g ⃗_ϑ}{ρ}→ \frac{g ⃗_ϑ}{ρ}•σ• \frac{ g ⃗_ϑ}{ρ}= 5ρ-8+ \frac{3}{ρ}[/math]
h=0.1;                       %Paso del muestreo
u=1:h:2;                    %u [1,2]
v=0:h:2*pi+h;               %v [0,2*pi]
[RO,THETA]=meshgrid(u,v);   %Se generan las matrices de u y v
figure(1)
X=RO.*cos(THETA);           %Parametrizamos el anillo
Y=RO.*sin(THETA);
hold on
mesh(X,Y,0*X)               %Dibujamos la malla
axis([-3,3,-3,3]) %Definimos los ejes con su dominio correspondiente
axis equal
f=(3./RO)+5.*RO-8;% Módulo de las tensiones en la dirección g_theta
surf(X,Y,f)
axis([-3,3,-3,3])
view(3)
hold off


Las gráficas se parecen a la de la divergencia, son troncos de cono. En el caso de [math]g ⃗_p• σ• g ⃗_p[/math] se trata de un cono invertido, los valores de x e y varían entre -2 y 2 y los de z entre 0 y 8, mientras que en [math]\frac{g ⃗_ϑ}{ρ}•σ•\frac{g ⃗_ϑ}{ρ}[/math] x e y varían entre -2 y 2 y z entre 0 y 3.5.

4.4 Tangenciales respecto al plano ortogonal a [math]g ⃗_ρ[/math]

[math]|σ• g ⃗_p- (g ⃗_p• σ• g ⃗_p)g ⃗_p |= 0[/math]

Por lo que no hay tensiones en el plano ortogonal a [math]g ⃗_p.[/math]

4.5 Tangenciales respecto al plano ortogonal a [math]g ⃗_ϑ/ρ[/math]

[math]|σ• \frac{g ⃗_ϑ}{ρ}-(\frac{g ⃗_ϑ}{ρ}• σ•\frac{g ⃗_ϑ}{ρ})\frac{g ⃗_ϑ}{ρ}|=0 [/math]

Por lo que no hay tensiones en el plano ortogonal a [math]\frac{g ⃗_ϑ}{ρ}[/math]

4.6 Tensión de Von Mises

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

Para obtenerlo, en nuestro caso concreto, se ha diseñado el programa que se muestra a continuación. Éste está basado en obtener un vector de dimensión 1 x n que contenga los valores de la tensión de Von Mises para cada valor de ρ, dependiendo n del paso de muestreo que se le dé a ρ. Hemos usado h=0.01: 

A=zeros(3,101);                     %se define una matriz vacía donde irán las triadas de autovalores para cada valor de RO
k=0;                                %se inicia un contador
tensvonmises=zeros(1,101); %se define un vector vacío que contendrá los valores de la tension de Von Mises para cada valor de RO
for i=1:0.01:2                      %paso de muestreo 0.01
    RO=i;
    sigmaij=[7.*RO-8+1./RO, 0, 0; 0, 5.*RO-8+3./RO, 0; 0, 0, 1./RO+3.*RO-4]; %matriz de tensor de tensiones sigma
    v=eig(sigmaij);                 %se obtienen los autovalores
    k=k+1;                          %el contador aumenta 1 unidad en cada entrada al bucle
    A(1:3,k)=v;                     %se meten las triadas de autovalores en las columnas de A
    tensvonmises(k)=sqrt(((v(1)-v(2)).^2+(v(2)-v(3)).^2+(v(3)-v(1)).^2)/2); 
    end                          %se van guardando los valores de von mises en el vector vacío


Una vez llegados a este punto, se trata de obtener una función [math]σ_{VM}=σ_{VM}(ρ)[/math] que nos permita pintar la superficie para la tensión de Von Mises, de la misma manera que hemos hecho previamente para otros elementos como el módulo de la divergencia, por ejemplo. Hemos utilizado un comando de Matlab llamado polyfit(x,y,n), que interpola un conjunto de puntos (x,y) mediante un polinomio de grado n, devolviendo un vector con los coeficientes de dicho polinomio en orden descendente de grado. Se trata de obtener el polinomio que pase por todos nuestros puntos formados por los vectores RO=1:0.01:2 y tensvonmises que calculó el programa anterior: 

RO=1:0.01:2;
p=polyfit(RO,tensvonmises,6); %se interpola una función de RO que pase por los puntos de la tensión de von mises, 
%es decir, se saca una expresión de sigma=f(RO). Esto nos da los coeficientes de un polinomio de grado n 
%(con grado 6 ya se obtiene uno aceptable en nuestro caso) en forma de vector.


Finalmente, pintamos la superficie de la tensión de Von Mises igual que siempre 

h=0.1;                      %Paso del muestreo
u=1:h:2;                    %u [1,2]
v=0:h:2*pi+h;               %v [0,2*pi]
[RO,THETA]=meshgrid(u,v);   %Se generan las matrices de u y v
figure(1)
X=RO.*cos(THETA);           %Parametrizamos el anillo
Y=RO.*sin(THETA);
mesh(X,Y,0*X)               %Dibujamos la malla
f=p(1).*(RO.^6)+p(2).*(RO.^5)+p(3).*(RO.^4)+p(4).*(RO.^3)+p(5).*(RO.^2)+p(6).*(RO)+p(7); %Utilizamos la funcion sigma=f(RO) antes interpolada
subplot(1,2,1)
surf(X,Y,f)
axis equal                  %Definimos los ejes con su dominio correspondiente
view(3)
subplot(1,2,2)
surf(X,Y,f)
axis equal                  %Definimos los ejes con su dominio correspondiente
view(2)


Por tanto, el resultado final de la tensión de Von Mises queda como:

Tensión de Von Mises

Debido a que todos los gráficos son muy parecidos, creamos un programa para verlos todos juntos y poder compararlos mejor: 

h=0.1;                       %Paso del muestreo
u=1:h:2;                    %u [1,2]
v=0:h:2*pi+h;               %v [0,2*pi]
[RO,THETA]=meshgrid(u,v);   %Se generan las matrices de u y v
figure(1)
X=RO.*cos(THETA);           %Parametrizamos el anillo
Y=RO.*sin(THETA);
mesh(X,Y,0*X)               %Dibujamos la malla
f1=(1./RO)-4+3.*RO;% Módulo de las tensiones en la dirección radial
subplot(1,3,1)
surf(X,Y,f1)
axis equal %Definimos los ejes con su dominio correspondiente
view(3)
title('DIVERGENCIA')
subplot(1,3,2)
f2=(1./RO)+7*RO-8;% Módulo de las tensiones en la dirección radial
surf(X,Y,f2)
axis([-3,3,-3,3])
view(3)
title('TENSION NORMAL RO')
subplot(1,3,3)
f3=(3./RO)+5.*RO-8;% Módulo de las tensiones en la dirección radial
surf(X,Y,f3)
axis([-3,3,-3,3])
view(3)
title('TENSION NORMAL THETA/RO')
Obtenemos lo siguiente:
Comparación de la divergencia y tensiones normales

5 Masa

Sabiendo que la densidad de la placa viene dada por la función:: [math] d(x,y,z)=xye^{\frac{-1}{x^2}}[/math]


[math]\int_{S}f dS = \int_{}\int_{[uxv]=[1,2]x[0,2pi]}f(u,v) [/math]


[math] =\int_{}\int_{[uxv]=[1,2]x[0,2pi]}\frac{1}{2}u^2 sin(2v)e^{-\frac{1}{u^2cos(v)}}|g_ρ×g_ϑ |dudv[/math]


[math]= \int_{}\int_{[uxv]=[1,2]x[0,2pi]}\frac{1}{2}u{^\frac{3}{2}} sin(2v)e^{-\frac{1}{u^2cos(v)}}dudv [/math] 


h=0.01;
u=1:h:2;
v=0:h:pi/2;
[U,V]=meshgrid(u,v);
f=(1/2).*U.^(3/2).*sin(2.*V).*exp(-1./(U.^2.*cos(V).^2));
a=h.^2.*f;
m=4*sum(sum(a));


Obtenemos una masa de 1.4345 


[[Categoría:Teoría de Campos]]
[[Categoría:TC14/15]]
Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Visualización_de_campos_escalares_y_vectoriales_en_elasticidad_(Grupo_3D)&oldid=23069»