Visualización de campos escalares y vectoriales (grupo A2)
| Trabajo realizado por estudiantes | |
|---|---|
| Título | Visualización de campos escalares y vectoriales (grupo A2) |
| Asignatura | Teoría de Campos |
| Curso | 2020-2021 |
| Autores | CHENGYI XU JINGJING MA YIJING ZHENG FANGZHENG YANG |
| Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura | |
Visualización de campos escalares y vectoriales en elasticidad. Consideramos una placa rectangular plana (en dimensión 2) que ocupa la región (x, y) ∈ [0, 10] × [−1, 1]. En ella vamos a suponer que tenemos definidas dos cantidades físicas: la temperatura [math]T(p,\theta)[/math], que viene dada por,
- [math]\vec T(p,\theta ) =p^2cos(\theta ) [/math]
y los desplazamientos ~u(x, y) producidos por la acción de una fuerza determinada. De esta forma, si definimos ~r0(x, y) 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}=\vec{ro}(x,y)+\vec{u}(x,y) [/math]
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
- [math]\vec{u}(x,y)=\frac{xy^2}{20}+f(y)\vec{j}.[/math]
donde \(f(y)\) es una cierta función.
Contenido
1 Ejercicio 1
Dibujar un mallado que represente los puntos interiores del sólido. Tomar los ejes (comando axis) en el rectángulo (x, y) ∈ [−0.5; 10.5] × [−1.5; 1.5] y como paso de muestreo h = 1/10 para las variables x e y.
En primer lugar representamos los puntos interiores de la placa con un mallado usando la función ‘mesh’, el mallado está definida entre el rectángulo (x, y) ∈ [−0.5; 10.5] × [−1.5; 1.5], con un paso de muestreo h=0.1 para los ejes de coordenadas, que usamos la función ‘axis’.
x=0:0.1:10;
y=-1:0.1:1;
% La region la placa y su muestreo
[xx,yy]=meshgrid(x,y);
%Genelacion del mallado
figure(1)
mesh(xx,yy,0*xx)
axis([-0.5,10.5,-1.5,1.5])
%Valor de los ejes de la grafica
title('Mallado de la placa plana rectangular');
% Título del gráfico
xlabel('EJE X');
% Título del eje x
ylabel('EJE Y');
% Título del eje y
view(2)
2 Ejercicio2
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.
Dibujamos las curvas de nivel de la temperatura que se define como [math]T\left ( \rho ,\theta \right )= \rho ^{2}cos(\theta ) [/math].[math] T\left ( \rho ,\theta \right )= \rho ^{2}cos(\theta )[/math] Como [math]x=\rho cos\theta e y=\rho sen\theta[/math] se deduce que [math]T=\rho x=x\sqrt{x^{2}+y^{2}}[/math] . Representando las curvas de nivel con la función ’contour’ y la función ‘colorbar’ para mostrar la escala de colores que nos define la temperatura, a mayor temperatura tiene un color más claro y a menor temperatura tiene un color más oscuro. En la gráfica se muestra que la temperatura se aumenta al desplazar a la derecha, con una temperatura máxima de 100.4988.
x=0:0.1:10;
y=-1:0.1:1;
[xx,yy]=meshgrid(x,y);
%Creacion del mallado
figure(1)
f=xx.*(sqrt(xx.^2+yy.^2));
%Funcion de la temperature
contour(xx,yy,f)
hold on
plot(x,-1+x-x,'k','linewidth',1);
plot(x,1+x-x,'k','linewidth',1);
plot(y-y,y,'k','linewidth',1);
plot(10+y-y,y,'k','linewidth',1);
title('FUNCIÓN DE LA TEMPERATURA');
%Título de la grafica
xlabel('EJE X');
%Título del eje x
ylabel('EJE Y');
%titulo del eje y
axis([-0.5,10.5,-1.5,1.5])
%Valor de los ejes de la grafica
view(2)
colorbar
Tmax=max(max(f))
hold off
3 Ejercicio3
Calcular ∇T y dibujarlo como campo vectorial. Observar gráficamente que ∇T es ortogonal a dichas curvas. (Nota: es muy importante que la escala de los ejes sea la misma para apreciar el ángulo que forman las curvas de nivel y el gradiente)
Primero calculamos el gradiente de la temperatura: [math]\triangledown [/math]T donde la función de la temperatura es[math] T\left ( \rho ,\theta \right )= \rho ^{2}cos(\theta )[/math] Como [math]x=\rho cos\theta e y=\rho sen\theta[/math] se deduce que [math]T=\rho x=x\sqrt{x^{2}+y^{2}}[/math] . Y el gradiente de un campo escalar T podemos definir como:[math]\triangledown T=\frac{\partial T}{\partial x}\overrightarrow{i}+\frac{\partial T}{\partial y}\overrightarrow{j}+\frac{\partial T}{\partial z}\overrightarrow{k} [/math]. Calculando las derivadas parciales[math]\frac{\partial T}{\partial x}=\frac{2x^2+y^2}{\sqrt{x^2+y^2}} y \frac{\partial T}{\partial y}=\frac{xy}{\sqrt{y^2+x^2}} [/math]. Con el grafico podemos ver que cumplen la propiedad del gradiente de que los:[math] \triangledown T [/math]son ortogonales a las superficies de nivel. Sobre la representación gráfica, la creación de las curvas de nivel es lo mismo que el apartado anterior, y añadimos los gradientes para hacer una comparación. Para ello usamos la función ‘quiver’ para representar el campo vectorial calculado.
x=0:0.1:10;
y=-1:0.1:1;
[xx,yy]=meshgrid(x,y);
figure(1)
f=xx.*(sqrt(xx.^2+yy.^2));
contour(xx,yy,f)
hold on
plot(x,-1+x-x,'k','linewidth',1);
plot(x,1+x-x,'k','linewidth',1);
plot(y-y,y,'k','linewidth',1);
plot(10+y-y,y,'k','linewidth',1);
title('GRADIENTE DE TEMPERATURA');
xlabel('EJE X');
ylabel('EJE Y');
axis([-0.5,10.5,-1.5,1.5])
view(2)
colorbar
x=0:0.1:10;
y=-1:0.1:1;
[xx,yy]=meshgrid(x,y);
u=((2.*xx.^2)+(yy.^2))./sqrt((xx.^2)+(yy.^2));
v=(xx.*yy)./sqrt((yy.^2)+(xx.^2));
%Gradientes de la función temperatura
quiver(xx,yy,u,v)
%Representacion del campo vectorial.
view(2)
hold off
4 Ejercicio4
Consideramos ahora el campo de desplazamientos [math] \vec{u} [/math]. Calcularlo sabiendo: El campo de desplazamiento dado es : [math] \overrightarrow{u}=\frac{xy^2}{20}\overrightarrow{i}+f(y)\overrightarrow{j} [/math],
Y el enunciado nos dicen:
1. Los puntos situados en el eje y = 0 no sufren desplazamiento en la dirección [math] \overrightarrow{i} [/math]. 2.[math] \triangledown \cdot u=\frac{y^2-2}{20} [/math]. Como la divergencia del campo vectorial es: [math]\triangledown \cdot \overrightarrow{u}=\frac{\partial u_{1}}{\partial x}+\frac{\partial u_{2}}{\partial y}+\frac{\partial u_{3}}{\partial z}[/math], aplicando esto igualando al [math] \triangledown\cdot u[/math] obteniendo el resultado de [math]\frac{y^2}{20}+{f}'(y)=\frac{y^2-2}{20} [/math], sacando que[math] f(y)=\frac{-y}{10} [/math]. Igualando las ecuaciones tenemos que: [math]\overrightarrow{u}=\frac{xy^2}{20}-\frac{y}{10}[/math]
5 Ejercicio5
Dibujar el campo de vectores en los puntos del mallado del sólido.
Disponiendo de la fórmula del campo de desplazamientos del apartado anterior, podemos dibujar el campo de vectores en los puntos del mallado del solido, primer con meshgrid se transforma a dos matrices asociado a cada coordenada y los dos vectores al función del resultado de meshgrid, representamos la grafica con la función de la biblioteca mesh y quiver para presentar el campo de vector
x=0:0.1:10;
y=-1:0.1:1;
[xx,yy]=meshgrid(x,y);
u=(xx.*yy.^2)/20;
v=-(yy/10)-((xx.*yy.^2)/20);
figure
mesh(xx,yy,0.*xx)
hold on
quiver(xx,yy,u,v,'k')
axis([-0.5,10.5,-1.5,1.5])
view(0,90)
hold off
6 Ejercicio6
Dibujar el sólido antes y después del desplazamiento dado por el campo de vectores [math] \vec{u} [/math]. Dibujar ambos en la misma figura usando el comando subplot.
Definimos como siempre el rectángulo de la malla utilizando la función meshgrid para convertir a dos matrices respectivamente asociada, la placa sin deformación lo representamos mediante la función con el valor cambiado al solución de la función meshgrid, utilizando subplot para dibujar en la misma figura y añadimos la placa influido con respecto la deformación al lado.
x=0:0.1:10;
y=-1:0.1:1;
[xx,yy]=meshgrid(x,y);
u=(xx.*yy.^2)./20;
v=-yy./10;
Rx=u+xx;
Ry=v+yy;
subplot(1,2,1);
mesh(xx,yy,0*xx);
axis([-0.5,10.5,-1.5,1.5])
title('PLACA SIN DEFORMACIÓN');
xlabel('Eje x');
ylabel('Eje y');
zlabel('Eje z');
subplot(1,2,2);
mesh(Rx,Ry,0*Rx)
axis([-0.5,10.5,-1.5,1.5])
title('PLACA CON DEFORMACIÓN');
xlabel('Eje x');
ylabel('Eje y');
zlabel('Eje z');
7 Ejercicio7
Dibujar ∇ · [math] \vec{u} [/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?
Primero obtenemos el resultado del calculo de la divergencia que es la suma de las derivadas en cada vector. Representamos gráficamente, nombrando a cada eje y con la función de la biblioteca pcolor,countour,surf mostramos la grafica y colorbar refrescamos su altura,con hold on y hold off representamos en la misma grafica y utilizamos el max y min para definir los valore máximos y mínimos.
x=0:0.1:10;
y=-1:0.1:1;
[xx,yy]=meshgrid(x,y);
divu=(yy.^2-2)/20;
surf(xx,yy,divu)
shading interp
colorbar
title('divergencia')
xlabel('x')
ylabel('y')
zlabel('z')
figure
pcolor(xx,yy,divu)
shading interp
hold on
contour(xx,yy,divu,'k')
xlabel('x')
ylabel('y')
hold off
figure
plot(yy,divu)
title('proyección sobre el eje YOZ')
xlabel('y')
ylabel('z')
grid on
maxi=max(max(divu))
mini=min(min(divu))maxi = -0.050000 mini = -0.10000
8 Ejercicio8
Calcular | ∇ × [math] \vec{u} [/math]| en todos los puntos del sólido y dibujarlo. ¿Qué puntos sufren un mayor rotacional?
Para calcular esto aplicamos la formula de rotacional, y a la hora de representar gráficamente usamos como siempre en primer lugar meshgrid, para devolver dos matrices con la coordenada X e Y de los puntos de la retícula obtenida y el rotacional en función con respecto la solución de meshgrid, representamos el rotacional a partir del mallado rectangular al final con el criterio max para calcular el valor máximo.
[math] \left | \bigtriangledown \times u \right |[/math] =[math] \begin{vmatrix}\vec{i} & \vec{j}& \vec{k}\\ \frac{\partial }{\partial x}& \frac{\partial }{\partial y}& \frac{\partial }{\partial z}\\ F_{1}& F_{2}& F_{3} \end{vmatrix}[/math] = [math] \begin{vmatrix} \vec{i} & \vec{j} & \vec{z}\\ \frac{\partial }{\partial x} &\frac{\partial }{\partial y} &\frac{\partial }{\partial z} \\ \frac{xy^{2}}{20} & \frac{-y}{10} & 0\end{vmatrix}[/math] =-[math]\frac{xy}{10}\vec{k} [/math].
x=0:0.1:10;
y=-1:0.1:1;
%Creación del mallado
[xx,yy]=meshgrid(x,y);
rot=-(xx.*yy)./10;
%representacion del rotacional a partir del mallado
surf(xx,yy,rot)
colorbar
axis([-0.5,10.5,-1.5,1.5])
view(2)
title('Rotacional')
xlabel('EJE X')
ylabel('EJE Y')
zlabel('EJE Z')
%valores maximos y minimos
rotmax=max(max(rot))rotacional máximo=1
9 Ejercicio9
Definamos [math] \vec{u} [/math] = (∇ [math] \vec{u} [/math] + ∇ [math] \vec{u} [/math]t)/2, 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 σij a través de la formula σ = λ∇ · [math] \vec{u} [/math] 1 + 2µ[math]\varepsilon [/math], donde λ y µ 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 λ = µ = 1, dibujar las tensiones normales en la dirección que marca el eje [math] \vec{i} [/math], es decir [math] \vec{i} [/math] · σ ·[math] \vec{i} [/math], las tensiones normales en la dirección que marca el eje [math] \vec{j} [/math], es decir [math] \vec{j} [/math] · σ · [math] \vec{j} [/math] y las correspondientes al eje [math] \vec{k} [/math], es decir [math] \vec{k} [/math] · σ ·[math] \vec{k} [/math].
- [math]\nabla \vec u = \frac{\partial u}{\partial x} \vec i + \frac{\partial u}{\partial y} \vec j = ( \frac{y^2}{20}\vec i )\cdot\vec i + (\frac{xy}{10}\vec i-\frac{1}{10}\vec j)\cdot\vec j = \frac{y^{2}}{20}\vec i\otimes \vec i +\frac{xy}{10}\vec i\otimes \vec j - \frac{1}{10}\vec j \otimes \vec j[/math]
- [math]\left ( \nabla \vec u ^{t} \right )= \frac{y^{2}}{20}\vec i\otimes \vec i +\frac{xy}{10}\vec j\otimes \vec i - \frac{1}{10}\vec j \otimes \vec j[/math]
- [math]\epsilon \left ( \vec{u}\right )=\frac{\nabla \vec u+ \nabla \vec u ^{t}}{2}= \frac{y^{2}}{20}\vec i\otimes \vec i - \frac{1}{10}\vec j \otimes \vec j+\frac{xy}{20}\left [ \left ( \vec{i}\otimes \vec{j} \right )+\left ( \vec{j}\otimes \vec{i} \right ) \right ][/math]
- [math]\sigma =\lambda \triangledown \cdot \vec{u}I+2\mu \epsilon [/math]
- [math]\lambda =\mu =1[/math]
- [math]\rightarrow \sigma =\triangledown \cdot \vec{u}\cdot I+2\cdot 1\cdot \epsilon [/math]
- [math]\triangledown \cdot \vec{u}\cdot I= \frac{y^{2}-2}{20}\cdot\begin{pmatrix} 1& 0 & 0\\ 0& 1& 0\\ 0& 0 & 1 \end{pmatrix}=\begin{pmatrix} \frac{y^{2}-2}{20}& 0 & 0\\ 0& \frac{y^{2}-2}{20} & 0\\ 0& 0& \frac{y^{2}-2}{20} \end{pmatrix}[/math]
- [math]\sigma =\begin{pmatrix} \frac{y^{2}-2}{20}& 0 & 0\\ 0& \frac{y^{2}-2}{20} & 0\\ 0& 0& \frac{y^{2}-2}{20} \end{pmatrix}+2\cdot \begin{pmatrix} \frac{y^{2}}{20} & \frac{xy}{20} & 0\\ \frac{xy}{20}& -\frac{1}{10} & 0\\ 0& 0 & 0 \end{pmatrix}=\begin{pmatrix} \frac{3y^{2}}{20} & \frac{xy}{10} & 0\\ \frac{xy}{10}& \frac{y^{2}-6}{20} & 0\\ 0 & 0 & \frac{y^{2}-2}{20} \end{pmatrix} [/math]
- [math]\vec{i}\cdot \sigma \cdot \vec{i}=\frac{3y^{2}-2}{20} [/math]
- [math]\vec{j}\cdot \sigma \cdot \vec{j}=\frac{y^{2}-6}{20}[/math]
- [math]\vec{k}\cdot \sigma \cdot \vec{k}=\frac{y^{2}-2}{20}[/math]
10 Ejercicio10
Calcular las tensiones tangenciales respecto al plano ortogonal a [math] \vec{i} [/math], es decir |σ · [math] \vec{i} [/math] − ( [math] \vec{i} [/math] · σ · [math] \vec{i} [/math]) [math] \vec{i} [/math]|
calculamos directamente con la solución obtenida de sigma y i*σ*i, y nos dice con respecto vector i(1,0,0),a través de hacer las operaciones producto de la matriz se puede obtener el resultado.
- [math]\left | \sigma \cdot \vec{i}-\left ( \vec{i}\cdot \sigma \cdot \vec{i} \right )\vec{i} \right |=\left | \begin{pmatrix} \frac{3y^{2}}{20} & \frac{xy}{10} & 0\\ \frac{xy}{10}& \frac{y^{2}-6}{20} & 0\\ 0 & 0 & \frac{y^{2}-2}{20} \end{pmatrix} \begin{pmatrix} 1\\ 0\\ 0 \end{pmatrix}-(\frac{3y^{2}-2}{20})\begin{pmatrix} 1\\ 0\\ 0 \end{pmatrix} \right |=\frac{xy}{10}[/math]
11 Ejercicio11
La tensión de Von Mises se define por la fórmula σvm=[math]\sqrt{\frac{(σ1-σ2)^2+(σ2-σ3)^2+(σ3-σ1)^2}{2}}[/math] donde σ1, σ2 y σ3 son los autovalores de σ (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. (Para calcular autovalores con OCTAVE/Matlab usar el comando eig.m)
primero determina los intervalos de vector x,y, obtención del mallado a través de meshgrid, y obtenemos los valores de la tensión a partir de los intervalores.Dibujamos y podemos saber que el valor máximo=1.7521
x=0:0.1:10;
y=-1:0.1:1;
%Creación del mallado
[xx,yy]=meshgrid(x,y);
Sigma=zeros(3,3);
%Obtenemos los valores de la tensión de Von Mises a partir de los autovalores
for i=1:21
for j=1:101
r=xx(i,j);
t=yy(i,j);
Sigma(1,1)=(3.*t.^2-2)./20;
Sigma(1,2)=(r.*t)/10;
Sigma(2,1)=(r.*t)/10;
Sigma(2,2)=(t.^2-6)./20;
Sigma(3,3)=(t.^2-2)./20;
%Obtención de los autovalores de la matriz sigma
[v,d]=eig(Sigma);
%Fórmula de Von Mises
M(i,j)=sqrt(((d(1,1)-d(2,2))^2+(d(2,2)-d(3,3))^2+(d(3,3)-d(1,1))^2)/2);
end
end
%Dibujamos la tensión de Von Mises
surf(xx,yy,M)
colorbar
view(2)
title('Tensión de Von Mises')
%Valor Máximo
ValorMaximo=max(max(M))
%Coordenadas del valor Máximo
[i,j]=find(M==ValorMaximo)
xx(i,j)
yy(i,j)valor maximo=1.7521
12 Ejercicio12
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
Dibujar [math] \vec{f} [/math] e interpretar la grafica.
- [math]\vec{F}=-\triangledown \cdot \sigma =-\frac{\partial \sigma ji}{\partial xj}\vec{ei}[/math]
- [math]\vec{F1}=-\frac{\partial \sigma j1}{\partial x}\vec{i}=-\left [ \frac{\partial \frac{3y^{2}}{20}}{\partial x}+\frac{\partial \frac{xy}{10}}{\partial y}+\frac{\partial 0}{\partial z}\right ]\vec{i}=-\frac{x}{10}\vec{i}.[/math]
- [math]\vec{F2}=-\frac{\partial \sigma j2}{\partial y}\vec{j}=-\left [ \frac{\partial \frac{xy}{10}}{\partial x} +\frac{\partial \frac{y^{2}-6}{20}}{\partial y}+\frac{\partial 0}{\partial z}\right ]\vec{j}=-\frac{y}{5}\vec{j}.[/math]
- [math]\vec{F3}=-\frac{\partial \sigma j3}{\partial z}\vec{k}=-\left [ \frac{\partial 0}{\partial x}+\frac{\partial 0}{\partial y}+\frac{\partial \frac{y^{2}-2}{20}}{\partial z} \right ]=0\vec{k}[/math]
13 Ejercicio14
Supongamos que la densidad de la placa es d(x, y) = 1 + xy log(1 + x + y2). Calcular la masa total aproximando la integral correspondiente numéricamente.
h=0.1;
x=0:0.1:10;
y=-1:0.1:1;
%Creación del mallado
[U,V]=meshgrid(x,y);
d=(exp(U.*V))./(U.^2+V.^2+2);
ms=abs(d.*h^2);
Masa=sum(sum(ms))Masa=52.2097
