Visualización de campos escalares y vectoriales en elasticidad (Grupo 5-A)
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Visualización de campos escalares y vectoriales en elasticidad (Grupo 5-A) |
| Asignatura | Teoría de Campos |
| Curso | 2016-17 |
| Autores | Antón Fernández Arriola, Bruno Campuzano, Antonio Manso |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
Contenido
- 1 Enunciado
- 2 Representación del sólido
- 3 Temperatura del sólido
- 4 Gradiente de temperatura
- 5 Cálculo del campo de desplazamientos
- 6 Representación del campo de vectores
- 7 Sólido antes y después del desplazamiento
- 8 Representación de [math] \nabla\cdot\vec{u} [/math]
- 9 Cálculo y representación de [math] \left | \nabla \times \vec{u} \right | [/math]
- 10 Representación de las tensiones normales en las direcciones [math] \vec{i} , \vec{j} , \vec{k} [/math]
- 11 Cálculo de tensiones tangenciales
- 12 Representación de la tensión de Von Mises
- 13 Dibujar el campo de fuerzas [math] \vec{F} [/math]
- 14 Masa total
1 Enunciado
En este trabajo vamos a analizar y visualizar campos escalares y vectoriales en el espacio en relación a su elasticidad.
Consideraremos una placa rectangular plana en dos dimensiones que ocupa la región (x, y) ∈ [−2, 2]×[0, 4].
Vamos a definir en ella dos cantidades físicas: la temperatura [math] T(\rho,\theta)= \rho^2\sin(\theta) [/math] y los desplazamientos [math]\vec{u}=xf(y)\vec{j} [/math] producidos por la acción de una fuerza determinada. Dado que trabajaremos en coordenadas cartesianas hemos pasado la temperatura de coordenadas polares a cartesianas
[math] T(\rho,\theta)= \rho^2\sin(\theta) [/math]
[math] \rho = \sqrt{x^2+y^2} [/math] ; [math] \sin\theta = \displaystyle\frac{y}{x^2+y^2} [/math]
[math] T(x,y) = y\sqrt{x^2+y^2} [/math]
Por otro lado, consideramos también la divergencia [math] \nabla\cdot\vec{u} = x/10 [/math], la densidad [math] d(x,y) = \displaystyle\frac{1}{x^2+y^2+2} [/math]
, que los puntos situados en el eje y = 0 no sufren desplazamiento en la
dirección [math] \vec{j} [/math], la fórmula para hallar la tensión a partir del desplazamiento: [math] \sigma = \lambda\nabla\cdot\vec{u}1 + 2\mu\epsilon [/math]
, la tensión de Von Mises [math] \sigma_{VM} = \sqrt{\displaystyle\frac{(\sigma_{1}-\sigma_{2})^2+(\sigma_{2}-\sigma_{3})^2+(\sigma_{1}-\sigma_{1})^2}{2}} [/math]
donde [math] \sigma_{1}, \sigma_{2} y \sigma_{3} [/math] son los autovalores de [math] \sigma [/math]
y el campo de fuerzas [math] \vec{F} = -\nabla \cdot \sigma = -\displaystyle\frac{{\partial\sigma_{ji}}}{{\partial x_{j}}}\vec{e}_i [/math]
2 Representación del sólido
En primer lugar dibujaremos un mallado representando todos los puntos del interior del sólido. Tomamos los ejes en el rectángulo [math] (x,y) \in -3≤x≤3;0≤y≤5 [/math], y utilizamos como paso de muestreo [math] h=1/10 [/math] para las variables x e y. El código es el siguiente:
x= (-2:0.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
n=size(x,2);
m=size(y,2);
Z= zeros(n,m); %Altura de los puntos del intervalo
surf(X,Y,Z); %Representación de la superficie
axis equal; %Caracterización de los ejes
axis([-3,3,0,5]); %Limitamos una región en los ejes
view(2)
3 Temperatura del sólido
En este apartado representamos las curvas de nivel de la temperatura definidas por la función de la misma, anteriormente pasada de polares a cartesianas [math] T(\rho,\theta)= \rho^2\sin(\theta) [/math] pasa a [math] T(x,y) = y\sqrt{x^2+y^2} [/math]
x= (-2:0.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
F=inline('y.*sqrt(x.^2+y.^2)','x','y'); %Función de temperatura
Z1=F(X,Y) %Elevación de los puntos de la función de temperatura
axis([-3,3,0,5]) %Caracterización de los ejes
subplot(1,3,1) %Creación de subventana en el espacio de representación
surf(X,Y,Z1) %Representación superficie en 3D
view([0,0,1]) %Vista desde arriba (paso a 2D)
subplot(1,3,2) %Cambio a subventana 2
mesh(X,Y,Z1) %Representación en mallado de la superficie
max(max(F(X,Y))) %Calculo del máximo de la matriz
subplot(1,3,3)
contour(X,Y,Z1)
view(2)
Podemos deducir a partir de los gráficos que la temperatura máxima se da en los extremos (-2,4) Y (2,4)
4 Gradiente de temperatura
Después de representar la temperatura, calcularemos su gradiente y lo representaremos como un campo vectorial
x= (-2:0.1:2); %Región de rango x entre los valores dados
y= (0:0.1:4); %Región de rango x entre los valores dados
[X,Y]= meshgrid(x,y) %Mallado de la placa
T=Y.*sqrt(X.^2+Y.^2); %Función de temperatura
[X,Y]= gradient(T,0.1,0.1); %Cálculo del gradiente
contour(x,y,T,100) %Representación de las curvas de nivel
hold on
quiver(x,y,X,Y); %Representación de los vectores del campo gradiente
axis equal
hold off
Colocamos los vectores encima de las curvas de nivel y comprobamos que son ortogonales
5 Cálculo del campo de desplazamientos
Tomamos las condiciones que nos dan en el enunciado
[math]\vec{u}=xf(y)\vec{j} [/math]
[math] \nabla\cdot\vec{u} = x/10 [/math]
Y que los puntos situados en el eje y = 0 no sufren desplazamiento en la dirección [math] \vec{j} [/math]
Sabemos que [math] \nabla\cdot\vec{u} = {\partial x_{1}}/{\partial x} + {\partial x_{2}}/{\partial y} + {\partial x_{3}}/{\partial z} [/math]
Por tanto, la divergencia de vector [math] \vec{u} [/math] será [math] \nabla\cdot\vec{u} = xf'(y) [/math]
Igualando esta con la divergencia del encunciado [math] xf'(y) = \displaystyle\frac{x}{10} [/math] deducimos que [math] f'(y) = \displaystyle\frac{1}{10} [/math] y finalmente la integral da como resultado [math] f(y) = \displaystyle\frac{y}{10} + C [/math] Tomando por último la condición según la cual en el eje [math] y = 0 [/math] la dirección [math] \vec{j} [/math] no sufre desplazamiento y por tanto, se anula, nos quedará C = 0 y por tanto
[math] \vec{u} = xf(y) = x\displaystyle\frac{y}{10} [/math]
6 Representación del campo de vectores
En este apartado se nos pide dibujar el campo de vectores en los puntos del mallado del sólido. El código para ello es:
x= (-2:0.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
Ux=zeros(size(X)); %solo se desplazan en el eje j
Uy=X.*Y./10;
quiver(X,Y,Ux,Uy); %Representación de los vectores de desplazamiento
axis equal
axis([-3,3,0,5]); %Limitamos una región en los ejes
view(2)
7 Sólido antes y después del desplazamiento
Vamos a representar el cambio en el sólido una vez tiene lugar el desplazamiento debido al campo de vectores [math] \vec{u} [/math]
x= (-2:0.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
n=size(x,2);
m=size(y,2);
subplot(1,2,1);
Z= zeros(n,m); %Altura de los puntos del intervalo
mesh(X,Y,Z); %Representación de la superficie
axis equal; %Caracterización de los ejes
axis([-3,3,0,5]); %Limitamos una región en los ejes
title('Antes')
view(2)
subplot(1,2,2);
Ux2=zeros(size(X)); %solo se desplazan en el eje j
Uy2=X.*Y./10;
Rx = X+ Ux2;
Ry = Y + Uy2;
mesh(Rx,Ry,0*Rx) %Representación de los vectores de desplazamiento
axis equal
axis([-3,3,0,5]); %Limitamos una región en los ejes
title('Después')
view(2)
8 Representación de [math] \nabla\cdot\vec{u} [/math]
x= (-2:0.1:2); %Región de rango x entre los valores dados
y= (0:0.1:4); %Región de rango x entre los valores dados
[X,Y]= meshgrid(x,y) %Mallado de la placa
Divergencia=X/10; %divergencia dada por el enunciado
surf(X,Y,Divergencia)
Se observa en la imagen que los valores máximos de la divergencia se dan en los puntos donde x = 2 y los mínimos allí donde x = -2. Es nula en x = 0
9 Cálculo y representación de [math] \left | \nabla \times \vec{u} \right | [/math]
Siguiendo la fórmula del rotacional nos queda que [math] \nabla \times \vec{u} = \displaystyle\frac{y}{10}\vec{k} [/math] Como nos piden el módulo, [math] \left | \nabla \times \vec{u} \right | = \sqrt{\displaystyle\frac{y^2}{100}} [/math] y por tanto,
[math] \left | \nabla \times \vec{u} \right | = \displaystyle\frac{y}{10} [/math]
El código queda
x=(-2:0.1:2); %Región de rango x entre los valores dados
y=(0:0.1:4); %Región de rango y
z=zeros(size(x));
[X,Y,Z]=meshgrid(x,y,z);
x2=X.*0;
y2=Y.*0;
z2=Y./10;
quiver3(X,Y,Z,x2,y2,z2)
Los puntos que sufren un mayor rotacional son aquellos en los que y = 4
10 Representación de las tensiones normales en las direcciones [math] \vec{i} , \vec{j} , \vec{k} [/math]
Definido [math] \epsilon (\vec{u}) = (\nabla \vec{u} + \nabla \vec{u}^{t} )/2 [/math] como la parte simétrica del tensor de deformaciones y conocido el tensor de tensiones [math] \sigma = \lambda\nabla\cdot\vec{u} 1 + 2\mu\epsilon [/math] calculamos la tensión en cada una de las direcciones analíticamente. En la dirección [math] \vec{k} [/math] concluimos que no hay tensiones. A continuación, las escribimos directamente en Matlab:
%Eje x
x=(-2:0.1:2);
y=(0:0.1:4);
[X,Y]=meshgrid(x,y);
Txx=X/10; %Tensiones en la dirección x
Tyx=Y*0;
subplot(1,2,1)
quiver(X,Y,Txx,Tyx)
axis equal
axis([-3,3,0,5]);
title('Tensiones normales eje i')
%Eje y
Txy=0*Y; %Tensiones en la dirección y
Tyy=3*X/10;
subplot(1,2,2)
quiver(X,Y,Txy,Tyy)
axis equal
axis([-3,3,0,5]);
title('Tensiones normales eje j')
11 Cálculo de tensiones tangenciales
Nos piden calcular [math] \left |\sigma\cdot\vec{i} - (\vec{i}\cdot\sigma\cdot\vec{i} )\vec{i}\right | [/math]
[math] (\vec{i}\cdot\sigma\cdot\vec{i} )\cdot\vec{i} = (x/10\cdot (1,0,0)) = \begin{pmatrix} x/10 \\ 0 \\ 0 \end{pmatrix} [/math]
[math] \sigma\cdot\vec{i} = \begin{pmatrix} x/10 & y/10 & 0 \\ y/10 & 3x/10 & 0 \\ 0 & 0 & x/10 \end{pmatrix}\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix} = \begin{pmatrix} x/10 \\ y/10 \\ 0 \end{pmatrix} [/math]
[math] \sigma\cdot\vec{i} - (\vec{i}\cdot\sigma\cdot\vec{i} )\vec{i} = \begin{pmatrix} 0 \\ -y/10 \\ 0 \end{pmatrix} [/math]
[math] \left | \sigma\cdot\vec{i} - (\vec{i}\cdot\sigma\cdot\vec{i} )\vec{i} \right | = y/10 [/math]
12 Representación de la tensión de Von Mises
Nos dan la fórmula [math] \sigma_{VM} = \sqrt{\displaystyle\frac{(\sigma_{1}-\sigma_{2})^2+(\sigma_{2}-\sigma_{3})^2+(\sigma_{1}-\sigma_{1})^2}{2}} [/math]
donde [math] \sigma_{1}, \sigma_{2} y \sigma_{3} [/math] son los autovalores de [math] \sigma [/math]
x=(-2:0.1:2);
y=(0:0.1:4);
[X,Y]=meshgrid(x,y);
sigma=[]; %Vector sigma vacío.
for i=1:length(x); %Bucle para crear sigmatotal
for j= 1: length(y);
sigma=[X(i,j)./10 Y(i,j)./10 0;Y(i,j)./10 (3.*X(i,j))./10 0;0 0 X(i,j)./10]; %Matriz sigma
EI=eig(sigma); %Función para sacar los autovalores de la matriz sigma
sigmavm= sqrt(((EI(1)-EI(2))^2+(EI(2)-EI(3))^2+(EI(3)-EI(1))^2)/2); %Ecuación de la tensión de Von Mises
Z(i,j)=sigmavm;
end
end
surf(X,Y,Z) %Dibujo de la tensión
13 Dibujar el campo de fuerzas [math] \vec{F} [/math]
Dada la fórmula del enunciado [math] \vec{F} = -\nabla \cdot \sigma = -\displaystyle\frac{{\partial\sigma_{ji}}}{{\partial x_{j}}}\vec{e}_i [/math] que nos indica el campo de fuerzas que actúan sobre la placa y que provocan el desplazamiento:
x=(-2:0.1:2); %Region de rango x entre los valores dados
y=(0:0.1:4); %Region y
z=zeros(size(x));
[X,Y,Z]=meshgrid(x,y,z);
X1=-X/10-Y/10;
Y1=X*0;
Z1=Z*0;
quiver3(X,Y,Z,X1,Y1,Z1)
14 Masa total
Para calcular la masa total de la placa calculamos la integral doble de la función de densidad [math] d(x,y) = \displaystyle\frac{1}{x^2+y^2+2} [/math] tomando como extremos los de la placa
x1=2;
x2=-2;
y1=4;
y2=0;
L1=x1-x2; %lado1 de la placa
L2=y1-y2; %lado2
AREA=L1*L2;
syms X Y; %funcion para nombrar las variables de la integral
D=inline( 1./(X.^2+Y.^2+2)); %funcion de densidad
f=int(int(D,X,x2,x1),Y,y2,y1); %Integral de la función de densidad
Masa=D*AREA
La masa total es M = 1088
