Visualización de campos escalares y vectoriales en elasticidad G13(A)
Contenido
- 1 Enunciado
- 2 Visualización de los campos en elasticidad
- 2.1 Mallado del sólido elástico
- 2.2 Representación de la temperatura
- 2.3 Representación del gradiente de la temperatura y de sus curvas de nivel
- 2.4 Representación del campo vectorial
- 2.5 Desplazamiento del mallado producido por el campo vectorial anterior
- 2.6 Divergencia del campo vectorial aplicado al sólido elástico
- 2.7 Rotacional del campo vectorial aplicado al sólido elástico
- 2.8 Tensiones normales a los ejes coordenados
- 2.9 Tensiones tangenciales respecto al plano ortogonal a \(\vec i\)
- 2.10 Tensiones tangenciales respecto al plano ortogonal a \(\vec j\)
- 2.11 Tensión de Von Mises
- 2.12 Masa total
1 Enunciado
Consideramos una placa rectangular plana (en dimensión 2) que ocupa la región \( [-1/2,1/2] \times [0,4]\). 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(x,y)\). De esta forma, si definimos \(r_0(x,y)\) el vector de posición de los puntos de la placa en reposo, la posición de cada punto \((x,y)\) viene dada por\[ \vec r (x,y)= \vec r_{0}(x,y)+\vec u(x,y). \] Vamos a suponer que la fuerza aplicada sobre la placa ha provocado un desplazamiento de los puntos de la misma dado por el vector de desplazamientos \[ \vec u(x,y) = \vec a (\vec b \cdot \vec r_0)^2, \] Donde \[ \vec a, \vec b\] son vectores dados: [math] \vec a=\frac{\vec j}{10}, \vec b=\vec j. [/math]
2 Visualización de los campos en elasticidad
En este trabajo vamos a observar el efecto que tienen los campos escalares y vectoriales sobre un sólido elástico, es decir, aquel sólido que es capaz de recuperar por completo su estado inicial al cesar las tensiones que se ejercen sobre él.
2.1 Mallado del sólido elástico
En este apartado nos piden dibujar un mallado que represente el sólido elástico del que nos habla el enunciado, tomando como paso de muestreo \(h=1/10\) para las variables x e y. El cógido y la gráfica son:
% Generamos los vectores
x=-0.5:0.1:0.5;
y=0:0.1:4;
% Generamos la malla
[Mx,My]=meshgrid(x,y);
figure(1)
% Pintamos la malla
mesh(Mx,My,0*Mx)
% Limitamos una región en los ejes
axis([-2,2,-1,5])
% Vista superior
view(2)
2.2 Representación de la temperatura
En el sólido elástico tenemos un foco de calor situado en el origen de coordenadas (0,0) y que no depende del tiempo. La función distribución de la temperatura en el sólido viene dada por la función [math]T(x,y)=(8-y^2+2y)e^{-x^2}[/math], en coordenadas cartesianas.
El gráfico de la temperatura muestra una gama de colores que van desde el rojo oscuro intenso al azul eléctrico, siendo el intervalo de temperaturas de la más elevada a la más baja respectivamente, comprobando que el centro de coordenadas tiene un color rojo intenso correspondiente a la mayor concentración en éste.
El código en MATLAB y los gráficos de la temperatura son los siguientes:
% generamos los vectores
x=-0.5:0.1:0.5;
y=0:0.1:4;
%hacemos la malla
[Mx,My]=meshgrid(x,y);
%escribimos la función temperatura
T=inline('(8 - y.^2 + 2*y).*exp(-(x).^2)','x','y');
%aplicamos la funcion a las matrices de la malla
TT=T(Mx,My);
%representamos la funcion
p=contour(Mx,My,TT,60);
axis([-1,1,-1,5])
%El maximo esta en el punto (0,1)
2.3 Representación del gradiente de la temperatura y de sus curvas de nivel
Dada la función de temperaturas del apartado anterior, vamos a calcular el gradiente para observar que es perpendicular a las curvas de nivel de dicha temperatura. El gradiente nos indica la menor distancia entre puntos de diferente temperatura, direccionado siempre hacia el foco de calor. Si nos movemos en la dirección del gradiente, aproximándonos al foco, la temperatura aumentará, ocurriendo lo contrario si nos alejamos de éste. Cuando intercalamos las curvas de nivel con el gradiente de la temperatura se puede observar que ambos son perpendiculares:
% generamos los vectores
u=-0.5:0.1:0.5; % generamos los vectores
v= 0:0.1:4;
[uu,vv]=meshgrid(u,v); %hacemos la malla
T=((8 - vv.^2 + 2.*vv).*exp(-(uu).^2));
%definimos las componentes del gradiente de la temperatura
figure(1)
Tx=(-2*(uu).*(8 - vv.^2 + 2*vv).*exp(-(uu).^2));
Ty=(-2*vv+2).*exp(-(uu).^2);
%representamos el campo vectorial generado por el gradiente
quiver(uu,vv,Tx,Ty)
hold on %unimos en unos ejes el gradiente con las lineas de nivel
contour(uu,vv,T) %dibujamos las lineas de nivel
%dibujamos un recuadro que limite la represenatacion
plot(u,u-u,'k','linewidth',1);
plot(u,6+u-u,'k','linewidth',1);
plot(-0.5+v-v,v,'k','linewidth',1);
plot(0.5+v-v,v,'k','linewidth',1);
axis([-1,1,-0.3,4.3]) %limitamos una region sobre los ejes
view(2)
%observamos que los vectores del gradiente son ortogonales a las
%curvas de nivel, ya que son vectores orientados al origen
%de las circunferencias concentricas que representan las
%estas curvas
2.4 Representación del campo vectorial
Ahora consideramos un campo vectorial definido como [math] u(x,y)=\displaystyle\frac{j}{10}(j·r_o)^2[/math],y desarrollando [math]r_o = xi + yj + zk[/math] al multiplicar escalarmente por j te queda simplificado el campo a:[math] u(x,y)=\displaystyle\frac{j}{10}(j·yj)^2[/math] equivalente a usar: [math] u(x,y)=\displaystyle\frac{y^2}{10}j[/math] que aplicado sobre el mallado producirá un desplazamiento de los puntos de éste. El campo vectorial, programado en MATLAB y dibujado, será:
%Intervalos de trabajo
x=-0.5:0.1:0.5;
y=0:0.1:4;
%Mallado
[X,Y]=meshgrid(x,y);
%Componentes del campo vectorial
fx=X.*0;
fy=(1/5).*Y;
%Dibujo del campo vectorial
quiver(X,Y,fx,fy)
axis([-1,1,-1,5.5])
view(2)
2.5 Desplazamiento del mallado producido por el campo vectorial anterior
Si aplicamos el campo vectorial \(\vec u\) al mallado que estamos estudiando, tenemos que los puntos de éste sufren un desplazamiento vertical, debido a que el campo vectorial solo tiene componente en el sentido del eje OY, es decir, solo tiene componente \(\vec j\). (CAMPO [math] u(x,y)=\displaystyle\frac{y^2}{10}j[/math]). Además, también dependerá de una única coordenada: y. Entonces, tendremos que el mallado sufre un desplazamiento de la siguiente forma:
%Intervalos del cuadrado
xc=-0.5:0.1:0.5;
yc=0:0.1:4;
%Mallado del cuadrado
[XC,YC]=meshgrid(xc,yc);
%Ecuaciones de los puntos desplazados
X=XC;
Y=YC+(1/5).*Y;
%Dibujo de los cuadrados
figure(1)
%Cuadrado sin desplazamiento
subplot(1,2,1);
mesh(XC,YC,0*XC)
axis([-1.5,2,-1,5])
%Cuadrado con desplazamiento
subplot(1,2,2)
mesh(X,Y,0*X);
axis([-1.5,1.5,-1,5]);
view(2)
2.6 Divergencia del campo vectorial aplicado al sólido elástico
A la hora de calcular la divergencia, estamos calculando el cambio de volumen local del sólido elástico cuando le aplicamos el campo vectorial \(\vec u\) en ese punto. Si ésta es positiva, tendremos una fuente (aumento de volumen) y si es negativa será un sumidero (disminución del volumen). El campo vectorial solo depende de la componente \(\vec j\) y de la coordenada y, al obtener la divergencia tenemos que la derivada parcial es [math]\displaystyle\frac{y}{5}j[/math], sacando como conclusión que el incremento de volumen del sólido elástico depende de y, es decir, que solo sufre desplazamientos verticales en el sentido de \(\vec j\), y si existe un incremento de volumen al ser la divergencia variable. Demostrado:
%Intervalos de trabajo
x=-0.5:0.1:0.5;
y=0:0.1:4;
%Mallado
[X,Y]=meshgrid(x,y);
%Componentes del campo vectorial
u=0.*X
v=(1/5).*Y;
%Divergencia
div=divergence(X,Y,u,v);
mesh(X,Y,div)
axis([-1,1,-0.5,5.5])
view(2)
2.7 Rotacional del campo vectorial aplicado al sólido elástico
Para hallar los puntos que sufren mayor rotacional debemos estudiar el módulo del producto \(\nabla \times \vec u\). El rotacional representa la tendencia de un campo vectorial a rotar alrededor de un punto, siendo mayor cuanto mayores son las derivadas en él. A la hora de calcularlo, es necesario utilizar los ejes coordenados, las componentes del campo vectorial y las derivadas parciales respecto a dichos ejes.
El rotacional en coordenadas cartesianas es:
[math]|\nabla\times\vec u|=\begin{vmatrix} \vec{i} & \vec{j} & \vec{k} \\ \frac{ \partial}{\partial x} & \frac{\partial }{\partial y} & \frac{\partial }{\partial z}\\ u_1 & u_2 & u_3 \end{vmatrix}= \begin{vmatrix} \vec{i} & \vec{j} & \vec{k} \\ \frac{ \partial}{\partial x} & \frac{\partial }{\partial y} & \frac{\partial }{\partial z}\\ \ 0 & \frac{y^2}{10} & 0 \end{vmatrix}=0[/math]
Por tanto es irrotacional.
2.8 Tensiones normales a los ejes coordenados
Las tensiones sobre un sólido vienen definidas por un tensor que depende de los coeficientes de Lamé (que son dos constantes elásticas que caracterizan por completo el comportamiento elástico de un sólido isótropo para pequeñas deformaciones).
[math]σ=λ \nabla uδ_{ij}+ μ·\varepsilon_{ij}[/math]
A continuación expondremos la gráfica en la que comparamos las tensiones normales con las del módulo de la divergencia (Figura 7).
Para hallar las tensiones normales en la dirección del eje i ,aplicamos [math] i·σ·i[/math] resultando [math]\displaystyle\frac{y}{5}[/math].
Para hallar las tensiones normales en la dirección del eje j aplicamos [math] j·σ·j[/math] también es [math]\displaystyle\frac{3y}{5}[/math].
x=[-1/2:0.1:1/2]; %Región de rango x entre los valores dados
y=[0,0.1,4]; %Región de rango y
[X,Y]= meshgrid(x,y); %Mallado de la placa
sigmax=Y/5 %Tensión eje x
sigmay=3.*Y./5 %Tensión eje y
subplot(1,2,1) %Creación de dos subventanas en el espacio de representación
mesh(sigmax) %Representación matriz tensiones en x
hold on
mesh(sigmay) %Representación matriz tensiones en y
hold off
subplot(1,2,2) %Cambio a la subventana 2
i=1; j=1;
Ux= [];
Uy=[];
for i=1:4
for j=1:1
Ux(i,j)=0; %Bucle para obtener matrices Ux, Uy con respecto a los valores i,j
Uy(i,j)=(Y(i,j).^2)/10;
end
end
div= divergence(X,Y,Ux,Uy) %divergencia del campo vectorial u de desplazamiento
mesh(div) %Representación de la divergencia
axis([-0.5,0.5,0,4]) %Caracterización de los ejes
%ojo hemos ampliado intervalos
2.9 Tensiones tangenciales respecto al plano ortogonal a \(\vec i\)
Las tensiones tangenciales son aquellas tangentes al plano de corte, y vienen definidas por la ecuación \(|\sigma \cdot \vec j-(\vec j \cdot \sigma \cdot \vec j) \vec j|\) respecto al plano ortogonal de \(\vec j\). Como ya se ha demostrado en el apartado 2.8, el tensor de tensiones es cero, haciendo cero el segundo sumando.
[math]|\left( \begin{array}{ll} \frac{y}{5}\ & 0 & 0 \\ 0 & \frac{3y}{5} & 0 \\ 0 & 0 & \frac{y}{5} \end{array} \right)\cdot\left( \begin{array}{l} 1 & 0 & 0 \end{array} \right)-( \left( \begin{array}{l} 1 \\ 0 \\ 0 \end{array} \right)\cdot\left( \begin{array}{ll} \frac{y}{5}\ & 0 & 0 \\ 0 & \frac{3y}{5} & 0 \\ 0 & 0 & \frac{y}{5} \end{array} \right)\cdot\left( \begin{array}{l} 1 & 0 & 0 \end{array} \right))\cdot\left( \begin{array}{l} 1 & 0 & 0 \end{array} \right)|=0 [/math]
2.10 Tensiones tangenciales respecto al plano ortogonal a \(\vec j\)
Las tensiones tangenciales son aquellas tangentes al plano de corte, y vienen definidas por la ecuación \(|\sigma \cdot \vec i-(\vec i \cdot \sigma \cdot \vec i) \vec i|\) respecto al plano ortogonal de \(\vec i\). Como ya se ha demostrado en el apartado anterior, el tensor de tensiones es cero, haciendo cero el segundo sumando.
[math]|\left( \begin{array}{ll} \frac{y}{5}\ & 0 & 0 \\ 0 & \frac{3y}{5} & 0 \\ 0 & 0 & \frac{y}{5} \end{array} \right)\cdot\left( \begin{array}{l} 0 & 1 & 0 \end{array} \right)-( \left( \begin{array}{l} 0 \\ 1 \\ 0 \end{array} \right)\cdot\left( \begin{array}{ll} \frac{y}{5}\ & 0 & 0 \\ 0 & \frac{3y}{5} & 0 \\ 0 & 0 & \frac{y}{5} \end{array} \right)\cdot\left( \begin{array}{l} 0 & 1 & 0 \end{array} \right))\cdot\left( \begin{array}{l} 0 & 1 & 0 \end{array} \right)|=0
[/math]
2.11 Tensión de Von Mises
La tensión de Von Mises se define como: [math] σ_{VM}=\sqrt[]{(σ_1-σ_2)^2+(σ_2-σ_3)^2+(σ_3-σ_1)^2}[/math]
donde σ1, σ2 y σ3 son los autovalores de σ (también conocidos como tensiones principales). Se trata de una magnitud escalar usada para indicar cuando un material se comporta de manera plástica, es decir cuando empieza a sufrir deformaciones no recuperables.
sigmatotal=[]; %Vector vacío para la primera interacción.
for y=linspace(0,0.1,4); %Bucle para crear el vector sigmatotal
sigma=[y/5 0 0;0 3*y/5 0;0 0 y/5]; %Matriz sigma
E=eig(sigma); %Función para sacar los autovalores de la matriz sigma
sigmavm= sqrt(((E(1)-E(2))^2+(E(2)-E(3))^2+(E(3)-E(1))^2)/2); %Ecuación de la tensión de Von Mises
sigmatotal=[sigmatotal,sigmavm] %Creación del vector sigmatotal a partir de las tensiones de las interacciones anteriores
end
plot(sigmatotal) %Dibujo de la tensión
2.12 Masa total
%Con el sigueinte programa calculamos la integral doble cuyo dominio es la placa e integrando la funcion densidad
syms x y;
f=input('Introducir funcion a integrar= ');
F=inline(char(f));
a=input('desde (x): ');
b=input('hasta (x): ');
a1=input('desde (y): ');
b1=input('hasta (y): ');
F=int(int(f,x,a,b),y,a1,b1)
%Ejecutando el programa, introduciendo la funcion densidad y los limites de integracion y se tiene que:
Introducir funcion a integrar= x*y*log(x+2)
desde (x): -0.5
hasta (x): 0.5
desde (y): 0
hasta (y): 4
F =
log(14348907/30517578125) + 8
>> log(14348907/30517578125) + 8
ans =
0.3376
