Diferencia entre revisiones de «Estudio de campos escalares y vectoriales»

Revisión del 01:10 5 dic 2015

1 Introducción

En este documento se estudia el comportamiento de una placa plana al someterla a un campo escalar y un campo vectorial. Para una mejor interpretación nos ayudamos de gráficas que facilitan la comprensión de los resultados obtenidos.

2 Placa plana

La placa plana objeto de estudio se trata de una corona circular centrada en el origen limitada por las circunferencias de radios 1 y 2 unidades.

mallado de la placa

%Paso de muestreo
h=0.1;
%Valor de u según el muestreo en el intervalo [1,2]
u=1:h:2;
%valor de v según el muestreo en el intervalo [0, pi]
v=0:h:pi;
%Matrices de las componentes polares de la placa
[rho,theta]=meshgrid(u,v);
%parametrización
xx=rho.*cos(theta);
yy=rho.*sin(theta);
%Representacion del mallado
mesh(xx,yy,0*xx);
axis equal
axis([-3,3,-3,3]);
colormap('autumn')

2.1 Temperatura sobre la placa

Temperatura sobre la placa

curvas de nivel sobre la placa

%Paso de muestreo
h=0.1;
%Valor de u según el muestreo en el intervalo[1,2]
u=1:h:2;
%Valor de v según el muestreo en el intervalo[0,2*pi]
v=0:h:pi;
%Matriz de las componentes polares de la placa
[rho,theta]=meshgrid(u,v);
%Parametrización
xx=rho.*cos(theta);         
yy=rho.*sin(theta);
%Se dibuja el mallado         
mesh(xx,yy,0*xx)
%Se establecen los límites de los ejes x e y        
axis([-3,3,-3,3])
%Función temperatura
T=1./(yy+2);
%Representación
surf(xx,yy,T)
contour(xx,yy,T,10,'r')
axis([-3,3,-3,3])
colorbar
axis tight
%máximo global aproximado
[vmaxcol,vfila]=max(T);
[maxf,j]=max(vmaxcol);
i=vfila(j);
x0=xx(j)
y0=yy(vfila(j))
maxf

El máximo global aproximado está en [math](x,y)=(1,0)[/math],(en realidad de x=1 a x=2 el valor es cte) siendo el valor de la función en ese punto de [math]0.5[/math] unidades

3 Campo de las temperaturas

Ahora que se ha definido la placa, que será el dominio de trabajo, se puede estudiar el efecto de la aplicación del campo temperatura en la placa. El campo temperatura es un campo escalar que depende de las coordenadas x e y del punto en que se aplique.: [math]T(x,y)=1/(y+2)[/math]Para entender como afecta este campo es útil calcular calcular su variación en el dominio, es decir, calcular el gradiente del campo (Campo vectorial).

3.1 Variación de la temperatura[math]∇T(x,y)[/math]

[math]∇T(x,y)=0\;\vec i -1/(y+2)^2\;\vec j [/math]El gradiente, en los puntos en los que se aplique, indica la dirección en que crece más la temperatura. Las curvas de nivel (también llamadas curvas equipotenciales o isocuantas) del campo escalar son aquellas que cumplen [math]T(x,y)=cte[/math]. Por definición [math]∇T(x,y)[/math] es ortogonal a las curvas de nivel en cada punto. En la gráfica se representa el campo vectorial [math]∇T(x,y)[/math] y las curvas equipotenciales. Se observa cómo varían la dirección (siempre ortogonal a la curva de nivel) y el módulo (que representa cuánto varía la función) del vector gradiente excepto en el máximo global.

Gradiente y Curvas de nivel sobre la placa

%Paso de muestreo
h=0.1;
%Valor de u según el muestreo en el intervalo[1,2]
u=1:h:2;
%Valor de v según el muestreo en el intervalo[0,2*pi]
v=0:h:pi;
%Matriz de las componentes polares de la placa
[rho,theta]=meshgrid(u,v);
%Parametrización
xx=rho.*cos(theta);         
yy=rho.*sin(theta);
%Se dibuja el mallado         
mesh(xx,yy,0*xx)
%Se establecen los límites de los ejes x e y     
axis([-3,3,-3,3]) 
hold on
%Campo temperatura
T=1./(yy+2);
%Representación del gradiente y las curvas de nivel
contour(xx,yy,T,10,'r')
dx= 0*xx;
dy= -1./(yy+2).^2;
quiver(xx,yy,dx,dy)
hold off
axis tight

4 Campo de los desplazamientos [math]\vec u[/math]

A continuación se estudia el efecto que produce el campo desplazamiento [math] \vec u(\rho,\theta)[/math] en la placa. El campo desplazamiento es el campo vectorial expresado matemáticamente por la siguiente función vectorial.:[math] \vec u(\rho,\theta)={(1-\rho)^2}\vec g_{\rho}[/math] Para entender cómo afecta este campo vectorial a la placa se analiza el cambio que provoca sobre la placa asi como las propiedades que posee el campo por su definición. Antes de seguir hay que decir que, para poder trabajar comodamente con el campo, este se define en el espacio, es decir, se trabajará en coordenadas cilíndricas asumiendo que el campo no tiene componente [math]\vec g_{z}[/math].:[math] \vec u(\rho,\theta,z)={(1-\rho)^2}\vec g_{\rho}[/math]También se define un concepto que ayudará en el análisis de las propiedades del campo, la derivada parcial covariante. Sea [math] \left \{\vec g_{\rho},\;\vec g_{\theta},\;\vec g_{z} \right \}[/math] la base natural asociada al sistema de coordenadas curvilíneas[math]\left ( x^1,\;x^2,\;x^3\right )=\left(\rho,\;\theta,\;z \right )[/math] y [math]\vec u=u^i\;\vec g_{i}[/math] un campo vectorial, las derivadas parciales covariantes de [math]\vec u[/math], [math]u_j^i[/math] vienen dadas por:[math]u_j^i=\displaystyle\frac{\partial\;u^i}{\partial\;x^j}+\Gamma_{kj}^i\;u^k[/math] que representadas en forma de matriz son:[math]u_j^i=\begin{bmatrix} 2(\rho-\;1) & 0 & 0 \\ 0 & \displaystyle\frac{1}{\rho}(1-\;\rho)^2 & 0 \\ 0 & 0 & 0 \end{bmatrix}[/math]

4.1 Desplazamiento de la placa

Campo de desplazamientos sobre la placa

%Paso de muestreo
h=0.1;
%Valor de u según el muestreo en el intervalo [1,2]
u=1:h:2;
%Valor de v según el muestreo en el intervalo[0,2*pi]
v=0:h:pi;
%Matriz de las componentes polares de la placa
[rho,theta]=meshgrid(u,v);
%Parametrización
xx=rho.*cos(theta);
yy=rho.*sin(theta);
%Componente x(rho,theta) del campo vectorial
fx=-(3./(20.*rho)).*sin(2.*pi.*rho/.3-pi/3).*rho.*sin(theta);
%Componente y(rho,theta) del campo vectorial
fy=(3./(20.*rho)).*sin(2.*pi.*rho/.3-pi/3).*rho.*cos(theta);
hold on;
%Representación del campo vectorial
quiver(xx,yy,fx,fy);
axis equal
axis([-3,3,-3,3]);
%Representación del sólido
mesh(xx,yy,0*xx);
axis tight
colormap('autumn')

Si se analiza la función vectorial es fácil comprobar que si [math]\rho\;\varepsilon\;[1,2][/math] esta función alcanza su máximo y mínimo valor en [math]\rho=2[/math] y [math]\rho=1[/math] respectivamente. En la gráfica se representa el cambio que sufre la placa al aplicar el campo desplazamiento y, como queda representado en la figura [math]\vec u(1,\theta)= 0\;\vec g_{\rho}[/math] y [math]\vec u(2,\theta)= 1\;\vec g_{\rho}[/math].

Mallado antes y después de aplicar el campo [math]\vec u[/math]

%Paso de muestreo
h=0.1;
%Valor de u según el muestreo en el intervalo [1,2]
u=1:h:2;
%Valor de v según el muestreo en el intervalo[0,2*pi]
v=0:h:pi;
%Matriz de las componentes polares de la placa
[rho,theta]=meshgrid(u,v);
%Parametrización
xx=rho.*cos(theta);
yy=rho.*sin(theta);
%Componente x(rho,theta) del campo vectorial
fx=-(3./(20.*rho)).*sin(2.*pi.*rho/.3-pi/3).*rho.*sin(theta);
%Componente y(rho,theta) del campo vectorial
fy=(3./(20.*rho)).*sin(2.*pi.*rho/.3-pi/3).*rho.*cos(theta);
subplot(1,2,1);
mesh(xx,yy,0*xx);
axis tight
%Representación del mallado desplazado
subplot(1,2,2);
mesh(xx+fx,yy+fy,0*xx);
axis tight
colormap('autumn')

4.2 Variación del volumen

La primera propiedad del campo en ser analizada es la divergencia. El sentido físico de la divergencia en un punto es la tasa de flujo neto hacia el exterior, en el punto, por unidad de volumen. Como el campo analizado provoca el desplazamiento relativo de los puntos que forman la placa, la divergencia expresará el cambio de área local debido al desplazamiento. La divergencia es nula por lo tanto la placa es incompresible.

Representación de la divergencia

%Paso de muestreo
h=0.1;
%Valor de u según el muestreo en el intervalo [1,2]
u=1:h:2;
%valor de v según el muestreo en el intervalo [0,2*pi]
v=0:h:pi;
%Matriz de las componentes polares de la placa
[rho, theta]=meshgrid(u,v);
%Parametrización
xx=rho.*cos(theta);
yy=rho.*sin(theta);
%Campo escalar de la divergencia
f=0.*xx;
%Representación de la divergencia sobre la placa
surf(xx,yy,f);
colorbar
axis equal

4.3 Tendencia a la rotación

Por último se analiza otra propiedad del campo, el rotacional. El rotacional de un campo vectorial en un punto es la tendencia del campo a inducir rotación alrededor de dicho punto. En este caso los puntos son los que pertenezcan al dominio de estudio (placa plana).:[math]\nabla \times \vec u = \vec 0[/math]

Código:

Rotacional

%Paso de muestreo
h=0.1;
%Valor de u según el muestreo en el intervalo [1,2]
u=1:h:2;
%valor de v según el muestreo en el intervalo [0,2*pi]
v=0:h:pi;
%Matriz de las componentes polares de la placa
[rho, theta]=meshgrid(u,v);
%Parametrización
xx=rho.*cos(theta);
yy=rho.*sin(theta);
% Rotacional
rot = rho.*(pi/5).*cos((4.*rho./3)-2*pi/3);
surf(xx, yy, rot)
title('Rotacional ')

.

4.4 Tensiones originadas

En un medio elástico lineal e isótropo los desplazamientos permiten escribir el tensor de tensiones [math]\sigma_j^i[/math] a través de la fórmula: [math]\sigma_j^i=\lambda\nabla\cdot\vec u\;\delta_j^i + 2\mu\;\epsilon_j^i\;[/math] donde [math]\lambda[/math] y [math]\mu[/math] son los conocidos como coeficientes de Lamé que dependen de las propiedades elásticas de cada material. En este caso dichos coeficientes tendrán el mismo valor [math]\lambda=\mu=1[/math] y [math]\epsilon_j^i[/math] la matriz 1-contravariante 1-covariante del tensor de deformaciones, definido como la parte simétrica del tensor gradiente de [math]\vec u[/math] siendo [math]\nabla \vec u=u_j^i\;\vec g_i\;\otimes\;\vec g^j[/math].:[math]\sigma_j^i=\begin{bmatrix} 6(\rho-1)+\frac{1}{\rho}(1-\rho)^2 & 0 & 0 \\ 0 & 2(\rho-1)+\frac{3}{\rho}(1-\rho)^2 & 0 \\ 0 & 0 & 2(\rho-1)+\frac{1}{\rho}(1-\rho)^2 \end{bmatrix}[/math]El esfuerzo interno es equivalente a las tensiones internas en la placa. Estas se dividen en tensiones normales, que provocan esfuerzos normales, y tensiones tangenciales, que provocan esfuerzos cortantes.

4.4.1 Tensiones normales en la direcciones que marcan los ejes [math]\vec g_\rho[/math] y [math]\vec g_\theta/\rho[/math]

La resultante de tensiones normales se puede proyectar en tres ejes generados por vectores que formen base del espacio [math]R^3[/math]. En este caso, como la placa es plana es interesante conocer la proyección de la resultante de tensiones normales en los ejes del espacio generados por [math]\vec g_\rho[/math] y [math]\vec g_\theta/\rho[/math]. Para hallar estas proyecciones hay que aplicar la defición de base recíproca de una base ortogonal y la representación diádica del tensor [math]\sigma[/math]:[math]\sigma=\sigma_j^i\;\vec g_i\;\otimes\;\vec g^j[/math]La resultante de tensiones normales en la dirección del eje generado por [math]\vec g_\rho[/math] se calcula mediante la operación [math]\vec g_{\rho} \cdot \sigma \cdot \vec g_\rho[/math]:[math]\vec g_\rho \cdot \sigma \cdot \vec g_\rho=6(\rho-1)+\frac{1}{\rho}(1-\rho)^2[/math]La resultante de tensiones normales en la dirección del eje generado por [math]\vec g_\theta/\rho[/math]:[math]\vec g_\theta/\rho\cdot \sigma \cdot \vec g_\theta/\rho=2(\rho-1)+\frac{3}{\rho}(1-\rho)^2[/math]Matemáticamente se obtiene que la proyección sobre el eje generado por [math]\vec g_\rho[/math] es mayor que la obtenida sobre el eje que genera [math]\vec g_\theta/\rho[/math], esto se puede visualizar en las siguientes gráficas obtenidas en MatLAB.

Representación de las tensiones normales en la dirección [math]\vec g_\rho[/math]

Representación de las tensiones normales en la dirección [math]\vec g_\theta/\rho[/math]

%Paso de muestreo
h=0.1;
%Valor de u según el muestreo en el intervalo [1,2]
u=1:h:2;
%Valor de v según el muestreo en el intervalo [0, 2*pi]
v=0:h:pi;
%Matrices de las componentes polares de la placa
[rho, theta]=meshgrid(u,v);
%Parametrización
xx=rho.*cos(theta);
yy=rho.*sin(theta);
figure(1)
%Representación de la función sigma en la direccion de gsubrho
sigmaunouno=0.*xx; 
surf(xx,yy,sigmaunouno)
axis tight
colorbar
%Representación de la función sigma en la direccion de gsubrthetaentrerho
figure(2)
sigmadosdos=0.*xx;
surf(xx,yy,sigmadosdos)
axis tight
colorbar

4.4.2 Tensiones tangenciales respecto del plano ortogonal al vector [math]\vec g_\rho[/math] y al vector [math]\vec g_\theta/\rho[/math]

Se calculan las tensiones tangenciales analíticamente y se dibujan como se ha hecho anteriormente

Representación de las tensiones normales en la dirección [math]\vec g_\rho[/math]

Representación de las tensiones normales en la dirección [math]\vec g_\theta/\rho[/math]

%Paso de muestreo
h=0.1;
%Valor de u según el muestreo en el intervalo [1,2]
u=1:h:2;
%Valor de v según el muestreo en el intervalo [0, 2*pi]
v=0:h:pi;
%Matrices de las componentes polares de la placa
[rho, theta]=meshgrid(u,v);
%Parametrización
xx=rho.*cos(theta);
yy=rho.*sin(theta);
figure(1)
%Representación de la función sigma en la direccion de gsubrho
sigmaunouno=0.*xx; 
surf(xx,yy,sigmaunouno)
axis tight
colorbar
%Representación de la función sigma en la direccion de gsubrthetaentrerho
figure(2)
sigmadosdos=0.*xx;
surf(xx,yy,sigmadosdos)
axis tight
colorbar

4.4.3 Tensión de Von Mises

Un elemento importante para comprender el comportamiento de la placa es la Tensión de Von Mises. Se obtiene con 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]Donde [math] \sigma_1 [/math], [math] \sigma_2 [/math] y [math] \sigma_3 [/math] son los autovalores de [math] \sigma [/math], también conocidos como tensiones principales. Esta magnitud escalar se usa como indicador en los ensayos de materiales para saber cuando un material inicia un comportamiento plástico (y no elástico puro). Con la siguente rutina en MatLAB se consigue representar la Tensión de Von Mises y su evolución, junto con la de los autovalores, a medida que aumenta [math]\rho[/math]

%Paso de muestreo
h=0.1;
%Valor de u según el muestreo en el intervalo [1,2]
u=1:h:2;
%valor de v según el muestreo en el intervalo [0, 2*pi]
v=0:h:2*pi+h;
%Matrices de las componentes polares de la placa
[ro, theta]=meshgrid(u,v);
%parametrización
xx=ro.*cos(theta);
yy=ro.*sin(theta);
%Creación de la matriz de autovalores V
V=[];
x=1;
for rho=1:0.1:2
    A=[6*(rho-1)+((1./rho).*((1-rho).^2)),0,0;0,2*(rho-1)+((3/rho).*((1-rho).^2)),0;0,0,2*(rho-1)+((1./rho).*((1-rho).^2))];
    [X,D]=eig(A);
    t=diag(D);
    V(:,x)=t;
    x=x+1;
end
% Vectores a,b y c que contienen los autovalores para cada valor rho
a=V(1,:); b=V(2,:); c=V(3,:);
%Tensión de Von Mises
VM=sqrt(((a-b).^2+(b-c).^2+(c-a).^2)./2);
%Generación de la matriz TVM
TVM=[]; 
for i=1:64
 TVM(i,:)=VM;
end
figure(1)
%Representación de la tensión de Von Mises sobre la placa
mesh(xx,yy,TVM);
colorbar
axis equal
figure(2)
%Representación de los autovalores junto con la tensión de Von Mises
hold on
plot(u,a,'r')
plot(u,b,'g')
plot(u,c)
plot(u,VM,'k')
legend('Primer autovalor en función de rho','Segundo autovalor en función de rho','Tercer autovalor en función de rho','Tensión de Von Mises según rho','location','best')
hold off

5 Masa de la placa

El cálculo de la masa de la placa, con la función de densidad [math]d(x,y,z)=xye^{-1/x^2}[/math], se realiza usando un método numérico para el cálculo aproximado de integrales definidas, el método del trapecio, con la siguente rutina en MatLAB. Además esta rutina permite representar la masa y su distribución en la placa según la función de densidad.

Distribución de masa en la placa

%número de puntos
N1=100;
N2=100;
%extremos de intervalos
a=1; b=2;
c=0;d=pi;
h1=(b-a)/N1; %variable en rho
h2=(d-c)/N1; %variable en phi
%definimos la partición
u=a:h1:b; v=c:h2:d;
%generamos el rectángulo
[uu,vv]=meshgrid (u,v);
x=uu.*cos(vv);
y=uu.*sin(vv);
%insertamos la función
f=abs(log(uu.^2+2));
w1=ones(N1+1,1);
w1(1)=1/2;  w1(N1+1)=1/2;
w2=ones (N2+1,1);
w2(1)=1/2; w2(N2+1)=1/2;
resultado=h1*h2*w2'*f*w1;
Masa=resultado;Masa
mesh(x,y,f)
axis image

La masa obtenida es [math]m=4.5431 unidades[/math]. Se observa cómo aumenta la densidad en los puntos que cumplen [math]x=y[/math] a medida que están más alejados del origen, debido a la componente exponencial de la función [math]d(x,y,z)[/math].