Campos escalares y vectoriales en elasticidad (Grupo 3B)

De MateWiki
Revisión del 20:42 14 dic 2023 de Mario.raya (Discusión | contribuciones) (→‎Cálculo de la velocidad de propagación)

Saltar a: navegación, buscar
Trabajo realizado por estudiantes
Título Campos escalares y vectoriales en elasticidad (Grupo 3B)
Asignatura Teoría de Campos
Curso 2023-24
Autores Eladio Rodríguez Rúa
Jorge Granadino Aranda
Mario Raya Sampere
Alejandro Villaverde Carrascosa
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura


El trabajo consiste en la visualización de campos escalares y vectoriales en elasticidad. Para ello, se ha utilizado principalmente el programa informático MATLAB que permite ver los cálculos de manera más visual.

Se considera una placa rectangular (en dimensión 2) que ocupa la región [math](x, y)∈[-1, 1] × [0, 12][/math].

En ella se supone que hay dos cantidades físicas definidas: La temperatura [math]T(x,y) [/math] que viene dada por:
[math]T(x,y) = 3log(1+(x-1)^2) + log(1+(y-8)^2)[/math]
y los desplazamientos [math] \vec {u}(x,y) [/math] producidos por la acción de una fuerza determinada. De esta forma, si se define [math] \vec{r_{0}}(x,y) = x\vec {i} + y\vec {j}[/math] como el vector de posición de los puntos de la placa antes de la deformación, la posición de cada punto (x,y) 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]

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(\pi k(\vec {d}· \vec {r_{0}}(x,y)-vt))[/math]

donde [math]\vec {a}[/math] se conoce como aplitud, [math] k\gt0 [/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 se detiene en [math] t = 0 [/math] en los primeros 10 apartados de este trabajo. De manera que, solo para los primeros apartados:

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

Se trata de una onda transversal en la que la dirección de propagación es ortogonal a la amplitud. En particular se toma:

[math]\vec {a}(x,y)=\frac{x}{3}\vec {i}[/math]  ; [math]\vec {d}=\frac{1}{12}\vec {j}[/math]  ; [math]\vec {k}= 1 [/math]



1 Dibujo del mallado

En primer lugar dibujamos el mallado de nuestra placa rectangular plana para representar los puntos interiores del sólido. Tomamos los ejes del rectángulo [-1; 1] × [0; 12] y como paso de muestreo h = 2/10 para las variables x e y. Se han utilizado los comandos axis() y meshgrid() en MATLAB.

REPRESENTACIÓN DEL MALLADO
% Mallado interior de la figura:
%definimos el paso de muestreo y las variables
h = 2/10;
x = -1:h:1;
y = 0:h:12;
%se realiza el mallado
[X,Y] = meshgrid(x,y);
mesh(X,Y,0*X);
%se define el rango de visión de la gráfica
axis([-10,10,-0.5,12.5]);
view(2)
%se pone título a la gráfica
title('Representación del mallado')
%se pone nombre a los ejes
xlabel('Eje x')
ylabel('Eje y')


2 Dibujo de las curvas de nivel de la temperatura y de [math]\nabla T[/math]

A continuación dibujamos mediante colores las curvas de nivel de la temperatura, que viene dada por el campo escalar [math]T(x,y) = 3log(1+(x-1)^2) + log(1+(y-8)^2)[/math]. Los colores azulados representan las zonas más frías, mientras que las zonas más cálidas tienen tonos más amarillentos y anaranjados. A partir de la gráfica podemos ver que el punto más próximo a la temperatura máxima es el (-1,0).

Se puede comprobar en la gráfica del gradiente que existe ortogonalidad entre el gradiente y las curvas de nivel del campo de temperatura.

CURVAS DE NIVEL Y GRADIENTE DE LA PLACA
% Curvas de nivel:
subplot(2,1,1)
%definimos el paso de muestreo y las variables
h=2/10;
x=-1:h:1;
y=0:h:12;
%se realiza el mallado
[X,Y]=meshgrid(x,y);
%se define la temperatura
T=3*log(1+(X-1).^2)+log(1+(Y-8).^2);
%dibujamos las curvas de nivel
contour(X,Y,T,40)
%barra de indicación de colores
colorbar
%se define el rango de visión de la gráfica
axis([-5,5,-0.5,12.5])
%titulamos la gráfica y los ejes
title('Curvas de nivel')
xlabel('Eje x')
ylabel('Eje y')

% Gradiente de la placa:
subplot(2,1,2)
%definimos el paso de muestreo y las variables
h=2/10;
x=-1:h:1;
y=0:h:12;
%se realiza el mallado
[X,Y]=meshgrid(x,y);
%se define la temperatura
T=3*log(1+(X-1).^2)+log(1+(Y-8).^2);
%dibujamos las curvas de nivel
contour(X,Y,T)
%se define el rango de visión de la gráfica
axis([-5,5,-0.5,12.5])
view(2)
%barra de indicación de colores
colorbar
%se calcula el gradiente en la misma gráfica
hold on
[Px,Py]=gradient(T)
quiver(X,Y,Px,Py)
%se define el rango de visión de la gráfica
axis([-5,5,-0.5,12.5])
view(2)
%titulamos la gráfica y los ejes
title('Gradiente de la placa')
xlabel('Eje x')
ylabel('Eje y')


3 Ley de Fourier

De acuerdo con la Ley de Fourier, la energía calorífica [math] \vec{Q} [/math] viaja de acuerdo con la formula [math]\vec Q [/math] [math] = −k \nabla T[/math]. Donde [math] k [/math] es la constante de conductividad térmica de la placa, que supondremos que es [math] k = 1 [/math].

Para la resolución de este apartado se reutilizará el código del apartado anterior y se hará uso de la función gradiente propia del programa MATLAB. Una vez calculado el gradiente, lo multiplicamos por la constante de conductividad térmica [math] k [/math] y representamos en una gráfica la energía calorífica [math] \vec{Q} [/math] como campo vectorial.

LEY DE FOURIER
% Cálculo de la energía calorífica:
%definimos el paso de muestreo y las variables
h=2/10;
k=1;
x=-1:h:1;
y=0:h:12;
%creación del mallado
[X,Y]=meshgrid(x,y);
%se define la temperatura
T=3*log(1+(X-1).^2)+log(1+(Y-8).^2);
hold on
%cálculo del gradiente
[Px,Py]=gradient(T);
%ley de Fourier
Qx=-k.*Px;
Qy=-k.*Py;
%campo de vectores
quiver(X,Y,Qx,Qy)
%se define el rango de visión de la gráfica
axis([-5,5,-0.5,12.5])
view(2)
%se pone título a la gráfica
title('Energía calorífica')
%se pone nombre a los ejes
xlabel('Eje x')
ylabel('Eje y')


4 Dibujo del campo de vectores

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

Comenzamos definiendo todas las variables de nuestro enunciado. A continuación,

CAMPO DE VECTORES
% Campo de vectores:
%definimos las variables
h=2/10;
k=1;
x=-1:h:1;
y=0:h:12;
%creación del mallado
[X,Y]=meshgrid(x,y);
hold on
Ux=(X./3).*sin((pi/12).*Y);
Uy=0.*Y;
%campo de vectores
quiver(X,Y,Ux,Uy)
%se define el rango de visión de la gráfica
axis([-5,5,-0.5,12.5])
view(2)
%titulamos la gráfica y los ejes
title('Campo de vectores en t=0')
xlabel('Eje x')
ylabel('Eje y')


5 Dibujo del sólido antes y después del desplazamiento

DESPLAZAMIENTO Y COMPARACIÓN
% Desplazamiento y comparación:
%antes del desplazamiento
h=2/10;
x=-1:h:1;
y=0:h:12;
subplot(2,2,1)
[X,Y] = meshgrid(x,y);
mesh(X,Y,0*X);
axis([-5,5,-0.5,12.5])
view(2)
title('Situación inicial')
xlabel('Eje x')
ylabel('Eje y')
axis equal
%después del desplazamiento
subplot(2,2,2)
Ux=(X./3).*sin((pi/12).*Y);
Uy=0.*Y;
mesh(X+Ux,Y+Uy,0.*X)
axis([-5,5,-0.5,12.5])
view(2)
title('Situación final')
xlabel('Eje x')
ylabel('Eje y')
axis equal
%comparación
subplot(2,2,3)
plot3(X,Y,0.*X,X+Ux,Y+Uy,0.*X)
axis([-5,5,-0.5,12.5])
view(2)
title('Comparación')
xlabel('Eje x')
ylabel('Eje y')
axis equal
%campo de vectores
subplot(2,2,4)
quiver(X,Y,Ux,Uy)
axis([-5,5,-0.5,12.5])
view(2)
title('Campo de Vectores')
xlabel('Eje x')
ylabel('Eje y')
axis equal


6 Cálculo y dibujo de la divergencia

A continuación calcularemos la divergencia del campo [math] \vec u [/math], que es una medida del cambio de volumen local debido al desplazamiento. La podemos calcular gracias a la siguiente expresión:
[math]\nabla \cdot \vec u(x,y) = \frac{\partial u_1}{\partial x}+\frac{\partial u_2}{\partial y}[/math]


DIVERGENCIA.
% Cálculo y dibujo de la divergencia
%definimos las variables
h=2/10;
k=1;
x=-1:h:1;
y=0:h:12;
%creación del mallado
[X,Y]=meshgrid(x,y);
Ux=(X./3).*sin((pi/12).*Y);
Uy=0.*Y;
%cálculo de la divergencia
div=divergence(X,Y,Ux,Uy);
MaxD=max(max(div));
MinD=min(min(div));
surf(X,Y,div)
%barra de indicación de colores
colorbar
%se define el rango de visión de la gráfica
axis([-2.5,2.5,-0.5,12.5,0,0.5])
%titulamos la gráfica y los ejes
title('Divergencia')
xlabel('Eje x')
ylabel('Eje y')
zlabel('Eje z')


7 Cálculo y dibujo del rotacional

En este apartado buscamos calcular el rotacional del campo [math]\vec u[/math] de los puntos del sólido en t=0. Para el calculo del rotacional utilizamos los siguientes cálculos:

[math] \nabla \times \vec u(x,y,z) = \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}[/math]

En nuestro caso tenemos:

[math] \nabla \times \vec u(x,y,z) = \begin{vmatrix} \vec{i} & \vec{j} & \vec{k}\\ \frac{\partial }{\partial x} & \frac{\partial }{\partial y} &\frac{\partial }{\partial z} \\ \frac{x}{3}sen(\frac{\pi}{12}y) & 0 & 0\end{vmatrix} =\frac{-1}{36}xπcos(\frac{π}{12}y) \vec{k}[/math]

Siendo el módulo:

[math]|\nabla \times \vec u(x,y,z)|= \frac{-1}{36}xπcos(\frac{π}{12}y)[/math]
Campo de vectores.
h=2/10;
k=1;
x=-1:h:1;
y=0:h:12;
[X,Y]=meshgrid(x,y);
Ux=(X./3).*sin((pi/12).*Y);
Uy=0.*Y;
[Rz,V]=curl(X,Y,Ux,Uy);
MaxR=max(max(Rz));


8 Dibujo de las tensiones normales

9 Cálculo de las tensiones tangenciales

En este apartado, calcularemos las tensiones tangenciales respecto al plano ortogonal a [math]\vec {i}[/math], es decir, [math]|\sigma\cdot \vec{i} - (\vec{i}\cdot \sigma \cdot \vec{i}) \cdot \vec{i}|.[/math]

10 Dibujo de la tensión de Von Mises

En este apartado vamos a trabajar con la fórmula de Von Mises, la tensión viene dada por la siguiente expresión:

[math] \sigma_{VM} =\sqrt{\frac{(\sigma_1-\sigma_2)^2+(\sigma_2-\sigma_3)^2+(\sigma_3-\sigma_1)^2}{2}} \,, [/math]

donde [math] \sigma_1 [/math], [math] \sigma_2 [/math] y [math] \sigma_3 [/math] son los autovalores de [math] \sigma [/math], también conocidos como tensiones principales. La tensión de Von Mises es una magnitud escalar que se suele usar como indicador para saber cuándo un material inicia un comportamiento plástico.

11 Cálculo de la velocidad de propagación

Para este apartado, vamos a usar la ecuación de la elasticidad lineal, definida por
[math]\vec{F}=\frac{\partial^2\vec{u}}{\partial t^2}-∇·σ[/math]

Para comenzar, definimos nuestra función [math]\vec{u}(x,y)=\frac{x}{3}sen(\frac{π}{12}y-πvt)·\vec{i}[/math].

Para continuar, definimos la función sigma que viene dada de la siguiente forma [math]σ=λ∇·\vec{u}1+2μє [/math].

Por último debemos definir algunas de las funciones que forman la funcion de tensiones sigma
[math]∇·\vec{u}=\frac{1}{3}sen(\frac{π}{12}y-πvt)[/math]


[math]є(\vec{u})=\begin{pmatrix} \frac{1}{3}sen(\frac{π}{12}y-πvt) & \frac{1}{72}πxcos(\frac{π}{12}y-πvt) & 0\\ \frac{1}{72}πxcos(\frac{π}{12}y-πvt) & 0 & 0\\ 0 & 0 & 0 \end{pmatrix} [/math]</center>

Una vez tenemos estas dos funciones podemos hallar ya la función sigma, que nos da la siguiente matriz
[math] σ= \begin{pmatrix} [λ+2μ]\frac{1}{3}sen(\frac{π}{12}y-πvt) & \frac{πx}{36}μcos(\frac{π}{12}y-πvt) & 0 \\ \frac{πx}{36}μcos(\frac{π}{12}y-πvt) & \frac{λ}{3}sen(\frac{π}{12}y-πvt) & 0 \\ 0 & 0 & \frac{λ}{3}sen(\frac{π}{12}y-πvt) \end{pmatrix} [/math]</center>

12 Dibujo de la función

En el último aparatado calcularemos el módulo del desplazamiento transversal (dada por la dirección [math] \vec i [/math]) a lo largo del tiempo en el intervalo t € [0,10].

FUNCIÓN U EN FUNCIÓN DEL TIEMPO
h=2/10;
mu=1;
lambda=-mu;
v=sqrt(mu/144);
t=0:h:10;
x=-1:h:1;
y=0:h:12;
[X,Y]=meshgrid(x,y);
x0=1/2;
y0=1;
u=x0/3*sin((pi*y0/12)-pi.*t.*v);
%creamos la gráfica
figure
%tazamos la función en función del tiempo
plot(t,u);
%Titulamos la gráfica y sus ejes
title('Función u en función del tiempo');
xlabel('Tiempo (s)');
ylabel('Valor de la función');
% Mostrar la cuadrícula
grid on;
Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Campos_escalares_y_vectoriales_en_elasticidad_(Grupo_3B)&oldid=64307»