Visualización de campos escalares y vectoriales en elasticidad (Grupo A1)
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Visualización de campos escalares y vectoriales en elasticidad |
| Asignatura | Teoría de Campos |
| Curso | 2014-15 |
| Autores | Inés Arana Bachelis, Javier Blanco Villarroel, Paula Botella Andreu, Marta Cavero Guillén |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
Contenido
1 Introducción
Consideramos una placa plana (en dimensión 2) que ocupa la región comprendida entre las parábolas P1: 18y-81x^2-1=0, y P2: 2y+x^2-1=0. Para representarla usaremos un sistema de coordenadas adaptado a su geometría:
- x=u*v
- y=0.5(u^2-v^2)
con u, v definidas en (u,v)∈ [1/3,1] X [-1,1].
En ella vamos a suponer que tenemos definidas dos cantidades físicas: la temperatura T(u,v), que depende de las dos coordenadas curvilíneas (u,v) y los desplazamientos vec(u)(x,y) producidos por la acción de una fuerza determinada. De esta forma, si definimos r0(u,v) el vetor 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: vec(r)(u,v)=vec(r0)(u,v)+vec(u)(u,v).
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: vec(u)(u,v)=vec(a)(vec(b)*vec(r0)), donde vec(a)=vec(gu)/mod(vec(gu)) y vec(b)=-4*vec(gu)/mod(vec(gu)).
Para la realización de este trabajo utilizaremos MATLAB como herramienta principal.
2 Dibujo de la gráfica
Para el dibujo de la gráfica dibujaremos un mallado con equidistancia 1/20 para las coordenadas u y v que represente los puntos interiores del sólido, todo esto, tomando como ejes x e y en el intervalo [-1,1] X [-1,1].
El programa utilizado ha sido:
h=1/20;
u=1/3:h:1;
v=-1:h:1;
[Mu,Mv]=meshgrid(u,v);
x=Mu.*Mv;
y=0.5.*((Mu.^2)-(Mv.^2));
z=zeros(size(Mu));
mesh(x,y,z);
axis([-1,1,-1,1]);
view(2)
Y hemos obtenido la siguiente gráfica:
3 Lineas coordenadas y vectores de la base natural
Para la línea coordenada de v:
- x=k*v
- y=0.5*(k^2-v^2)/-> x/k=(k^2-2*y)^(-1/2)1
Operando las dos ecuaciones llegamos a una función para las curvas coordenadas que define una familia de curvas dependiente de la constante k:
- y=(k^2/2)-(x^2/2k^2)
Para la curva coordenada de u consideramos la v como constante de la misma manera:
- x=u*k
- y=0.5*(k^2-v^2)
Tras operar se llega a la ecuación para la familia de curvas:
- y=(x^2/2k^2)-k^2
Para visualizar los vectores de la base natural hemos utilizado el siguiente programa:
h=1/20;
u=1/3:h:1;
v=-1:h:1;
[Mu,Mv]=meshgrid(u,v);
x=Mu.*Mv;
y=0.5.*((Mu.^2)-(Mv.^2));
subplot(1,2,1)
U=y;
V=x;
quiver(x,y,U,V);
axis([-1,1,-1,1]);
title('Vec(gu) de la base natural');
subplot(1,2,2)
A=x;
B=-y;
quiver(x,y,A,B);
axis([-1,1,-1,1]);
title('Vec(gv) de la base natural');
Para visualizar las líneas coordenadas hemos implementado el siguiente programa:
u=1/3:1/20:1;
v=-1:1/20:1;
[Mu,Mv]=meshgrid(u,v);
x=Mu.*Mv;
y=0.5.*((Mu.^2)-(Mv.^2));
z=zeros(size(Mu));
subplot(1,2,1)
mesh(x,y,z)
axis([-1,1,-1,1]);
view(2)
subplot(1,2,2)
hold on
a=-1:0.05:1;
b=(0.125-((a.^2)/0.5));
plot(a,b,'r')
c=-1:0.05:1;
d=(c.^2/0.5)-0.125;
plot(c,d,'b')
hold off
4 Campo escalar de la temperatura
Dado el campo escalar función de x e y, se pide representar las curvas de nivel sobre el mallado de la placa y determinar el punto de mayor temperatura, así como representar el gradiente del mismo y observar gráficamente que éstas son ortogonales a las curvas de nivel.
4.1 Curvas de nivel
Para representarlas se implementó el siguiente programa:
:h=1/20;
:u= 1/3:h:1;
:v= -1:h:1;
:%Mallado
:[uu,vv]= meshgrid(u,v);
:%Parametrización
:x=uu.*vv;
:y=1/2*(uu.^2-vv.^2);
:hold on
:plot(x,y,'k')
:%Definición de función f=Temperatura
:f=inline('(8-y.^2+2.*y).*(e.^(-(x.^2)))','x','y');
:%Curvas de nivel
:z=f(x,y);
:contour(x,y,z)
:hold off
4.2 Gradiente del campo de temperatura
Definimos el concepto matemático de gradiente como la diferencia del valor del campo entre dos puntos infinitamente próximos mediante un campo escalar, en este caso, nuestro campo escalar es la temperatura,que sería máxima, si lleva la dirección del gradiente, y mínima en caso contrario. Así la expresión matemática a la que nos referimos es:
> Campo escalar de Temperaturas: T=T(x,y)
> Gradiente de Temperaturas: ∇T=dT(x,y)/dx+dT(x,y)/dy
A continuación representamos dicho resultado en Matlab con el siguiente programa, teniendo en cuenta que las coordenadas (x,y) van dadas en función de las coordenadas parabólicas (u,v), y que la temperatura viene dada por la expresión del apartado anterior:
h=1/20;
u=1/3:h:1;
v=-1:h:1;
%Mallado
[Mu,Mv]=meshgrid(u,v);
%Parametrización
x= Mu.*Mv;
y= 1/2*(Mu.^2-Mv.^2);
%Definición de función f=Temperatura.
f= inline('(8-y.^2+2.*y)*(e.^(-(x.^2)))','x','y');
%Derivas parciales de la Temperatura respecto a x--->fx y respecto a y--->fy
fx= inline('(-2.*x.*e.^(-(x.^2))).*(8-y.^2+2.*y)','x','y');
fy= inline('(-2.*y+2).*e.^(-(x.^2))','x','y');
%Definición de las funciones mediante matrices para su visualización.
T=(8-y.^2+2.*y).*(exp(-(x.^2)));
U=(-2.*x.*exp(-(x.^2))).*(8-y.^2+2.*y);
V=(-2.*y+2).*exp(-(x.^2));
%Creación de gráficas.
hold on
contour(x,y,T)
plot(Mu.*(-1), 1/2.*(Mu.^2-(-1).^2),'b')
plot(Mu.*(1), 1/2.*(Mu.^2-(1).^2),'b')
plot((1/3).*Mv, 1/2.*((1/3).^2-Mv.^2),'b')
plot((1).*Mv, 1/2.*((1).^2-Mv.^2),'b')
axis([-1,1,-1,1])
hold on
quiver(x,y,U,V); colorbar
hold off
axis([-1,1,-1,1]);
title('Gradiente de temperaturas')
Se puede observar gráficamente, como el gradiente de temperatura (debido a la propia definición de gradiente) es ortogonal a las curvas de nivel en cada punto (uo,vo), y a su vez, cómo el campo aumenta en módulo en la dirección de dicho gradiente.
5 Campo de desplazamientos
Dado el campo de desplazamiento, u(u, v) =a(b · r0) en función de las coordenadas parabólicas (u,v), siendo: a=gu/|gu|
b = −4gu/|gu|
r0= vector de posición de los puntos de la placa antes de la deformación.
Estudiaremos la representación gráfica del campo de vectores, veremos qué desplazamiento produce este campo en el sólido,así mismo, calcularemos y representaremos su divergencia y su gradiente.
5.1 Representación del campo de vectores
Para representar el campo de desplazamiento en los puntos del mallado del sólido hemos utilizado el siguiente programa:
h=0.05;
u=1/3:h:1;
v=-1:h:1;
%Mallado
[Mu,Mv]= meshgrid(u,v);
%Parametrización
x= Mu.*Mv;
y= 1/2*(Mu.^2-Mv.^2);
%Definición de la función desplazamiento en x e y
gx= -2.*Mu.*Mv;
gy= -2.*Mu.^2;
%Mallado de flechas
quiver(x,y,gx,gy);
axis([-1,1,-1,1]);
%Se observa cómo el vector desplazamiento lleva la dirección tangente a las curvas de nivel.
%Se diferencia claramente el sentido del vector desplazamiento hacia el centro de la superficie.
title('Campo de desplazamiento')
Se observa en la visualización la dirección y sentido en la que se desplaza cada punto del sólido, mediante el vector correspondiente aplicado en dicho punto.
5.2 Desplazamiento producido en el sólido
Para observar el efecto que produce el campo de desplazamiento sobre el sólido dibujaremos el sólido antes y después del desplazamiento dado por el campo de vectores vec(u). Para ello hemos utilizado el siguiente programa:
h=1/20;
u=1/3:h:1;
v=-1:h:1;
[Mu,Mv]=meshgrid(u,v);
subplot(2,2,1)
x=Mu.*Mv;
y=0.5*((Mu.^2)-(Mv.^2));
z=zeros(size(Mu));
mesh(x,y,z);
axis([-1,1,-1,1]);
view(2)
title('Sólido Original');
subplot(2,2,2)
x1=(-2*Mu.*Mv)+x;
y1=(-2*(Mu.^2))+y;
mesh(x1,y1,z);
axis([-1,1,-2.5,0.5]);
view(2)
title('Sólido Desplazado');
subplot(2,2,3)
x=Mu.*Mv;
y=0.5*((Mu.^2)-(Mv.^2));
z=zeros(size(Mu));
mesh(x,y,z);
axis([-1,1,-1,1]);
view(3)
title('Sólido Original 3D');
subplot(2,2,4)
x1=(-2*Mu.*Mv)+x;
y1=(-2*(Mu.^2))+y;
mesh(x1,y1,z);
axis([-1,1,-2.5,0.5]);
view(3)
title('Sólido Desplazado 3D');
5.3 Divergencia del campo vectorial
Se define la divergencia del campo vectorial vec(u) como la traza de su gradiente, siendo la divergencia de vec(U) : div(v)=-2v. Una vez calculada, procedemos a dibujarla en todos los puntos del sólido. En la gráfica observamos en 2 dimensiones que la zona con mayor divergencia aparece en color rojo, mientras que en la gráfica en 3 dimensiones observamos que esta zona se corresponde con la más elevada, es decir, la que se ha expandido en mayor proporción, que al final es la interpretación de la divergencia. Para obtener dichas gráficas hemos utilizado el siguiente programa:
h=1/20;
u=1/3:h:1;
v=-1:h:1;
[Mu,Mv]=meshgrid(u,v);
x=Mu.*Mv;
y=0.5.*((Mu.^2)-(Mv.^2));
div=-2*Mv; %Campo escalar resultado del módulo del rotacional.
subplot(1,2,1)
surf(x,y,div);
axis([-1,1,-1,1]);
view(2)
title('Divergencia vec(u) 2D');
subplot(1,2,2)
mesh(x,y,div);
axis([-1,1,-1,1]);
view(3)
title('Divergencia vec(u) 3D');
5.4 Rotacional del vector desplazamiento
El módulo del rotacional de un campo vectorial es una medida de “cuanto gira” el campo, siendo su dirección a del eje en torno al cual gira dicho campo.
Tras realizar los cálculos a mano, aplicando la fórmula del rotacional, obtenemos que el módulo del vector rotacional mod(Vxvec(u))=2u.
Mediante el siguiente programa MATLAB, dibujamos el módulo:
h=1/20;
u=1/3:h:1;
v=-1:h:1;
[Mu,Mv]=meshgrid(u,v);
x=Mu.*Mv;
y=0.5.*((Mu.^2)-(Mv.^2));
Vx=2*Mu; %Campo escalar resultado del módulo del rotacional.
subplot(1,2,1);
surf(x,y,Vx);
axis([-1,1,-1,1]);
view(2)
title('Módulo rotacional vec(u) 2D');
subplot(1,2,2);
mesh(x,y,Vx);
axis([-1,1,-1,1]);
view(3)
title('Módulo rotacional vec(u) 3D');
Y obtenemos la siguiente gráfica:
En la gráfica observamos en 2 dimensiones que la zona con mayor módulo del rotacional aparece en color rojo, mientras que en la gráfica en 3 dimensiones observamos que esta zona se corresponde con la más elevada, es decir, se trata de la región del campo que sufre mayor rotación.
5.5 Tensiones normales a una dirección dada
En un medio elástico lineal, isótropo y homogéneo los desplazamientos permiten escribir el tensor de tensiones a través de la fórmula:
- σ=λ*div(vec(u))*1+2με
- Siendo ε(vec(u)) la parte simétrica del gradiente de vect(u).
- Siendo λ y μ los coeficientes de Lamé que dependen de las propiedades elásticas de cada material, para este sólido consideraremos λ=μ=1
A partir del siguiente programa de MATLAB obtenemos las gráficas de las tensiones normales en las direcciones vec(gu) y vec(gv).
h=1/20;
u=1/3:h:1;
v=-1:h:1;
[Mu,Mv]=meshgrid(u,v);
x=Mu.*Mv;
y=0.5.*((Mu.^2)-(Mv.^2));
gu=(-14*Mv.*(Mu.^2)-6*Mv.^3)./(Mu.^2+Mv.^2); %Campo escalar
subplot(1,2,1)
surf(x,y,gu);
axis([-1,1,-1,1]);
view(2)
title('Tensiones normales a la dirección vec(gu) 2D');
subplot(1,2,2)
mesh(x,y,gu);
axis([-1,1,-1,1]);
view(3)
title('Tensiones normales a la dirección vec(gu) 3D');
h=1/20;
u=1/3:h:1;
v=-1:h:1;
[Mu,Mv]=meshgrid(u,v);
x=Mu.*Mv;
y=0.5.*((Mu.^2)-(Mv.^2));
gv=(6*(Mu.^2).*Mv-2*Mv.*3)./(Mu.^2+Mv.^2);
subplot (1,2,1);
surf(x,y,gv);
axis([-1,1,-1,1]);
view(2)
title('Tensiones normales a la dirección vec(gv) 2D');
subplot(1,2,2)
mesh(x,y,gv);
axis([-1,1,-1,1]);
view(3)
title('Tensiones normales a la dirección vec(gv) 3D');
Observamos que la gráfica del tensor normal en la dirección de vec(gu) es muy parecida a la de la divergencia, lo que nos indica que debe haber una relación entre estas dos gráficas; es decir que, la tensión superficial que soporta cada punto de la placa es muy similar al módulo de la divergencia en este mismo punto. Si interpretamos el módulo de la divergencia como la deformación que experimenta el cuerpo en dicho punto, tiene sentido relacionar estas dos magnitudes, defendiendo que la deformación que experimenta el punto es proporcional a la tensión que soporta, y que, por tanto, de esta manera se justifica la similitud entre sus gráficas. Sin embargo, no tenemos evidencias matemáticas para demostrar este hipótesis, por lo cual, esto no es más que una interpretación.
Sin embargo, la siguiente cuestión sería analizar el por qué aparece esta relación tensión-deformación de la dirección del vec(gu) y no en la del vec(gv). Una posible interpretación sería el que las tensiones son mucho menores, casi despreciables en esta dirección, y por ello prácticamente existe una correlación entre tensión-deformación en la dirección vec(gu), lo que explica la gran semejanza entre sus gráficas. Al igual que en el párrafo anterior, sólo contamos con el apoyo visual de las gráficas para defender este supuesto, y por tanto las conclusiones nunca serán definitiva; pero el color verde que aparece en la mayoría de gráfica de la tensión en la dirección vec(gv), lo que corresponde con una tensión en torno al cero, es decir tensión nula, sería un buen principio en la defensa de esta interpretación.
Por último, al comparar las gráficas del módulo del rotacional con la de las tensiones en ambas direcciones, vemos que no existe ninguna relación, al menos no a nivel visual, entre ellas. Lo cual a nivel teórico tiene sentido puesto que, la placa, como hemos indicado en los dos párrafos anteriores, se desplaza debido a las tensiones que soporta en las direcciones del vec(gu) y el vec(gv) (principalmente vec(gu)), estando ambas direcciones contenidas en el plano de la placa, es decir, la placa se desplaza en el plano que la contiene, en horizontal, y por lo tanto no gira, lo cual serviría como explicación a la falta se similitud entre las gráficas.
6 Tensión de Von Mises
La tensión de Von Mises (Richard von Mises, físico austrohúngaro,1913) es un indicativo o criterio de plasticidad, según el cuál esta se alcanza cuando las componentes de la tensión en un punto del sólido satisfacen la siguiente relación:
σVM =\sqrt{(((σ1 − σ2)^2 + (σ2 − σ3)^2 + (σ3 − σ1)^2)/2)}
Donde σ1,σ2,σ3,también llamados tensiones principales (de las infinitas tensiones que puede haber en un punto del sólido,las tensiones principales son aquellas que adoptan los valors máximos y mínimos) son los autovalores de σ, tensor de tensiones.
Procedemos al cálculo de la tensión de Von Mises en cada punto del sólido mediante el siguiente programa:
h=1/20;
u=1/3:h:1;
v=-1:h:1;
%Bucle para optimizar dimensiones
for i= 1:length(u)
for j= 1:length(v)
S=[-6.*v(j) -6.*u(i) 0.*v(j); -6.*u(i) -2.*v(j) 0.*v(j); 0.*v(j) 0.*v(j) -2.*v(j)];
%Calculo de autovalores del tensor de tensiones
autov=eig(S);
%Definición de la tensión de Von Mises
VM(i,j)= sqrt(((autov(1)-autov(2)).^2+(autov(2)-autov(3)).^2+(autov(3)-autov(1)).^2)/2);
end
end
%Mallado
[Mu,Mv]= meshgrid(u,v);
%Parametrización
x= Mu.*Mv;
y= 1/2*(Mu.^2-Mv.^2);
%Representación Tensión de Von Mises
surf(x,y,VM')
colorbar
axis([-1,1,-1,1]);
view(2)
title('Tensión de Von Mises')donde hemos representado,a su vez, la tensión de Von Mises sobre la superficie del sólido. Puede apreciarse que dicha tensión es máxima en los extremos superiores del sólido.
7 Cálculo de la masa de la placa
Para calcular la masa de la placa utilizamos el método del paralelogramo, quedando en MATLAB el siguiente programa:
N1=200; N2=100;
a=1/3; b=1; c=-1; d=1;
h1=(b-a)/N1; h2=(d-c)/N2;
u=a:h1:b; v=c:h2:d;
[Mu,Mv]=meshgrid(u,v);
f=0.5*(Mu.*Mv).*(Mu.^2-Mv.^2).*exp(-1./((Mu.*Mv).^2)).*(Mu.*2+Mv.^2);
w1=ones(N1+1,1);
w(1)=1/2; w(N1+1)=1/2;
w2=ones(N2+1,1);
w(1)=1/2; w(N2+1)=1/2;
result=h1*h2*w2'*f*w1
Como resultado a la integral, hemos obtenido el valor: -2.3180e-20. En un principio, el que la masa de la placa sea negativa parece erróneo, pero si comprobamos los valores que toma la función densidad en el dominio estudiado, observamos que en varios puntos este es negativo, lo cual, aunque no deja de ser físicamente contradictorio, explica que la masa de nuestra placa sea negativa.
