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

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. El foco está situado en el origen y por eso es la zona más caliente.

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([math]\nabla u[/math]) derivando respecto de las dos variables de las que depende. El gradiente indica la dirección en la cual un campo cambia más rápidamente, en este caso, en función de la temperatura.

%Como ya tenemos definida la malla y el campo, 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 es ortogonal a las curvas de nivel.

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

[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 incremento 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] .

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.

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 y=0,y=1 y en y=2 el rotacional es máximo.

8 Tensor de deformaciones

Definimos la parte simétrica del tensor gradiente de [math]\vec u[/math], que se denomina tensor de defrmaciones : [math]\epsilon(\vec u)=(\nabla \vec u + \nabla \vec u^t)/2[/math] 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 i*\epsilon*i, que da nulo pues es la proyección del vector \epsilon*i que solo depende de la coordenada j . Para hallar las tensiones normales en la dirección del eje j obtenemos el mismo resultado pues la proyección del vector \epsilon*j solo depende de la coordenada i.

Las tensiones tangenciales respecto al plano ortogonal a i vienen dadas por la ecuación [math]|σ⋅i→−(i→⋅σ⋅i→)i→|[/math]y respecto al plano ortogonal a j por la ecuación [math]|σ⋅j→−(j→⋅σ⋅j→)j→|[/math]

Así queda la gráfica de las tensiones tangenciales respecto al plano ortogonoal a i

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

Así queda la gráfica de las tensiones tangenciales respecto al plano ortogonal a j

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

figura 8: Tensiones tangenciales respecto al plano ortogonal a j

Comparando las dos gráficas observamos que en la figura 5 y en la figura 8 existe un máximo común en y=2; sin embargo las deformaciones mínimas son diferentes: Figura 5: y=0; y=1,5 // Figura 8: mínimo y=1.