Visualizacion de campos escalares y vectoriales (grupo 3A)

Visualización de campos escalares y vectoriales en elasticidad.

Marta de Castro Pérez/ Alejandra García-Page Acevedo/ Silvia Pinedo Gil/ Ana María Ragolta Villarroya/

1 Enunciado

Consideramos una placa rectangular plana (en dimensión 2) que ocupa la región [math] [-1/2,1/2] \times [0,2] [/math]. En ella vamos a suponer que tenemos definidas dos cantidades físicas: la temperatura [math]T(x,y,t)[/math], que depende de las dos variables espaciales [math](x,y)[/math] y el tiempo [math]t[/math], y los desplazamientos [math]\vec u(x,y,t)[/math]. De esta forma, si definimos [math]r_0(x,y)[/math] el vector de posición de los puntos de la placa en reposo, la posición de cada punto [math](x,y)[/math] de la placa en un instante de tiempo [math]t[/math] viene dada por: [math] \vec r (x,y,t)= \vec r_{0}(x,y)+\vec u(x,y,t). [/math] Vamos a suponer que sobre la placa se ha aplicado una fuerza que ha provocado una vibración de manera que los desplazamientos vienen dados por la onda: [math] \vec u(x,y,t) = \vec a \sin(\vec b \cdot \vec r_0-ct), [/math] donde [math]\vec a[/math] se conoce como amplitud, [math]\vec b[/math] es la fase que indica la dirección de propagación y [math]c/|\vec b|[/math] es la velocidad de propagación.

Si [math]\vec a [/math] es paralelo a [math]\vec b[/math] diremos que la onda es longitudinal mientras que si es perpendicular hablaremos de onda transversal.

En este trabajo vamos a centrarnos en las ondas transversales. Supondremos lo siguiente: [math] \vec a=\frac{\vec i}{10}, \qquad \vec b= \pi \vec j, \qquad t=0. [/math] En este caso, [math]\vec u(x,y)=\frac{\sin(\pi y)}{10}\vec i[/math].

2 Mallado

Definimos el mallado representado por los puntos interiores del sólido.

x=-0.5:0.1:0.5;       % Creamos los intervalos x= [-1/2,1/2] y=[0,2]
y=0:0.1:2;             
[xx,yy]=meshgrid(x,y); %Creamos la matriz a partir del intervalo 
mesh(xx,yy,0*xx)       % Dibujamos malla
axis([-2,2,-1,3])      % seleccionamos región

Figura 1 Mallado de la placa

3 Temperatura del sólido

Como ya hemos definido previamente las variables x e y creamos directamente el campo temperatura ([math]T(\rho,\theta)=-\log(\rho+0.1)[/math] ) a partir de esas variables. Sin embargo esta función viene dada en coordenadas polares por lo que realizamos el cambio de variable a coordenadas cartesianas: [math] \rho =\sqrt{x^2 + y^2} [/math] y obtenemos la función:

[math] T(x,y)=-\log(\sqrt{x^2+y^2}+0.1)[/math]

T=-log(sqrt(xx.^2+yy.^2)+0.1); 
surf(xx,yy,f)          % Dibujamos la superficie 
axis([-2,2,-1,3])      % Seleccionamos la región

En el gráfico las altas temperaturas aparecen en rojo y las más bajas en azul. En el origen se encuentra el foco, y por eso es la zona con mayor temperatura.

Figura 2 Distribución de la temperatura en la malla

Figura 3 Distribución de la temperatura en la malla vista en planta

4 Gradiente y curvas del campo

Calculamos el gradiente del campo vectorial T ([math]\nabla T[/math]) derivando respecto de las dos variables de las que depende (x,y). El gradiente indica la dirección en la que el campo T cambia más rápidamente en función de la temperatura. Es decir, indica la mínima distancia entre puntos de diferente temperatura, en dirección al foco.

%Como ya tenemos definida la malla y el campo T, calculamos directamente el gradiente
Tx=-xx./((sqrt(xx.^2+yy.^2)).*(sqrt(xx.^2+yy.^2)+0.1));
Ty=-yy./((sqrt(xx.^2+yy.^2)).*(sqrt(xx.^2+yy.^2)+0.1));
%Dibujamos el Vector gradiente (figura 3)
quiver(xx,yy,Tx,Ty)
hold on 
contour(xx,yy,T)       % lineas de nivel

plot(x,x-x,'m','linewidth',1);
plot(x,2+x-x,'m','linewidth',1);
plot(-0.5+y-y,y,'m','linewidth',1);
plot(0.5+y-y,y,'m','linewidth',1);
axis([-1,1,-0.5,2.5]) %Numeramos los ejes

Figura 4 Representación del gradiente y las curvas de nivel del campo

Se puede observar gráficamente que el gradiente de T([math]\nabla T[/math]) es ortogonal a las curvas de nivel del campo.

5 Campo de vectores adicional

Consideremos ahora el campo de vectores [math]\vec u(x,y)=\frac{\sin(\pi y)}{10}\vec i[/math] y lo dibujamos en los puntos del mallado del sólido para determinar el desplazamiento de los puntos del mallado.

%U(xx,yy)=(0.1*sin(pi.*yy))i/10 + 0.*xx j
%definimos las componentes del campo U
Ux=(0.1*sin(pi.*yy))/10;
Uy=0.*xx;
quiver(xx,yy,Ux,Uy)
axis([-1,1,-0.5,2.5])
view([0,0,1])

Figura 4 Campo de vectores en los puntos del mallado del sólido

6 Desplazamiento sobre el mallado

A continuación consideramos el vector [math]\vec u(x,y)=\frac{\sin(\pi y)}{10}\vec i[/math] como un vector de desplazamiento sobre la malla y comparamos con la figura 1(antes del desplazamiento).

[xx,yy]=meshgrid(x,y); %dibujamos la malla inicial 
subplot(1,2,1)
mesh(xx,yy,0*xx)
axis([-1,1,-0.5,2.5]) %malla con el vector de desplazamiento aplicado
subplot(1,2,2)
Ux=0.1*sin(pi*yy);
Uy=0*xx;
mesh(xx+Ux,yy+Uy,0*xx)
axis([-1,1,-0.5,2.5]) %Numeración de ejes

Figura 5 Comparación de la malla sin y con el campo de desplazamiento

7 Cálculo de la divergencia y el rotacional

La divergencia resulta nula puesto que el vector u sólo tiene componente en x, y esa componente depende de la variable y. La divergencia representa en este caso el cambio de volumen en cada punto del sólido bajo la acción del campo vectorial [math]\vec u(x,y)=\frac{\sin(\pi y)}{10}\vec i[/math] (debido al desplazamiento).

En nuestra malla la divergencia es nula porque al tratarse de un plano todos los puntos se mueven de la misma forma,lo que hace que el cambio (desplazamiento) sea igual en todos los puntos.

El caso del rotacional es distinto. El rotacional de un vector representa la tendencia de un campo vectorial a rotar alrededor de un punto o eje. Los puntos con mayor rotacional serán los que tengan mayores derivadas parciales.

R=abs((0.1*pi).*cos(pi.*yy)); %función rotacional de U
surf(xx,yy,R)
axis([-1,1,-0.5,2.5]) %numeramos los ejes
view([0,0,1])
colorbar %añadimos la leyenda de colores

Figura 6 Representación del rotacional

Se puede observar que en las rectas y=0,y=1 y en y=2 el rotacional es máximo. Los puntos contenidos en estas rectas son más susceptibles a sufrir una rotación en torno a un eje perpendicular a estas rectas.

8 Tensor de deformaciones

Definimos el tensor [math]\epsilon(\vec u)=(\nabla \vec u + \nabla \vec u^t)/2[/math], que es la parte simétrica del tensor gradiente de [math]\vec u[/math], que se denomina tensor de deformaciones. Todo esto en un medio elástico lineal, isótropo y homogéneo donde los desplazamientos permiten escribir el tensor de tensiones [math]\sigma_{ij}[/math] a través de la fórmula: [math] \sigma_{ij}=\lambda \nabla \cdot \vec u \delta_{ij} + 2\mu \epsilon_{ij}, [/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.

Tomando [math]\lambda=\mu=1[/math], dibujamos las tensiones normales en la dirección que marca el eje [math]\vec i[/math], es decir [math]\vec i \cdot \sigma \cdot \vec i[/math] y las tensiones normales en la dirección que marca el eje [math]\vec j[/math], es decir [math]\vec j \cdot \sigma \cdot \vec j[/math].

Para hallar las tensiones normales en la dirección del eje i ,aplicamos [math]\vec i \cdot \sigma \cdot \vec i[/math] , que da nulo pues es la proyección del vector [math]\sigma \cdot \vec i[/math], que solo depende de la coordenada [math]\vec j[/math]. Las tensiones normales en la dirección del eje [math]\vec j[/math] son también nulas, pues la proyección del vector [math]\sigma \cdot \vec j[/math], solo depende de la coordenada [math]\vec i[/math].

Hallamos las tensiones tangenciales:

Las tensiones tangenciales respecto al plano ortogonal a [math]\vec i[/math] vienen definidas por la expresión [math]|\sigma \cdot \vec i-(\vec i \cdot \sigma \cdot \vec i) \vec i|[/math] y las tensiones tangenciales respecto al plano ortogonal a [math]\vec j[/math] , [math]|\sigma \cdot \vec j-(\vec j \cdot \sigma \cdot \vec j) \vec j|[/math].

Así queda la gráfica de las tensiones tangenciales respecto al plano ortogonoal a [math]\vec i[/math] :

Ux=0;
Uy=(pi*cos(pi*yy))/10;
quiver(xx,yy,xx*0,Uy)
axis([-1,1,-0.5,2.5])
view([0,0,1])

figura 7: Tensiones tangenciales respecto al plano ortogonal a i

En comparación con la figura 5, observamos que las zonas de mayor deformación coinciden con las zonas de menor tensión y viceversa.

Así queda la gráfica de las tensiones tangenciales respecto al plano ortogonal a [math]\vec j[/math] :

Uy=0;
Ux=(pi*cos(pi*yy))/10;
quiver(xx,yy,Ux,0*yy)
axis([-1,1,-0.5,2.5])
view([0,0,1])

figura 8: Tensiones tangenciales respecto al plano ortogonal a [math]\vec j[/math]

En comparación con la figura 5, observamos que las zonas máximas y mínimas son las mismas.