Observación de campos escalares y vectoriales en elasticidad

De MateWiki
Revisión del 02:31 3 dic 2022 de Alberto.munozf (Discusión | contribuciones) (→‎Tensor de tensiones)

Saltar a: navegación, buscar
Trabajo realizado por estudiantes
Título Observación de campos escalares y vectoriales en elasticidad
Asignatura Teoría de Campos
Curso 2022-23
Autores Daniel Casas
Pablo Moreno
Alberto Muñoz
Alberto Núñez
Juan Utrilla
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

Consideramos una placa rectangular plana (en dimensión 2) que ocupa la región [math](x, y) ∈ [0, 10]×[0, 2][/math]. En ella vamos a suponer que tenemos definidas dos cantidades físicas: la temperatura [math]T(x, y)[/math], que viene dada por,

[math]T(x, y) = (x − 3)2 + (10(y − 1/2))2[/math]

y los desplazamientos [math]\vec{u}(x, y)[/math] producidos por la acción de una fuerza determinada. De esta forma, si definimos [math]\vec{r_{0}}(x, y)[/math] el vector de posición de los puntos de la placa antes de la deformación, la posición de cada punto [math](x, y)[/math] de la placa después de la deformación viene dada por

[math]\vec{r_{d}}(x, y)=\vec{r_{0}}(x, y)+\vec{u}(x, y).[/math]

Vamos a suponer que la fuerza aplicada sobre la placa ha provocado un desplazamiento ondulatorio de los puntos de la misma dado por el vector

[math]\vec{u}(x, y, t)=\vec{a}sin(k(\vec{d}·\vec{r_{0}}(x,y)-vt)),[/math]

donde [math]\vec{a}[/math] se conoce como amplitud, [math]k \gt 0[/math] es el número de onda, [math]\vec{d}[/math] es un vector unitario que marca la dirección de propagación y [math]v[/math] es la velocidad de propagación. La variable [math]t[/math] representa el tiempo que congelaremos en [math]t = 0[/math] en los primeros 10 apartados de este trabajo de manera que supondremos, para los primeros apartados,

[math]\vec{u}(x, y, t)=\vec{a}sin(k(\vec{d}·\vec{r_{0}}(x,y))).[/math]

Supondremos que se trata de una onda transversal en la que la dirección de propagación es ortogonal a la amplitud. Tomaremos en particular

[math]\vec{a}=2/5\vec{j}, k=1, \vec{d}=\vec{i}[/math]

1 Representación placa rectangular

Dibujar un mallado que represente los puntos interiores del sólido. Tomar los ejes (comando axis) en el rectángulo [math](x, y) ∈ [−0.5; 10.5] × [−0.5; 2.5][/math] y como paso de muestreo [math]h = 2/10[/math] para las variables [math]x[/math] e [math]y[/math].

Comenzamos el programa utilizando los programas de limpieza para borrar todos los anteriores. Creamos el paso de muestreo [math]h[/math] y definimos los parámetros de la placa [math]x[/math] e [math]y[/math]. Creamos el mallado con el comando meshgrid() y, a continuación, dibujamos el mallado con mesh() habiendo definido ya las matrices [math]X, Y, Z[/math] (siendo esta última 0). Finalmente, ajustamos los ejes a los límites establecidos, los nombramos y ponemos un título al gráfico.


Mallado
clc
clear
%Definimos las regiones
h=2/10; 
x=[0:h:10];       
y=[0:h:2];
%Creamos el mallado
[X,Y]=meshgrid(x,y); 
mesh(X,Y,0*Y)
%Establecemos los límites
axis([-0.5,10.5,-0.5,2.5]);
%Escribimos el título y nombramos los ejes
title('Mallado');
xlabel('X');
ylabel('Y');

view(2);


2 Curvas de nivel de la temperatura y punto en el que es máxima

Dibujar las curvas de nivel de la temperatura (comando contour) y decidir en qué punto la temperatura es máxima a partir de la gráfica.

Repetimos el proceso anterior para definir las matrices [math]X[/math] e [math]Y[/math], y definimos la función dada de la temperatura [math]T[/math] en función de esos parámetros. Representamos la gráfica mediante el comando [math]mesh()[/math] con los parámetros [math]X[/math], [math]Y[/math] y [math]T[/math], y utilizamos [math]contour()[/math] para representar las curvas de nivel que queramos. De nuevo ajustaremos los ejes y les pondremos nombre, incluyendo además una barra de color a un lado del gráfico que facilite la visualización de este, usando [math]colorbar[/math]. Observando el gráfico podemos decir que la temperatura máxima es 274ºC (x=10; y=2).

Curvas de nivel
clc
clear
%Definimos regiones
h = 2/10;
x=[0:h:10];          
y=[0:h:2];
%Creamos el mallado
[X,Y]= meshgrid(x,y);
%Función temperatura
T = (X-3).^2 + (10*(Y-1/2)).^2;
%Escribimos el título y ejes
title('Curvas de nivel de temperatura');
xlabel('X');
ylabel('Y');
%Establecemos el límite de los ejes
axis([-0.5,10.5,-0.5,2.5]);
%Representación de la temperatura y las curvas de nivel
subplot(1,2,1);
mesh(X,Y,T);
subplot(1,2,2);
contour(X,Y,T,50);
colorbar


3 Gradiente de la temperatura

Calcular [math]∇T[/math] y dibujarlo como campo vectorial. Observar gráficamente que [math]∇T[/math] es ortogonal a dichas curvas.

Realizaremos el calculo del gradiente:

[math]\frac{\partial T}{\partial x}\vec{i} + \frac{\partial T}{\partial y}\vec{j} = (2x+6)\vec{i}+(200y-100)\vec{j}[/math]

Para la realización del código de matlab, definiremos las mismas variables y regiones anteriores, y representaremos de nuevo las curvas de nivel de la temperatura. Crearemos las variables [math]dx[/math] y [math]dy[/math] en función de las matrices [math]X[/math] y [math]Y[/math]. Por último, para poder representarlo como campo vectorial, utilizaremos el comando [math]quiver[/math].

Gradiente
clc
clear
%Definimos regiones
h = 2/10;
x=[0:h:10];          
y=[0:h:2];
%Creamos el mallado
[X,Y]= meshgrid(x,y);
%Función temperatura
T = (X-3).^2 + (10*Y).^2;
contour(X,Y,T,30);
%Gradiente de T
dx = 2*X-6;
dy = 200*Y-100;
%Título y ejes
title('Gradiente de temperatura');
xlabel('X');
ylabel('Y');
%Representación de la temperatura y las curvas de nivel
hold on
quiver(x,y,dx,dy);
axis equal
colorbar


4 Campo de vectores de desplazamiento

Dibujar el campo de vectores en los puntos del mallado del sólido, en [math]t=0[/math].

Calcularemos [math]\vec{u}[/math], y de nuevo estableceremos la región y los parámetros anteriores. Definiremos el campo de desplazamientos con las variables [math]ux[/math] y [math]uy[/math] en función de las matrices [math]X[/math] e [math]Y[/math], y utilizaremos el comando [math]quiver()[/math] para representar el campo de vectores.


Campo de vectores
clc
clear
%Definimos regiones
h=1/10; 
x=[0:h:10];          
y=[0:h:2];
%Creación del mallado
[X,Y]=meshgrid(x,y);
%Definimos el campo de desplazamientos
ux = X*0;
uy = (2/5).*sin(X);
%Pintamos el campo de vectores en los puntos del mallado
quiver(x,y,ux,uy)
%Establecemos límites
axis([-0.5,10.5,-0.5,2.5]);
axis equal
%Escribimos el título del gráfico
title('Campo de desplazamientos');
%Escribimos los nombres de los ejes
xlabel('X');
ylabel('Y');
view(2);


5 Desplazamiento del sólido

Dibujar el sólido antes y después del desplazamiento dado por el campo de vectores [math]\vec{u}[/math] (en [math]t=0[/math]). Dibujar ambos en la misma figura usando el comando [math]subplot[/math].

Definiremos las regiones y los parámetros del sólido, distinguiendo ahora los datos antes y después del desplazamiento que sufre el sólido. Representaremos ambas situaciones, y por último, representaremos una mezcla de las dos, en la que se puedan comparar.

Desplazamiento
clc
clear

%Definición de regiones
h=2/10; 
x=[0:h:10];          
y=[0:h:2];

%Matriz de x e y
[Mx,My]=meshgrid(x,y);

%Vector desplazamiento
Ux= 0;
Uy=(2/5).*sin(Mx);

figure

%Antes del desplazamiento
subplot(2,2,1)
mesh(Mx,My,0*My)
axis([-0.5,10.5,-0.5,2.5]);
axis equal
view(2)
title('Situación inicial');
xlabel('X');
ylabel('Y');

%Después del desplazamiento
subplot(2,2,2)
mesh(Mx+Ux,My+Uy,0*Uy)
axis([-0.5,10.5,-0.5,2.5]);
axis equal
view(2)
title('Situación final');
xlabel('X');
ylabel('Y');

%Comparación
subplot(2,2,3)
plot3(Mx,My,0*My)
hold on
plot3(Mx+Ux,My+Uy,0*Uy)
hold off
axis([-0.5,10.5,-0.5,2.5]);
axis equal
view(2)
title('Comparación');
xlabel('X');
ylabel('Y')


6 Divergencia del campo vectorial

Dibujar [math]∇·\vec{u}[/math] en [math]t=0[/math]. Determinar analíticamente los puntos en los que la divergencia de [math]\vec{u}[/math] es máxima, mínima y nula. La divergencia es una medida del cambio de volumen local debido al desplazamiento. ¿Se puede apreciar esto en la gráfica?

Conocido [math]\vec{u}[/math], calcularemos la divergencia:

[math]∇·\vec{u} = \frac{∂}{∂x}(0)+\frac{∂}{∂y}(2/5senx) = 0[/math]

Al ser el resultado observaremos un color uniforme en la gráfica, sin ninguna variación en el volumen local. Para el desarrollo del código, estableceremos las regiones y los parámetros utilizados anteriormente, y definiremos la variable de la divergencia [math]div[/math], la cual será nula, y representaremos el gráfico.

Divergencia
clc
clear

%Definición de regiones
h=2/10; 
x=[0:h:10];          
y=[0:h:2];

%Matriz de x e y
[Mx,My]=meshgrid(x,y);

%Divergencia
div=Mx.*0+My.*0;

%Pintamos la divergencia que será de un único color al ser nula
surf(Mx,My,div)
shading flat
axis([-0.5,10.5,-0.5,2.5]);
axis equal
view(2);
title('Divergencia del vector desplazamiento (nula)');
xlabel('X');
ylabel('Y');
colorbar


7 Rotacional del sólido

Calcular [math]|∇ × \vec{u}|[/math] en todos los puntos del sólido en [math]t = 0[/math] y dibujarlo. ¿Qué puntos sufren un mayor rotacional?

Calcularemos el rotacional:

[math]∇ × \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 & 2/5senx & 0\end{vmatrix} = 2/5cosx \vec{k}[/math]

El módulo será

[math]|∇ × \vec{u}|= 2/5cosx[/math]

Si observamos la gráfica, podemos decir que los puntos que sufren un mayor rotacional son los que se encuentran en las posiciones [math]x=0[/math] y [math]x=6,4[/math].

Para realizar el código de matlab, repetiremos los procesos anteriores iniciales, pero esta vez definiremos la variable [math]rot[/math] (el módulo del rotacional), y representaremos la gráfica en función de ella y de las matrices de [math]x[/math] y de [math]y[/math].


Rotacional
clc
clear

%Definición de regiones
h=2/10; 
x=[0:h:10];          
y=[0:h:2];

%Matriz de x e y
[Mx,My]=meshgrid(x,y);

%Módulo del rotacional
rot = (2/5)*cos(Mx);
%Representación gráfica
surf(Mx,My,rot)
shading flat
axis equal
colorbar
view(2);
axis([-0.5,10.5,-0.5,2.5]);
title('Módulo del rotacional');
xlabel('X');
ylabel('Y')


8 Tensor de tensiones

Definamos [math]ϵ(\vec{u}) = (∇\vec{y} + ∇\vec{u}^t)/2[/math], la parte simétrica del tensor gradiente de [math]\vec{u}[/math] conocido como tensor de deformaciones. En un medio elástico lineal, isótropo y homogéneo los desplazamientos permiten escribir el tensor de tensiones [math]σij[/math] a través de la fórmula
[math]σ = λ∇ · \vec{u} 1 + 2µϵ[/math]
,

donde [math]1[/math] es el tensor identidad en el conjunto de vectores libres del espacio [math]R^3[/math] y [math]λ[/math], [math]µ[/math] son los conocidos como coeficientes de Lamé que dependen de las propiedades elásticas de cada material. A pesar de que los desplazamientos son planos (es decir [math]\vec{u}[/math] no tiene componente en la dirección de [math]\vec{k}[/math]) las tensiones no tienen por qué ser planas y puede haber tensiones en la dirección ortogonal al plano de la placa. Tomando [math]λ = µ = 1[/math], dibujar las tensiones normales en la dirección que marca el eje [math]\vec{i}[/math], es decir [math]\vec{i}· σ ·\vec{i}[/math], las tensiones normales en la dirección que marca el eje [math]\vec{j}[/math], es decir [math]\vec{j} · σ · \vec{j}[/math] y las correspondientes al eje [math]\vec{k}[/math], es decir [math]\vec{k}· σ · \vec{k}[/math] (dibujar las que no son nulas).


Calculamos [math]ϵ(\vec{u})[/math]:

[math]Ԑ(\vec{u}) = \begin{pmatrix} 0 & \frac{1}{5}cosx & 0 \\ \frac{1}{5}cosx & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}[/math]

Conocido ya [math]∇ · \vec{u][/math], [math]σ[/math] será:

[math]σ = λ∇ · \vec{u} 1 + 2µϵ = 2·Ԑ(\vec{u})= \begin{pmatrix} 0 & \frac{2}{5}cosx & 0\\ \frac{2}{5}cosx & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}[/math]

Aplicando esto a las tensiones en las direcciones [math]i[/math], [math]j[/math] y [math]k[/math], obtenemos:

[math]\vec{1}· σ · \vec{k} = 0[/math]

[math]\vec{j}· σ · \vec{k} = 0[/math]

[math]\vec{k}· σ · \vec{k} = 0[/math]

Al ser todas nulas, no se representarán.

9 Tensiones tangenciales

Calcular las tensiones tangenciales respecto al plano ortogonal a [math]\vec{i}[/math], es decir [math]|σ ·\vec{i} − (\vec{i} · σ ·\vec{i})\vec{i}|[/math], en [math]t = 0[/math]. Dibujar sólo las que no son nulas.

[math]|σ·\vec{i} − (\vec{i}·σ·\vec{i})\vec{i}|= |σ·\vec{i}| = |\begin{pmatrix} 0 & \frac{2}{5}cosx & 0 \\ \frac{2}{5}cosx & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix}\begin{pmatrix} 1\\0\\0 \end{pmatrix}| = |\begin{pmatrix} 0\\\frac{2}{5}cosx\\0 \end{pmatrix}| = \frac{2}{5}cosx [/math]

10 Tensión de Von Mises

La tensión de Von Mises se define por la fórmula

[math]σ_{VM}=\sqrt{\frac{(σ_{1}-σ_{2})^2+(σ_{2}-σ_{3})^2+(σ_{3}-σ_{1})^2}{2}}[/math]

donde [math]σ_{1}[/math], [math]σ_{2}[/math] y [math]σ_{3}[/math] son los autovalores de [math]σ[/math] (también conocidos como tensiones principales). Se trata de una magnitud escalar que se suele usar como indicador para saber cuando un material inicia un comportamiento plástico (y no elástico puro). Pintar la tensión de Von Mises y señalar en qué punto se alcanza el mayor valor.

Como podemos ver en la gráfica, los puntos de mayor valor se alcanzan en [math]x=0[/math], [math]x=3,2[/math], [math]x=6,4[/math] y [math]x=9,6[/math].

Para el desarrollo del código, repetimos los procesos anteriores, crenado ahora laa variable [math]MVonM[/math] y definiendo la función de Von Mises ([math]VonMises[/math]). creamos un bucle en el que asignamos los valores de la tensión en cada punto y, por último, representamos la gráfica.


Tensión de Von Mises
clear
clc

%Definición de regiones
h=2/10; 
x=[0:h:10];          
y=[0:h:2];

%Matriz de X e Y
[X,Y]=meshgrid(x,y);
MVonM=0.*X;

%Definimos la funcion de Von mises donde tp1,2,3 son las tensiones principales
VonMises=inline('(((tp1-tp2)^2+(tp2-tp3)^2+(tp3-tp1)^2)/2)^(1/2)','tp1','tp2','tp3');
[f,c]=size(X);

%Le asignamos a la matriz MVonM los valores de la tension de Von Mises en cada punto
for i=1:f
 for j=1:c
  deformaciones=[[0;(2/5).*cos(X(i,j));0],[(2/5).*cos(X(i,j));0;0],[0;0;0]];   
  lamdas=eig(deformaciones);
  tp1=lamdas(1,1);
  tp2=lamdas(2,1);
  tp3=lamdas(3,1);
  MVonM(i,j)=VonMises(tp1,tp2,tp3);
 end
end

%Pintamos
surf(X,Y,MVonM)
shading flat
axis([-0.5,10.5,-0.5,2.5]);
axis equal
title('Tension de Von Mises');
xlabel('X');
ylabel('Y');
view(2);
colorbar


11 Campo de fuerzas que actúa sobre la placa

El campo de fuerzas [math]\vec{F}[/math] que actúa sobre la placa (y que son las causantes del desplazamiento observado) se aproxima usando la ecuación de la elasticidad lineal

[math]\vec{F}=\frac{∂^2\vec{u}}{∂t^2}-∇· σ[/math]

donde [math]∇ · σ[/math] es el campo vectorial que se obtiene al hacer la divergencia de los vectores cuyas componentes son las filas de la matriz [math]σ[/math]. Calcular la velocidad de propagaci´on de las ondas [math]v[/math] en términos de las constantes de Lamé, suponiendo que [math]\vec{F} = 0[/math]. Si la onda fuera longitudinal, es decir, tomando [math]\vec{a} = 2/5\vec{i}[/math], ¿cuál sería la velocidad de propagación? Comprobar que sobre un mismo medio las ondas transversales y longitudinales no viajan a la misma velocidad, tal y como se observa en la transmisión de ondas sísmicas.

12 Módulo del desplazamiento vertical

Fijado ahora el punto [math](1/2, 1)[/math], calcular el módulo del desplazamiento vertical (dirección [math]\vec{j}[/math]) a lo largo del tiempo en el intervalo [math]t ∈ [0, 10][/math]

Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Observación_de_campos_escalares_y_vectoriales_en_elasticidad&oldid=49429»