Visualización de campos escalares y vectoriales en elasticidad grupo 14A
Contenido
1 Enunciado
Consideramos una placa plana (en dimensión 2) que ocupa el anillo comprendido entre las circunferencias centradas en el origen de radios 2 y 1. En ella vamos a suponer que tenemos definidas dos cantidades físicas: la temperatura [math]T(\rho,\theta,t)[/math], que depende de las dos coordenadas polares [math](\rho,\theta)[/math] y el tiempo [math]t[/math], y los desplazamientos [math]\vec u(\rho,\theta,t)[/math]. De esta forma, si definimos [math]r_0(\rho,\theta)[/math] el vector de posición de los puntos de la placa en reposo, la posición de cada punto [math](\rho,\theta)[/math] de la placa en un instante de tiempo [math]t[/math] viene dada por: [math] \vec r (\rho,\theta,t)= \vec r_{0}(\rho,\theta)+\vec u(\rho,\theta,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 en un tiempo [math]t_0[/math] dado vienen dados por: [math] \vec u(\rho,\theta)=\frac{\sin(\pi \theta/2)}{30\rho}\vec g_{\rho}. [/math]
Se pide:
- Dibujar un mallado que represente los puntos interiores del sólido (ver figura). Tomar como paso de muestreo [math]h=1/10[/math] para las variables [math]x[/math] e [math]y[/math].
- La temperatura del sólido proviene de un foco de calor muy concentrado en el origen. Supongamos que la conocemos y viene dada por [math]T(\rho,\theta)=-\log(\rho+0.1)[/math]. Dibujarla.
- Calcular [math]\nabla T[/math] y dibujarlo como campo vectorial. Calcular las curvas de nivel de T y observar gráficamente que [math]\nabla T[/math] es ortogonal a dichas curvas.
- Consideramos ahora el campo de vectores [math]\vec u=\frac{\sin(\pi \theta/2)}{30\rho}\vec g_{\rho}[/math]. Dibujar el campo de vectores en los puntos del mallado del sólido.
- Si [math]\vec u[/math] determina el desplazamiento que sufre cada punto del sólido, dibujar el sólido antes y después del desplazamiento. Dibujar ambos en la misma figura usando el comando subplot.
- Calcular la [math]\nabla \cdot \vec u[/math] en todos los puntos del sólido y dibujarla. ¿Qué puntos tienen mayor divergencia? La divergencia es una medida del cambio de volumen local debido al desplazamiento.
- Calcular [math]|\nabla \times \vec u|[/math] en todos los puntos del sólido y dibujarlo. ¿Qué puntos sufren un mayor rotacional?
- Definamos [math]\epsilon(\vec u)=(\nabla \vec u + \nabla \vec u^t)/2[/math], la parte simétrica del tensor gradiente de [math]\vec u[/math], que se denomina tensor de deformaciones. En un medio elástico lineal, isótropo y homogéneo 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], dibujar las tensiones normales en la dirección [math]\vec g_{\rho}[/math], es decir [math]\vec g_{\rho} \cdot \sigma \cdot \vec g_\rho[/math] y las tensiones normales en la dirección [math]\vec g_\theta/\rho[/math], es decir [math]\vec g_\theta/\rho \cdot \sigma \cdot \vec g_\theta/\rho[/math].
- Dibujar las tensiones tangenciales respecto al plano ortogonal a [math]\vec g_\rho[/math], es decir [math]|\sigma \cdot \vec g_\rho-(\vec g_\rho \cdot \sigma \cdot \vec g_\rho) \vec g_\rho|[/math]. ¿Donde son mayores? Compararlas con los puntos de mayor deformación de la malla.
- Dibujar las tensiones tangenciales respecto al plano ortogonal a [math]\vec g_\theta/\rho[/math], es decir [math]|\sigma \cdot \vec g_\theta/\rho-(\vec g_\theta/\rho \cdot \sigma \cdot \vec g_\theta/\rho) \vec g_\theta/\rho|[/math]. ¿Donde son mayores? Compararla con los puntos de mayor deformación de la malla.
2 Mallado del sólido
El mallado del sólido reprsenta el dominio en el cual se encuentran los puntos que componen el sólido. Por lo tanto, para realizar el mallado de deben crear dos matrices cuyos componentes serán distintos valores que pueden tomar cada una de las variables (ro, theta). Después por interpolación hallaremos el cojunto total de los puntos que componen el sólido, y procedemos a dibujar el mallado.
h=0.1; %Establecemos paso de muestra
u=1:h:2; %Valor de u según el muestreo del intervalo
v=0:h:2*pi+h; %Valor de v según el muestreo del intervalo
[uu,vv]=meshgrid(u,v); %Matrices de la placa en coordenadas polares
figure(1) %Abrimos la ventana donde dibujaremos el mallado
xx=uu.*cos(vv); %Parametrización
yy=uu.*sin(vv); %Parametrización
mesh(xx,yy,0*xx) %Dibujamos el mallado
axis([-3,3,-3,3]) %Establecemos los límites de los ejes x e y en la representación de la placa
view(2) %visualización
3 Campo de temperaturas
La temperatura del sólido varía a lo largo de este según la función -log(ro+0.1).
h=0.1; %Establecemos paso de muestra
u=1:h:2; %Valor de u (ro) según el muestreo del intervalo
v=0:h:2*pi+h; %Valor de v (theta) según el muestreo del interval
[uu,vv]=meshgrid(u,v); %Matrices de la placa en coordenadas polares
figure(1) %Abrimos la ventana donde dibujaremos el mallado
xx=uu.*cos(vv); %Parametrización
yy=uu.*sin(vv); %Parametrización
f=-log(uu+0.1); %Campo escalar
surf(xx,yy,f) %Dibujamos la superficie
axis([-3,3,-3,3]) %Establecemos los límites de los ejes x e y en la representación de la placa
colorbar %Hacemos aparecer en la figura la leyenda de los colores
view(2) %Visualización
Como se observa en la figura la temperatura del sólido es mayor en su interior y menor en su exterior. Todos los puntos pertenecientes a una misma circunferencia tendrán la misma temperatura. Es decir, las circunferencias que componen el sólido de mayor radio tendrán menor temperatura que las de menor radio. Las circunferencias de los extremos contendrán a los puntos con que tienen la mínima y máxima temperatura (radio 1, máxima temperatura; radio 2 mínima)
4 Gradiente del campo de temperaturas
El gradiente es el campo vectorial cuyas funciones coordenadas son las derivadas parciales del campo escalar dado, que representa la temperatura en cada punto de nuestro sólido. En este caso, el gradiente en cada uno de los puntos apunta en la dirección en la cual la temperatura aumenta más rápido (es decir, la dirección en que la derivada direccional es máxima). El módulo del gradiente indica la rapidez con que aumenta la temperatura en esa dirección.
En este caso el gradiente del campo es:
[math]\nabla T = (\frac{-x}{x^2+y^2+0.1\sqrt{x^2+y^2}},\frac{-y}{x^2+y^2+0.1\sqrt{x^2+y^2}})[/math]
Como sabemos por las propiedades del gradiente, el vector gradiente será perpendicular a las líneas de contorno (líneas "equiescalares", definidas por ϕ =cte) del sólido. Podemos observarlo fácilmente en la figura en 3D del sólido que se muestra más abajo.Vemos claramente en los dibujos que estas líneas equiescalares son circunferencias concéntricas que muestran la disminución de la temperatura a medida que los puntos del sólido se alejan del foco.
h=0.1; %Establecemos paso de muestra
u=1:h:2; %Valor de u según el muestreo del intervalo
v=0:h:2*pi+h; %Valor de v según el muestreo del intervalo
[uu,vv]=meshgrid(u,v); %Matrices de la placa en coordenadas polares
figure(1) %Abrimos la ventana donde dibujaremos el mallado
xx=uu.*cos(vv); %Parametrización
yy=uu.*sin(vv); %Parametrización
mesh(xx,yy,0*xx) %Dibujamos el mallado
axis([-3,3,-3,3]) %Establecemos los límites de los ejes x e y en la representación de la placa
view(2) %visualización
f=-log(uu+0.1); %Campo escalar
surf(xx,yy,f) %Dibujamos la superficie
axis([-3,3,-3,3]) %Establecemos los límites de los ejes x e y en la representación de la placa
colorbar %Hacemos aparecer en la figura la leyenda de los colores
view(2) %Visualización
hold on
grad=-log(uu+0.1); % Introducimos la nueva función “grad” que representa el gradiente del campo de temperaturas
contour(xx,yy,grad,10,'g') % Dibujamos 10 líneas de nivel en color verde
dx=-xx./((xx.^2+yy.^2)+0.1.*sqrt(xx.^2+yy.^2)); %Introducimos la función de la derivada parcial respecto de x que hemos calculado previamente
dy=-yy./((xx.^2+yy.^2)+0.1.*sqrt(xx.^2+yy.^2)); %Introducimos la función de la derivada parcial respecto de y que hemos calculado previamente
quiver(xx,yy,dx,dy) %Dibujamos el campo vectorial del gradiente de temperaturas
xlabel('x'); %Nombramos el eje x
ylabel('y'); %Nombramos el eje y
zlabel('Temperatura'); %Nombramos el eje z
5 Campo de vectores del desplazamiento del sólido
Supongamos que cada punto del sólido es desplazado en función de u. El campo vectorial que afecta al sólido viene representado sobre el mallado en la siguiente figura.
h=0.1; %Establecemos paso de muestra
u=1:h:2; %Valor de u según el muestreo del intervalo
v=0:h:2*pi+h; %Valor de v según el muestreo del intervalo
[uu,vv]=meshgrid(u,v); %Matrices de la placa en coordenadas polares
figure(1) %Abrimos la ventana donde dibujaremos el mallado
xx=uu.*cos(vv); %Parametrización
yy=uu.*sin(vv); %Parametrización
mesh(xx,yy,0*xx) %Dibujamos el mallado
axis([-3,3,-3,3]) %Establecemos los límites de los ejes x e y en la representación de la placa
hold on
[rho,theta]=meshgrid(u,v); % matrices en componentes polares de la placa sobre la que se trabaja
xx=rho.*cos(theta); % Parametrizamos
yy=rho.*sin(theta); % Parametrizamos
fx=(sin((pi*theta)/2))./(30*rho); % Componente x(rho,theta) de mi campo vectorial
fy=0*yy; % Componente y(rho,theta) del campo vectorial
quiver(xx,yy,fx,fy) % Representamos el campo vectorial sobre la placa
axis([-3,3,-3,3]) % Establecemos los límites de los ejes x e y en la representación de la placa
view(2) % Visualización
xlabel('x'); %Nombramos el eje x
ylabel('y'); %Nombramos el eje y
title('Mallado y campo de vectores u') %Damos un título al gráfico
6 Desplazamiento del sólido
La posición de nuestra placa circular sólida viene dada en cada punto por la ecuación:
r(x,y,t)= r0(x,y,t) + [math] \vec u[/math]. Siendo [math] \vec u[/math]. nuestro campo [math] \vec u(\rho,\theta)=\frac{\sin(\pi \theta/2)}{30\rho}\vec g_{\rho}. [/math]
h=0.1; % establecemos un paso de muestra
u=1:h:2; % valor u según el paso de muestreo en el intervalo [1,2]
v=0:h:2*pi+h; % valor v según el muestreo en el intervalo [0,2*pi]
[rho,theta]=meshgrid(u,v); % matrices en componentes polares de la placa sobre la que se trabaja
xx=rho.*cos(theta); % parametrizamos
yy=rho.*sin(theta);
fx=(sin((pi*theta)/2))./(30*rho); % componente x(rho,theta) de mi campo vectorial
fy=0*yy; % componente y(rho,theta) del campo vectorial
quiver(xx,yy,fx,fy) % representamos el campo vectorial sobre la placa
axis([-3,3,-3,3]) % establecemos los límites de los ejes x e y en la representación de la placa
view % visualización
La figura muestra la visualización tras aplicar el campo a los puntos. En ella se aprecia la dirección que cada uno de los puntos de nuestra placa sufre frente a su posición inicial.
Una vez aplicado el desplazamiento podremos comparar posiciones iniciales y finales.
h=0.1 % paso de muestreo
u=1:h:2; % valor u según el muestreo en el intervalo [1,2]
v=0:h:2*pi+h; % valor v según el muestreo en el intervalo [0,2*pi]
[rho,theta]=meshgrid(u,v); % matrices de las componentes polares de la placa
xx=rho.*cos(theta); % parametrización
yy=rho.*sin(theta);
subplot(1,2,1)
mesh(xx,yy,0*yy)
axis([-4,4,-4,4])
subplot(1,2,2)
fx=(sin((pi*theta)/2))./(30*rho);
fy=0*yy;
mesh(xx+fx,yy+fy,yy*0)
axis([-4,4,-4,4])
Como puede apreciarse, la variación es mínima. Ésto es debido a la ecuación del campo aplicada. En el punto de mayor deformación es a penas 1/30 para el caso del círculo interior de radio 1, y 1/60 en el caso del exterior por ser de radio 2. (Siendo por tanto evidente que el anillo se deforma más cuanto más cerca nos encontremos de su centro).
Sólo con una vista en detalle podremos observar que el desplazamiento sufrido rompe el anillo ya que desplaza la posición de la variable theta=0º, frente a theta= 360º.
7 Divergencia del Campo
La divergencia es una medida del cambio de volumen local debido al desplazamiento:
[math]\bigtriangledown{}\cdot{}\overrightarrow{u}
=\frac{1}{\rho}[(\frac{\partial}{\partial\rho})(\rho\cdot{}\frac{\sin(\frac{\pi\theta}{2})}{30\rho})] =0[/math]
El volumen de control se comprime, expande o queda igual. Si [math]\bigtriangledown{}\cdot{}\overrightarrow{u} =0[/math] para todo punto del dominio el campo se llama incompresible.
Al ser la divergencia 0 , todos los puntos del solido tienen la misma y el cambio de volumen experimentado por el sólido es nulo, ya que experimenta cambios transversales que no alteran su volumen.
8 Rotacional del Campo
El rotacional es un operador vectorial que muestra la tendencia de un campo vectorial a inducir rotación alrededor de un punto.
En nuestro caso el rotacional es el siguiente:
[math]\overrightarrow{\bigtriangledown{}}\times{}\overrightarrow{u}
=\frac{1}{\rho}\begin{vmatrix} \overrightarrow{g_\rho} & \overrightarrow{g_\theta} & \overrightarrow{g_z} \\ \frac{\partial}{\partial\rho} & \frac{\partial}{\partial\theta} & \frac{\partial}{\partial z} \\ \frac{sin\frac{\pi\theta}{2}}{30\rho} & 0 & 0 \end{vmatrix}=-(\frac{\pi}{60\rho^2})(cos\frac{\pi\theta}{2})\overrightarrow{g_z}[/math]
Y el módulo del rotacional es:
[math]|\overrightarrow{\bigtriangledown{}}\times{}\overrightarrow{u}|
=(\frac{\pi}{60\rho^2})(cos\frac{\pi\theta}{2})[/math]
h=0.1; % paso de muestreo
u=1:h:2; % valor u según el muestreo en el intervalo [1,2]
v=0:h:2*pi+h; % valor v según el muestreo en el intervalo [0,2*pi]
[rho,theta]=meshgrid(u,v); % matrices de las componentes polares de la placa
xx=rho.*cos(theta); % parametrización
yy=rho.*sin(theta);
f=abs((pi./((rho.^2)*60)).*(cos((pi.*theta)./2))); % campo escalar del rotacional
subplot(1,2,1)
surf(xx,yy,f) % representación del campo escalar sobre la placa
axis([-3,3,-3,3]) % límites de los ejes x e y en la representación de la placa
view % visualización del rotacional sobre la placa
subplot(1,2,2)
surf(xx,yy,f) % representación del campo escalar sobre la placa
axis([-3,3,-3,3]) % límites de los ejes x e y en la representación de la placa
view % visualización del rotacional sobre la placa
En la imagen de la representación del rotacional podemos observar que los puntos que tienen mayor rotacional son los más cercanos al radio menor de la placa ya que el rotacional depende de ρ en el denominador, y por tanto cuanto menor es ρ mayor es el rotacional.
9 Tensiones
En un medio elástico lineal, isótropo y homogéneo los desplazamientos permiten escribir el tensor de tensiones σij a través de la fórmula: σij = λ [math]\bigtriangledown{}\cdot{}\overrightarrow{u}[/math]δij + 2 μ εij. Donde εij es la parte simétrica del vector [math]\bar{u}[/math]. Tomamos λ=μ=1 y los vectores [math]\bar{i}[/math] = [math]\bar{g}[/math]ρ y [math]\bar{j}[/math] = [math]\bar{g}[/math]θ/ρ. En este apartado se nos pide dibujar las tensiones normales en la dirección que marca el eje [math]\bar{i}[/math], es decir [math]\bar{g}[/math]ρ · σ · [math]\bar{g}[/math]ρ = [math]\frac{-2sen(\pi \theta/2)}{30\rho^2} [/math] y las tensiones normales en la dirección que marca el eje [math]\bar{j}[/math], es decir [math]\bar{g}[/math]θ/ρ · σ · [math]\bar{g}[/math]θ/ρ = [math] \frac{\sin(\pi\theta/2)}{30\rho^2}[/math].
% Tensiones normales en la dirección que marca el eje [math]\bar{i}[/math].
h=0.1; %Establecemos paso de muestra
u=1:h:2; %Valor de u según el muestreo del intervalo
v=0:h:2*pi+h; %Valor de v según el muestreo del intervalo
[uu,vv]=meshgrid(u,v); %Matrices de la placa en coordenadas polares
figure(1) %Abrimos la ventana donde dibujaremos el mallado
xx=uu.*cos(vv); %Parametrización
yy=uu.*sin(vv); %Parametrización
f=((-2.*sin(pi*vv./2))./(30*(xx.^2+yy.^2))); % Campo [math]\bar{g}[/math]ρ
surf(xx,yy,f) % Dibujamos la malla
axis([-2,2,-2,2]) % seleccionamos región
axis equal
colorbar
title('Tensión normal en la dirección gsubρ')
view(1)
% Tensiones normales en la dirección que marca el eje [math]\bar{j}[/math].
h=0.1; %Establecemos paso de muestra
u=1:h:2; %Valor de u según el muestreo del intervalo
v=0:h:2*pi+h; %Valor de v según el muestreo del intervalo
[uu,vv]=meshgrid(u,v); %Matrices de la placa en coordenadas polares
figure(1) %Abrimos la ventana donde dibujaremos el mallado
xx=uu.*cos(vv); %Parametrización
yy=uu.*sin(vv); %Parametrización
subplot(1,2,2)
f=((sin(pi*vv./2))./(30*(xx.^2+yy.^2))); % Campo [math]\bar{g}[/math]θ/ρ
surf(xx,yy,f) % Dibujamos la malla
axis([-2,2,-2,2]) % Seleccionamos la región
axis equal
colorbar
title('Tensión normal en la dirección gsubθ/ρ')
view(2)
10 Tensiones tangenciales
Para calcular las tensiones tangenciales con respecto al plano ortogonal a [math]\bar{i}[/math] = [math]\bar{g}[/math]ρ utilizaremos el tensor de teniones σij y εij del apartado anterior. Siendo | σ · [math]\bar{g}[/math]ρ - ([math]\bar{g}[/math]ρ · σ · [math]\bar{g}[/math]ρ)· [math]\bar{g}[/math]ρ| = [math]\frac{πcos(\pi \theta/2)}{60\rho^3} [/math].
f=((pi.*cos(pi*vv./2))./(30*(uu.^3))); % tensiones tangenciales con respecto al plano ortogonal a
surf(xx,yy,f) % Dibujamos la malla
axis([-2,2,-2,2]) % seleccionamos región
axis equal
colorbar
title('Tensiones tangenciales con respecto al plano ortogonal a g¯ρ')
view(1)
Las tensiones tangenciales con respecto al plano ortogonal a [math]\bar{j}[/math] = [math]\bar{g}[/math]θ/ρ se calculan de la misma manera siendo | σ · [math]\bar{g}[/math]θ/ρ - ([math]\bar{g}[/math]θ/ρ · σ · [math]\bar{g}[/math]θ/ρ)· [math]\bar{g}[/math]θ/ρ| = [math]\frac{-πsen(\pi \theta/2)}{30\rho^5} [/math].
f=((sin(pi*vv./2))./(30*(uu.^5))); % tangenciales con respecto al plano ortogonal a g¯θ/ρ
surf(xx,yy,f) % Dibujamos la malla
axis([-2,2,-2,2]) % seleccionamos región
axis equal
colorbar
title('Tensiones tangenciales con respecto al plano ortogonal a g¯θ/ρ')
view(1)
La deformación es directamente proporcional a la tensión siendo la constante de young dicho grado de proporcionalidad. En comparación con la malla sin exposición a tensiones hay mayor deformación en la zona interna del anillo en los dos casos ya que en esa zona la tensión también es mayor.
