Campo de desplazamiento de una placa semicircular A4
Consideramos una placa plana que ocupa el anillo circular centrado en el origen y comprendido entre los radios 1 y 2 en el semiplano y≥0. En ella vamos a suponer que tenemos definidas dos cantidades físicas: la temperatura \(T(x,y)\), que depende de las dos variables espaciales \((x,y)\), y los desplazamientos \(\vec u\) producidos por la acción de una fuerza determinada.
A partir de esto, estudiaremos el efecto de la temperatura dada por \(T(x,y)\) para después plantear la acción de las vibraciones \(\vec u(\rho,\theta)\) y las tensiones que éstas producen.
Contenido
- 1 Visualización de la placa
- 2 Distribución de temperaturas del sólido
- 3 Estudio del gradiente de temperaturas
- 4 Campo vectorial en los puntos del sólido
- 5 Efecto de los desplazamientos
- 6 Estudio de la divergencia
- 7 Estudio del rotacional
- 8 Tensiones normales
- 9 Tensiones tangenciales
- 10 Tensión de Von Mises
- 11 Masa total de la placa
1 Visualización de la placa
La placa es un semianillo comprendido entre los radios 1 y 2.
h= 0.1; %paso de muestreo
%usamos coordenadas polares
r= 1:h:2;
tetha= 0:h:pi;
[rr,tt]= meshgrid(r,tetha); %mallado
%parametrizamos en cartesianas
x=rr.*cos(tt);
y=rr.*sin(tt);
clf %borramos las posibles gráficas que hubiera
mesh(x,y,0*x) %hacemos la visualización de la placa
view(2) %vemos la placa en 2D
axis ([-3,3,-1,3])
2 Distribución de temperaturas del sólido
La distribución de temperaturas en cada punto del sólido viene dada por el campo escalar:
\(T(x,y)=log({y+2})\)
h= 0.1; %paso de muestreo
%usamos coordenadas polares
r= 1:h:2;
tetha= 0:h:pi;
[rr,tt]= meshgrid(r,tetha); %mallado
%parametrizamos en cartesianas
x=rr.*cos(tt);
y=rr.*sin(tt);
clf
%campo temeperatura (en cartesianas)
T=log(y+2); %campo escalar de temperatura
subplot(1,2,1) %Dividimos la pantalla en dos
surf(x,y,T) %representamos el campo escalar de temperaturas
view(2)
axis ([-3,3,-1,3])
colorbar %mostramos la escala
subplot(1,2,2) %escribimos en la segunda pantalla
contour(x,y,T,40) %lineas de nivel
colorbar %mostramos las escala
axis ([-3,3,-1,3])Según podemos observar en el gráfico obtenido, el campo de temperaturas aumenta en dirección de y, y es máximo en el punto y=2
3 Estudio del gradiente de temperaturas
Al calcular el gradiente de nuestro campo de temperaturas obtenemos el siguiente campo vectorial:
grad T = \(1./(y+2)\vec j\)
Este campo indica la dirección en la cual la temperatura del sólido aumenta a mayor velocidad, y su módulo nos indica dicha velocidad.
h= 0.1; %paso de muestreo
%usamos coordenadas polares
r= 1:h:2;
tetha= 0:h:pi;
[rr,tt]= meshgrid(r,tetha); %mallado
%parametrizamos en cartesianas
x=rr.*cos(tt);
y=rr.*sin(tt);
% lineas nivel campo temeperatura (en cartesianas)
T=log(y+2); %campo escalar de temperatura
subplot(1,2,1) %dividimos la pantalla en dos
contour(x,y,T,40) %lineas de nivel
colorbar % mostramos las escala
axis ([-3,3,-1,3])
%gradiente temperatura (en cartesianas)
tx=inline('0','x','y'); %derivada parcial respecto de x
ty=inline('1./(y+2)','x','y'); %derivada parcial respecto de y
TX=tx(x,y);
TY=ty(x,y);
subplot(1,2,2)
h= quiver(x,y,TX,TY); %representamos el campo vectorial
axis ([-3,3,-1,3])
set(h,'maxheadsize',0.33) %cambiamos formato flechas
Como vemos en las dos gráficas, las curvas de nivel son perpendiculares al gradiente de T.
4 Campo vectorial en los puntos del sólido
A continuación vamos a observar la deformación de la placa producida por el siguiente campo de desplazamientos:
\(\vec u(\rho,\theta)=1/5*sin(2θ)\vec g_{\theta}\)
h= 0.1; %paso de muestreo
%usamos coordenadas polares
r= 1:h:2;
theta= 0:h:pi;
[rr,tt]= meshgrid(r,theta); %mallado
%parametrizamos en cartesianas
x=rr.*cos(tt);
y=rr.*sin(tt);
clf %Borramos las posibles gráficas anteriores
%escribimos el campo de desplazamientos en coordenadas cartesianas
%T=ui+vj
u=inline('-y/5.*sin(2.*atan(y./x))','x','y');
v=inline('x/5.*sin(2.*atan(y./x))','x','y');
%aplicamos el campo de desplazamientos a todos los puntos del mallado de la placa
U=u(x,y);
V=v(x,y);
quiver(x,y,U,V) %Hacemos la gráfica del gradiente
axis ([-3,3,-1,3]) %Ajustamos los ejes
El efecto del campo de desplazamientos es girar los puntos de la placa hacia la parte superior del semianillo. Su efecto es mayor en los puntos que están a mitad de altura (puntos próximos a θ=±π/4).
5 Efecto de los desplazamientos
En este apartado vamos a comparar los puntos de la placa antes y después de los desplazamientos para observar los efectos de estos.
h= 0.1; %paso de muestreo
%usamos coordenadas polares
r= 1:h:2;
tetha= 0:h:pi;
[rr,tt]= meshgrid(r,tetha); %mallado
%parametrizamos en cartesianas
x=rr.*cos(tt);
y=rr.*sin(tt);
clf
subplot(1,2,1) %sólido antes de los desplazamientos
mesh(x,y,0*x)
view(2)
axis ([-3,3,-1,3])
subplot(1,2,2) %sólido después de los desplazamientos
u=inline('-y/5.*sin(2.*atan(y./x))','x','y');
v=inline('x/5.*sin(2.*atan(y./x))','x','y');
U=u(x,y);
V=v(x,y);
X=x+U;
Y=y+V;
mesh(X,Y,0*X)
view(2)
axis ([-3,3,-1,3])Vemos como los puntos de los extremos están más separados y que los del centro están más juntos
6 Estudio de la divergencia
Aquí estudiaremos las modificaciones de volumen local producidas por el campo de desplazamientos.
El campo es:
\(\vec u(\rho,\theta)=1/5*sin(2*theta)\vec g_{\theta}\)
Como lo tenemos en cilindricas hacemos la divergencia en cilindricas y luego pasamos a cartesianas para poder representarlo con matlab.
h= 0.1; %paso de muestreo
%usamos coordenadas polares
r= 1:h:2;
tetha= 0:h:pi;
[rr,tt]= meshgrid(r,tetha); %mallado
%parametrizamos en cartesianas
x=rr.*cos(tt);
y=rr.*sin(tt);
clf
%Divergencia en cilíndricas= 2/5*cos(2*tetha)
%pasamos a cartesianas: Divergencia en cartesianas= 2/5*cos(2*atan(y/x))
div=inline('2.*cos(2.*atan(y./x))./5','x','y');
DIV=div(x,y);
surf(x,y,DIV)
colorbar
axis ([-3,3,-1,3])
view(3)
La divergencia es mayor en los extremos del semianillo. Según se observa en la gráfica, el valor de dicha divergencia en estos puntos es 0.4 (se ve en color rojo y más alto). No se aprecia mucho el cambio de volumen porque el campo no se caracteriza por estar expandiendose, si no por su rotación. El campo tiende a comprimir los puntos en el medio del anillo
7 Estudio del rotacional
El rotacional de un campo nos muestra como gira la superficie al aplicarle el campo u. En este apartado vamos a estudiar el giro del campo:
\(\vec u(\rho,\theta)=1/5*sin(2θ)\vec g_{\theta}\)
Al igual que en la divergencia, como el campo dado está en coordenadas cilíndricas, operaremos así y luego las pasaremos a cartesianas para hacer la gráfica en matlab.
h= 0.1; %paso de muestreo
%usamos coordenadas polares
r= 1:h:2;
tetha= 0:h:pi;
[rr,tt]= meshgrid(r,tetha); %mallado
%parametrizamos en cartesianas
x=rr.*cos(tt);
y=rr.*sin(tt);
%Dibujo
subplot(1,2,1);
rot=2/5*sin(2*tt);
surf(x,y,rot)
axis equal
view(2) %rotacional 2D
colorbar
subplot(1,2,2);
surf(x,y,rot) %rotacional 3D
colorbar
Maximo=max(max(rot))Cómo vemos en la gráfica, el rotacional es máximo para los puntos con θ=π/4
8 Tensiones normales
Si consideramos que estamos en un medio elástico lineal, isótropo y homogéneo, podemos definir el tensor de tensiones σ como:
σ(x,y)=λ∇.u.1 + 2με
siendo λ y μ los coeficientes de Lamé, que dependen de cada material. En nuestro caso λ y μ son iguales a 1. Además, ε es la parte simétrica del vector gradiente u, es decir, (∇u + ∇ut)/2.
Aunque los desplazamientos se realizan en el plano, las tensiones no tienen por qué ser planas, y se puede representar las tensiones en las tres direcciones del espacio. En nuestro caso, como estamos usando coordenadas cilíndricas, representamos las tensiones en las direcciones de la base ortogonal cilíndrica {\(\vec g_{ρ}\), \(\vec g_{\theta}\), \(\vec g_{z}\)}. Para ello, multiplicamos escalarmente la matriz sigma por el vector \(\vec g_{i}\) por ambos lados: tensión normales en dirección i = \(\vec g_{i}\) * σ * \(\vec g_{i}\)
a) dirección \(\vec g_{\rho}\): (1,0,0)*σ*(1,0,0) = 2/5*cos(2θ)
b) dirección \(\vec g_{\theta}\)/ρ: (0,1,0)*σ*(0,1,0) = 2/5*cos(2θ)*(1+2/ρ^2)
c) dirección \(\vec g_{z}\): (0,0,1)*σ*(0,0,1) = 2/5*cos(2θ)
Representamos en 2D, y en 3D para mejor comprensión del efecto de las tensiones.
2D
h=0.1;
r=1:h:2;
t=0:h:pi;
[rr,tt]=meshgrid(r,t);
clf
subplot(1,3,1);
xx=rr.*cos(tt);
yy=rr.*sin(tt);
a=2/5.*cos(2*tt); %elemento (1,1) de la matriz sigma
b=a.*(1+2./rr.^2);
surf(xx,yy,a)
axis equal
view(2)
colorbar
subplot(1,3,2); %elemento (2,2) de la matriz sigma
surf(xx,yy,b)
axis equal
view(2)
colorbar
subplot(1,3,3)
c=a; %elemento (3,3) de la matriz sigma
surf(xx,yy,c)
axis equal
view(2)
colorbar
3D
h=0.1;
r=1:h:2;
t=0:h:pi;
[rr,tt]=meshgrid(r,t);
clf
subplot(1,3,1);
xx=rr.*cos(tt);
yy=rr.*sin(tt);
a=2/5.*cos(2*tt); %elemento (1,1) de la matriz sigma
b=a.*(1+2./rr.^2);
surf(xx,yy,a)
colorbar
subplot(1,3,2); %elemento (2,2) de la matriz sigma
surf(xx,yy,b)
colorbar
subplot(1,3,3)
c=a; %elemento (3,3) de la matriz sigma
surf(xx,yy,c)
colorbar
9 Tensiones tangenciales
Las tensiones del apartado son σ.*\(\vec g_{\rho}\)+(\(\vec g_{\rho}\)*σ*\(\vec g_{\rho}\))\(\vec g_{\rho}\), y como \(\vec g_{\rho}\)tiene de coordenadas cilíndricas (1,0,0), esto da 1/ρ*(1/5*sen(2θ)-1)
Respecto al plano ortogonal a \(\vec g_{\rho}\)
h=0.1;
r=1:h:2;
t=0:h:pi;
[rr,tt]=meshgrid(r,t);
xx=rr.*cos(tt);
yy=rr.*sin(tt);
a=1./rr.*(1/5.*sin(2*tt)-1);
subplot(2,1,1)
surf(xx,yy,a)
colorbar
subplot(2,1,2)
surf(xx,yy,a)
colorbar
view(2)
Plano ortogonal a \(\vec g_{\theta}/ρ\)
h=0.1;
r=1:h:2;
t=0:h:pi;
[rr,tt]=meshgrid(r,t);
xx=rr.*cos(tt);
yy=rr.*sin(tt);
a=1./rr.*(1/5.*sin(2*tt)-1);
subplot(2,1,1)
surf(xx,yy,a)
colorbar
subplot(2,1,2)
surf(xx,yy,a)
colorbar
view(2)
Como se puede ver, ambas tensiones tangenciales son iguales. Esto se debe a que las tensiones de la primera parte del apartado dan una función de tensiones igual al elemento (2,1) de la matriz σ, y las tensiones del segundo apartado dan lugar a auna función de tensiones igual al elemento (1,2). Como la matriz σ es simétrica, estos dos elementos coinciden.
10 Tensión de Von Mises
La tensión de Von Mises se define como:
siendo \(\ σ_{i}\) con i=1,2,3 los autovalores de σ (también son conocidos como tensiones principales)
La tensión de von Mises es una magnitud escalar que es un indicador para saber cuando un material inicia un comportamiento plástico y no elástico puro
h=0.1;
r=[1:h:2];
t=[0:h:pi];
[rr,tt]=meshgrid(r,t);
%creación de la matriz sigma
s=zeros(3,3);
%creamos la matriz de von mises
vm=zeros(32,11); %tiene un tamaño 32x11 porque es el tamaño del mallado
for i=1:32*11
rrr=rr(i)';
ttt=tt(i)';
s(1,1)=(2/5).*cos(2.*ttt);
s(1,2)=(1./rrr).*(0.2.*sin(2*ttt)-1);
s(2,1)=s(1,2);
s(2,2)=s(1,1)+0.8.*cos((2*ttt)./(rrr.^2));
s(3,3)=s(1,1);
[v,d]=eig(s); %autovalores de la matriz sigma
vm(i)=sqrt(((d(1,1)-d(2,2))^2+(d(2,2)-d(3,3))^2+(d(3,3)-d(1,1))^2)/2); %fórmula de Von Mises
end
%cambiamos a cartesianas
xx=rr.*cos(tt);
yy=rr.*sin(tt);
surf(xx,yy,vm)
colorbar
view(2)
11 Masa total de la placa
La densidad de la placa viene dada en coordenadas cartesianas.
función de densidad = log(x^2+y^2+2)
Pero en nuestro caso es más conveniente pasar la densidad a coordenadas cilíndricas.
función de densidad = log(ρ^2+2)
Como no depende del ángulo θ, podemos transformar la integral doble en una integral simple, porque el dθ sale fuera de la integral:
Masa = π * integral simple que solo depende de ρ
Esta integral la realizamos utilizando métodos numéricos. En concreto, usamos la integración siguiendo el método del trapecio (que consiste en aproximar el área bajo la gráfica de la función, por trapecios muy pequeños de altura h)
Por ello primero creamos la función integralTrap que realiza la integral por el método del trapecio
%entradas a la función: f=función a integrar; a,b=límites de integración; m=número de subdivisiones (cuanto mayor sea, más preciso es el resultado)
%salida de la función: I=valor de la integral
function I=integralTrap(f,a,b,m)
h=(b-a)/m; %altura de los trapecios diferenciales
x=linspace(a,b,m+1);
I=0;
for i=1:m
A=(feval(f,x(i))+feval(f,x(i+1)))*h/2; %área del trapecio (evaluamos la función en los vértices superiores)
I=I+A;
end
Ahora ya podemos aplicar el método del trapecio
f=inline('r*log(r^2+2)','r') %creamos la función densidad en coordenadas cilíndricas
%los límites de integración son r=1, r=2
%introducimos 100 subdivisiones
masa=pi.*integralTrap(f,1,2,100)Que da un resultado de:
masa=6.9975
