Visualización de campos escalares y vectoriales en elasticidad: Grupo 3-A

Trabajo realizado por estudiantes
Título	Deformaciones de una placa plana rectangular. Grupo 3-A
Asignatura	Teoría de Campos
Curso	2016-17
Autores	Bryan Barrantes Pérez, María de Tomás Pardo de Andrade, Paula Vázquez Bustos
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 \((x,y) ∈ [−2,2]×[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 [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 \((x,y)\) de la placa después de la deformación viene dada por: [math]\vec r(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 de los puntos de la misma dado por el vector de desplazamientos: [math]\vec u(x,y)=-\frac{x^3}{20}\vec i+f(y)\vec j[/math] donde \(f(y)\) es una cierta función.

1 Mallado

En primer lugar, es necesario dibujar un mallado que represente los puntos interiores del sólido, dando como resultado una visualización de la placa rectangular mediante sus puntos. Las dimensiones de esta malla; es decir, el número de puntos del sólido resultante de la realización del mallado, está claramente condicionado por el paso de muestreo o incremento seleccionado seleccionado; y por supuesto, condicionado por la región o dominio de la placa. Esta condición está marcada por el siguiente hecho: De cada componente del intervalo \(x\) e \(y\), parten rectas verticales y horizontales, respectivamente. Claramente, estas rectas se van a intersecar en diferentes puntos, dando lugar a la malla de la que hablamos. Por lo tanto, con este mallado obtenemos dos matrices que contienen las coordenadas \(x\) e \(y\) de todos los puntos de la malla y por tanto, de la placa. Dicho sea de paso, estas matrices presentarán la misma dimensión; es decir, compartirán el mismo número de filas y de columnas.

En este caso, hemos considerado los ejes en el rectángulo \((x,y) ∈ [−3,3] × [0,5]\) y como paso de muestreo \(h=1/10\) para las variables \(x\) e \(y\). De esta manera, para la obtención de la representación gráfica del mallado se ha empleado el código de Matlab que se puede ver a continuación:

right

%Definición del paso de muestreo como variable
pm=1/10; 
%Definición de parámetros de la superficie
x=[-2:pm:2];          
y=[0:pm:4];
%Creaci?n del mallado
[X,Y]=meshgrid(x,y);
%Dibujo de la malla, al ser mesh() un comando con tres elementos de entrada en vez de crear una matriz Z nula, se toma una de las creadas y se multiplica por el elemento nulo 
mesh(X,Y,0*X)
axis([-3,3,0,5]);
view(2);

2 Visualización del campo escalar: Temperatura

Tal como se ha mencionado al principio de este artículo; una de las cantidades físicas definidas en la placa rectangular es la temperatura del sólido \(T(x,y)\), que depende de las dos variables espaciales \((x,y)\). Ésta viene dada por la función \(T(x,y)=(x+2)^2+(y+2)^2\), representando un campo escalar que a cada punto \((x,y) ∈ [−2,2]×[0,4]\) le asigna un valor escalar; el cual, como interpretación física, se corresponde con la temperatura del sólido en dicho punto \((x,y)\).

Ahora bien; si consideramos el valor escalar, resultante de analizar nuestro campo de temperaturas para cada punto \((x,y)\) de nuestro sólido, como una altura o coordenada \(z\), obtendremos un punto \((x,y,z)\). De esta manera; el lugar geométrico, definido por este conjunto de puntos \((x,y,z)\), será una superficie representada en el espacio tridimensional.

Podemos visualizar lo dicho anteriormente en la imagen de los gráficos proporcionados por Matlab. La superficie de la que hablamos la podemos observar en el primer modelo dotada de una graduación del color condicionada por el valor de la temperatura en cada uno de sus puntos. De esta manera, observamos que el comando surfc y shading flat (responsable de la graduación e interpolación de los colores de relleno de la retícula, proporcionada por surfc), designan como predeterminado un rango de colores que va desde el azul oscuro (temperatura mínima) hasta el amarillo (temperatura máxima). Por lo tanto, los colores que van a apareciendo a medida que aumenta el valor de la función de temperaturas, se corresponden con temperaturas intermedias. Del mismo modo, la superficie de temperaturas está acompañada de sus correspondientes curvas de nivel proyectadas en el plano horizonal \(XY\) o en el plano \(Z=0\). Dicho sea de paso, esta proyección de las curvas de nivel se consigue con el comando surfc, anteriormente mencionado, tratándose de una extensión del comando surf; que añade únicamente las curvas de nivel de la superficie, dotadas del color de temperatura que le corresponde en función del valor de temperatura que representa.

Estas curvas de nivel se definen como el lugar geométrico de los puntos \((x,y)\) en los que el campo de temperaturas toma el mismo valor. Por lo tanto, para conocer cuál es el lugar geométrico que representa las curvas de nivel de nivel de nuestro campo de temperaturas, basta con igualar la función a un valor escalar determinado o bien, genérico. Geométricamente, lo que estamos haciendo es intersecar la superficie del campo de temperaturas con planos paralelos al plano \(XY\) o planos \(Z=k\); de tal manera que la intersección de dos superficies nos da como resultado las curvas de nivel del campo. Analizando esta igualdad, sabremos la forma geométrica de las curvas de nivel del campo escalar. El procedimiento sería el siguiente:

[math]T(x,y) = (x+2)^2+(y+2)^2=k [/math]

De esta manera, las curvas de nivel quedan definidas por la siguiente expresión: [math](x+2)^2+(y+2)^2=k [/math] Por lo tanto, analizando la estética o forma de la igualdad conseguimos saber, sin ninguna dificultad, qué tipo de curva representa las curvas de nivel en cuestión: circunferencias de centro \((-2,-2)\) y radio [math] R = + \sqrt x [/math].

Por otro lado, en el segundo modelo de la imagen, podemos visualizar la superficie del primer modelo proyectada en el plano horizontal, acompañada tambien de sus curvas de nivel a distintos valores. Mediante esta proyección conseguimos visualizar la placa rectangular dotada, además, de curvas de nivel y de una graduación e interpolación del color condicionada por el valor de la temperatura en cada uno de sus puntos, tal como se ha mencionado anteriormente.

Por lo tanto, en ambos modelos podemos observar perfectamente que la temperatura máxima de la placa se corresponde con el punto \((2,4)\); es decir, el extremo superior derecho. Esto se deduce fácilmente a partir de la visualización del campo y de sus curvas de nivel.

centro

Los gráficos se han obtenido empleando el código de Matlab a continuación:

pm=1/10; 
x=[-2:pm:2];          
y=[0:pm:4];
[X,Y]=meshgrid(x,y);
%Definición de la función temperatura especificando como incógnitas x e y
T=inline('((x+2).^2)+((y+2).^2)','x','y');
%Obtención de los valores de la temperatura para cada uno de los puntos del mallado
Z=T(X,Y);
%Representación gráfica de los resultados obtenidos
subplot(1,2,1);
surfc(X,Y,Z);
shading flat
xlabel('eje X')
ylabel('eje Y')
zlabel('eje Z')
subplot(1,2,2);
pcolor(X,Y,Z);
axis equal
shading flat
hold on
contour(X,Y,Z,'k')
colorbar
title('La temperatura máxima se alcanza en el punto (2,4)')
hold off

3 Visualización del campo vectorial: Gradiente de Temperatura

En esta sección, comprobaremos la ortogonalidad que existe entre el gradiente \(∇T\) del campo escalar \(T(x,y)\) y las curvas de nivel del mismo. Teóricamente, el gradiente se define como la variación de de una magnitud en función de los valores que toma, dotada de una dirección y un sentido. Por lo tanto, el gradiente de un campo tiene estructura de campo vectorial.

Para comprobar teóricamente la ortogonalidad entre el campo gradiente y las curvas de nivel es necesario introducir el concepto de derivada direccional, cuya expresión matemática es la siguiente: [math]D_{\vec{v}}{f} = \nabla f \cdot \vec{v}[/math]

Donde el vector [math]\vec{v}[/math] es un vector unitario. Esta expresión se define como la pendiente del campo en la dirección marcada por el vector [math]\vec{v}[/math]. Por lo tanto, la dirección para la cual la pendiente es nula, es la marcada por las curvas de nivel; ya que en dichas direcciones, por definición, el valor del campo no presenta ninguna variación; y por lo tanto, la pendiente es nula en dichas direcciones. Como la expresión anterior está definida por un producto escalar de vectores, podemos afirmar que obligatoriamente el campo gradiente y las curvas de nivel son ortogonales.

Teniendo en cuenta que el gradiente del campo escalar de temperaturas viene dada por la siguiente expresión:

[math]\nabla{T}=\frac{\partial{T}}{\partial {x}_{i}}{\vec {e}}_{i}=\frac{\partial{T}}{\partial{x}}{\vec{e}}_{1}+\frac{\partial{T}}{\partial{y}}{\vec{e}}_{2}={2·(x+2)}{\vec{e}}_{1}+{2·(y+2)}{\vec{e}}_{2}[/math]

elaboramos un código Matlab, mediante el cual representamos el gradiente del campo como un campo vectorial y las curvas de nivel del campo escalar; tal como podemos observar en la imagen proporcionada. De esta manera, comprobamos que efectivamente los vectores del campo gradiente son ortogonales a las curvas de nivel del campo de temperaturas.

• Dicho sea de paso, los vectores del campo gradiente no están aplicados en todos los puntos de la malla obtenida inicialmente, con la finalidad de conseguir menos vectores representados y por tanto, una mejor visualización.

El código Matlab utilizado para la visualización de la ortogonalidad del campo gradiente y las curvas de nivel del campo de temperaturas, es el siguiente:

%Definición la función temperatura tomando como incógnitas x e y
T=inline('((x+2).^2)+((y+2).^2)','x','y');
%Definición las derivadas parciales de la temperatura tomando como incógnitas en ambos casos x e y

right

TX=inline('2*(x+2)','x','y');
TY=inline('2*(y+2)','x','y');
pm=1/10;
x=[-2:pm:2];          
y=[0:pm:4];
[X,Y]=meshgrid(x,y);
Z=T(X,Y);
%Representación gráfica del gradiente de la temperatura
contour(X,Y,Z);
colorbar
hold on
%Definición de dos vectores que toman valores cada 0.25 en los intervalos definidos para x e y
xx=[-2:0.25:2];
yy=[0:0.25:4];
%Creación de un segundo mallado a partir de esos vectores
[XX,YY]=meshgrid(xx,yy);
%Obtención de valores de las derivadas parciales de la temperatura para cada punto del mallado que se acaba de definir
U=TX(XX,YY);
V=TY(XX,YY);
quiver(XX,YY,U,V);
axis equal
hold off

4 Campo de desplazamiento

Consideramos ahora el campo de desplazamientos [math]\vec u[/math]. Inicialmente, sabemos que este campo presenta la siguiente expresión: [math]\vec u(x,y)=-\frac{x^3}{20}\vec i+f(y)\vec j[/math] donde \(f(y)\) es una cierta función.

La expresión del campo de desplazamientos muestra que los diferentes puntos de la placa solo tendrán desplazamiento en la dirección de [math]\vec i[/math] y [math]\vec j[/math], tratándose de un campo vectorial definido en [math]R^2[/math].

Por otro lado, la función \(f(y)\) es desconocida en principio. Sin embargo, disponemos de los siguientes datos y características de nuestro campo de desplazamientos:

• Los puntos situados en el eje [math]y=0[/math] no sufren desplazamiento en la dirección [math]\vec j[/math]
• [math]\nabla\cdot\vec u= -\frac{3x^2}{20} + \frac{1}{10}[/math]

Al disponer de la divergencia del campo [math]\vec u(x,y)[/math] como dato, basta con calcular la divergencia del campo de desplazamientos. De esta manera, al igualar ambos resultados obtendremos la expresión de la función \(f(y)\). Como nuestro campo está definido en el sistema de coordenadas cartesianas con base ortonormal canónica [math]\{\vec i ;\vec j ;\vec k \}\,[/math], la expresión genérica de nuestro campo es la siguiente:

[math]\vec u(x,y) = u_x(x,y)\vec i + u_y(x,y)\vec j[/math]

Teniendo en cuenta la expresión anterior, la divergencia del campo [math]\vec u(x,y)[/math] viene dada por:

[math]\nabla\cdot\vec u = \frac{\partial u_x}{\partial x}+ \frac{\partial u_y}{\partial y}[/math]

De esta manera, realizando los diferentes cálculos: [math]\nabla\cdot\vec u = -\frac{2x^2}{20} + f'(y)[/math]

[math]-\frac{3x^2}{20} + \frac{1}{10} = -\frac{2x^2}{20} + f'(y)[/math]

Por lo tanto, atendiendo a la igualdad: [math]f'(y)= \frac{1}{10}[/math]

Obtenemos la función \(f(y)\) integrando la expresión anterior: [math]f'(y)= \frac{dy}{dx}= \frac{1}{10}[/math]

[math]dy= \frac{1}{10}\, dx[/math]

[math]\int\frac{1}{10}\,\text{d}x= \frac{y}{10} + C [/math] donde \(C\) es una constante de integración, tal que [math]C∈R[/math]

De momento, la expresión de nuestro campo de desplazamientos [math]\vec u(x,y)[/math] es la siguiente: [math]\vec u(x,y)=-\frac{x^3}{20}\vec i + (\frac{y}{10}+C)\vec j[/math]

Para conocer esta constante de integración, basta con utilizar el primer dato sobre los desplazamientos: Los puntos situados en el eje [math]y=0[/math] no sufren desplazamiento en la dirección [math]\vec j[/math]. Por lo tanto, analizando nuestro campo en un punto genérico que represente la condición anterior; es decir \((x,0)\), obtenemos el valor de \(C\): [math]\vec u(x,0)=-\frac{x^3}{20}\vec i + C\vec j[/math]

como sabemos que en dichos puntos no existe desplazamiento en la dirección de [math]\vec j[/math], obligatoriamente \(C=0\).

Por lo tanto, nuestro campo de desplazamientos tiene la siguiente expresión: [math]\vec u=-\frac{x^3}{20}\vec i + \frac{y}{10}\vec j[/math]

5 Visualización del campo de vectores de desplazamiento

Ahora que disponemos de de la expresión del campo de desplazamientos, podemos representar dicho campo vectorial en los puntos del mallado del sólido. De esta manera, al observar el campo de vectores podremos intuir a priori cómo será la forma final de nuestra placa tras sufrir la deformación producida por el campo de fuerzas [math] \vec F [/math] que actúa sobre ella, y cuya expresión conoceremos más adelante.

A continuación, tenemos la imagen de la representación del campo [math]\vec u(x,y)[/math] proporcionada por el código de Matlab expuesto al final de esta sección. En esta imagen, observamos en color verde el mallado de la placa y en color rojo, el campo de vectores de desplazamiento aplicados en cada punto de la malla del sólido. Tal y como se ha mencionado anteriormente, ya nos podemos hacer una idea de la deformación que sufrirá la placa. Del mismo modo, comprobamos una de las condiciones del campo de desplazamientos; ya que observamos que los vectores aplicados en los puntos \((x,0)\); es decir, los vectores aplicados en el eje \(x\) no sufren ningún desplazamiento en el sentido de [math]\vec j[/math].

El código de Matlab empleado para la representación del campo vectorial de desplazamientos es el siguiente:

right

pm=1/10;
x=[-2:pm:2];
y=[0:pm:4];
%Definición de una nueva función como la primera componente del gradiente del campo de desplazamiento.
UX=inline('(-1*x.^3)./20','x','y');
%Definición de una nueva función como la segunda componente del gradiente del campo de desplazamiento.
UY=inline('y./10','x','y');
%Creación del mallado.
[X,Y]=meshgrid(x,y);
%Obtención de valores de las funciones anteriormente definidas para cada punto del mallado.
U=UX(X,Y);
V=UY(X,Y);
%Representación gráfica.
hold on
mesh(X,Y,0*X)
quiver(X,Y,U,V,'r');
axis([-3 3 0 5]);
hold off

6 Comparación antes y después del desplazamiento

En esta sección, visualizaremos la forma final de la placa tras haber sufrido la deformación. Ésto se consigue mediante la siguiente expresión: [math]\vec r(x,y) = \vec r_{0}(x,y) + \vec u(x,y)[/math] donde [math]\vec r_{0}(x,y)[/math] representa el vector de posición de los puntos de la placa antes de sufrir la deformación y [math]\vec u(x,y)[/math], es el campo de desplazamientos definido anteriormente.

Por lo tanto, el método para conocer el estado final de la placa consiste únicamente en aplicar a los diferentes puntos del mallado del sólido el desplazamiento definido por [math]\vec u(x,y)[/math]; siguiendo la fórmula antes descrita sobre el vector de posición final de los puntos de la placa [math]\vec r(x,y)[/math]. De esta manera, si analizamos los puntos del mallado inicial de la placa en esta expresión, obtendremos sus puntos finales tras sufrir la deformación. Si describimos el vector de posición final de los puntos de la placa de la siguiente manera: [math]\vec r(x,y) = r_x\vec i + r_y\vec j[/math] podemos afirmar que las coordenadas \(x\) de los puntos finales se corresponden con [math]r_x[/math]; mientras que sus coordenadas \(y\) se corresponden con [math]r_y[/math]. De esta manera, habremos obtenido las matrices \(POSFX\) Y \(POSFY\), con las coordenadas \(x\) e \(y\) respectivamente, de los puntos de la placa final.

A continuación podemos observar una imagen de la representación del estado de la placa antes y después de la deformación. La imagen consta de tres representaciones de la placa etiquetadas con el estado analizado.

centro

Como puede apreciarse en la representación de la placa ya deformada, hay un desplazamiento lateral de compresión, compuesta, además, por un desplazamiento de expansión vertical; únicamente en el sentido positivo del eje \(y\). Esta apreciación es razonable y lógica si se comprueba la representación del campo de vectores de desplazamiento, y se consideran los signos y direcciones, definidos por [math]\vec u(x,y)[/math].

El código de Matlab empleado para obtener la representación del estado de la placa es el que se puede ver a continuación:

pm=1/10;
x=[-2:pm:2];
y=[0:pm:4];
%Definición de las funciones que describen los desplazamientos que tienen lugar en la placa tomando como incógnitas x e y 
RX=inline('x-(x.^3)./20','x','y');
RY=inline('y+y./10','x','y');
%Creación del mallado 
[X,Y]=meshgrid(x,y);
%Obtención de valores de las funciones descriptoras del desplazamiento para cada punto del mallado
POSFX=RX(X,Y);
POSFY=RY(X,Y);
%Representación gráfica
subplot(1,3,1)
surf(X,Y,0*X);
view(2)
axis([-3 3 0 5]);
title('Antes del desplazamiento');
subplot(1,3,2)
surf(POSFX,POSFY,0*POSFX);
view(2)
axis([-3 3 0 5]);
title('Después del desplazamiento');
subplot(1,3,3)
hold on
surf(X,Y,0*X);
surf(POSFX,POSFY,0*POSFX);
view(2)
axis([-3 3 0 5]);
title ('Comparación')
hold off

7 Visualización de la divergencia del campo de desplazamientos

La divergencia del campo de desplazamientos fue dada como una de las condiciones que debía cumplir el campo [math]\vec u[/math], para la obtención de la función desconocida [math]f(y)[/math]. Su expresión viende dada por:[math]\nabla\cdot\vec u= -\frac{3x^2}{20} + \frac{1}{10}[/math]

Sabemos que la divergencia es una medida del cambio de volumen local debido al desplazamiento; por lo que es necesario representar la divergencia del campo [math]\vec u[/math], como campo escalar; así como también, aquellos puntos en los que la divergencia sea mínima, nula y máxima. De esta manera, obtendremos información muy importante sobre la intensidad de la deformación en los diferentes puntos y zonas de la placa.

A continuación, tenemos la representación de la divergencia del campo, su proyección sobre el plano \(XY\) y acto seguido, la proyección de la superficie del campo sobre el plano \(ZX\); debido a la forma de la superficie y con objeto de facilitar la visualización de los puntos máximo, mínimo y nulo; de la placa.

centro

En la imagen, podemos observar que la superficie presenta de sección una parábola, mediante planos paralelos al plano \(ZX\) o planos \(y=k\), donde \(k ∈ [0,4]\). Ésto da lugar a la superficie que vemos en la imagen.

Por lo tanto, podemos afirmar que los puntos máximo, mínimo y nulo; de la divergencia del campo, se tratan de las proyecciones de las rectas máximo, mínimo y nulo; de la divergencia. Estas rectas horizontales en [math]R^3[/math] presentan como valor fijo \(x\) y \(z\); mientras que [math]y ∈ [0,4][/math].

En definitiva, visualizando tanto la superficie de la divergencia como su proyección sobre el plano horizontal, comprobamos la premisa inicial: la divergencia es una medida del cambio de volumen local debido al desplazamiento; centrándose principalmente en la zona media de la placa. Si echamos la vista atrás hacia el estado final de la placa, comprobamos lo dicho anteriormente.

La representación de la divergencia se han obtenido con el código de Matlab siguiente:

pm=1/10; 
x=[-2:pm:2];          
y=[0:pm:4];
[X,Y]=meshgrid(x,y);
%Definición de la divergencia del campo de desplazamientos u tomando como incógnitas x e y 
DIVU=inline('((-3*x.^2)/20)+1/10','x','y');
Z=DIVU(X,Y);
m=DIVU(x,0);
%Obtención del valor máximo de la divergencia del campo de desplazamientos
divmax=max(m);
%Obtención del valor mínimo de la divergencia del campo de desplazamientos
divmin=min(m);
%Obtención de las raices del polinomio de la divergencia para conocer los valores nulos de la divergencia del campo de desplazamientos
t=[-3/20,0,1/10];
raices=roots(t);
%Representación gráfica de la divergencia
subplot(1,3,1)
surf(X,Y,Z);
shading flat
xlabel('Eje x')
ylabel('Eje y')
zlabel('Eje z')
subplot(1,3,2)
pcolor(X,Y,Z)
xlabel('Eje x')
ylabel('Eje y')
%Representación gr?fica de los valores máximos, mínimos y nulos de la divergencia. 
subplot(1,3,3)
hold on
plot(x,m,'b');
axis([-3 3 -2 2]);
for x=-2:pm:2;
    n=DIVU(x,0);
    if n==divmax
        plot(x,n,'xr','markersize',10)
    end
    if n==divmin
        plot(x,n,'xg','markersize',10)
    end
end
    plot(raices(1),0,'xk','markersize',10)
    plot(raices(2),0,'xk','markersize',10)
    xlabel('Eje x')
    ylabel('Eje z')
legend('gráfica de la diverencia','punto mínimo','punto máximo','punto mínimo','punto nulo','punto nulo');
hold off

8 Obtención y visualización del módulo del rotacional

Sabemos que el rotacional de un campo vectorial presenta interpretación física; de manera que el rotacional se define como la tendencia o capacidad de rotación o giro por parte de un cuerpo. Ésto tiene sentido si consideramos un cuerpo sometido a algún esfuerzo mecánico, que provocará una deformación en el mismo. Por lo tanto, si aplicamos el rotacional del campo que define la deformación del cuerpo y éste es distinto de cero; podemos afirmar que la deformación del sólido consistirá en un giro o una rotación; como puede ser el caso de un fluido. De lo contrario, afirmaremos que nuestro sólido no sufre tal deformación.

Sin realizar ningún cálculo, si observamos el estado de nuestra placa después de la deformación no vemos ningún tipo de rotación posible. A priori, podemos afirmar que el rotacional de nuestro campo de desplazamientos [math]\vec u[/math] será nulo. Vamos a comprobar esta suposición.

Si definimos la expresión de [math]\vec u[/math] como: [math]\vec u = u_x \vec i + u_y \vec j + u_z \vec k[/math]

la expresión del rotacional de un campo vectorial definido en el sistema de coordenadas cartesianas con base ortonormal canónica [math]{\vec i;\vec j;\vec k}[/math], viene dada por:

[math] \nabla\times \mathbf{u}=\left| \begin{matrix} \vec i & \vec j & \vec k \\ & & \\ \cfrac{\partial}{\partial x} & \cfrac{\partial}{\partial y} & \cfrac{\partial}{\partial z} \\ & & \\ u_x & u_y & u_z \end{matrix}\right| [/math]

[math] \nabla\times \mathbf{u} =\left( \frac{\partial u_z}{\partial y} - \frac{\partial u_y}{\partial z}\right)\vec i + \left(\frac{\partial u_x}{\partial z} - \frac{\partial u_z}{\partial x}\right)\vec j + \left(\frac{\partial u_y}{\partial x} - \frac{\partial u_x}{\partial y}\right)\vec k [/math]

Siguiendo la expresión anterior, calculamos el rotacional: [math]\nabla\times \mathbf{u} =\left( \frac{\partial 0}{\partial y} - \frac{\partial \frac{y}{10}}{\partial z}\right)\vec i + \left(\frac{\partial (-\frac{x^3}{20})}{\partial z} - \frac{\partial 0}{\partial x}\right)\vec j + \left(\frac{\partial \frac{y}{10}}{\partial x} - \frac{\partial (-\frac{x^3}{20})}{\partial y}\right)\vec k [/math]

Por lo tanto, comprobamos que: [math]\nabla\times \mathbf{u} = 0[/math]

De esta manera, comprobamos la suposición realizada al principio.

Por otro lado, el hecho de que el rotacional de un campo vectorial sea nulo nos puede aportar información muy importante sobre el campo vectorial en cuestión. Según el teorema del potencial, si el rotacional de un campo vectorial es nulo, podemos afirmar que dicho campo es el gradiente de un campo escalar; su función potencial. Además, podemos asegurar que dicho campo vectorial es un campo conservativo; de tal manera que, teóricamente, el cálculo de la circulación de dicho campo sobre una curva \(α\) se reduce analizando la función potencial del campo en el punto final e inicial de la curva, y realizando la diferencia de ambos valores. Por lo tanto:

si parametrizamos la curva \(α\) definida en [math]R^2[/math]: [math]α(t)= \begin{cases} x = x(t) \\ y = y(t) \end{cases}[/math] donde: [math]t ∈ [a,b][/math]

La circulación del campo [math]\vec u[/math] sobre la curva \(α\) viene dada por

[math]\int_α \mathbf{\vec u}\cdot\,d\mathbf{\vec r} = \int_a^b \mathbf{\vec u}(\mathbf{\vec r}(t))\cdot\mathbf{r}'(t)\,dt.= f(a_f,b_f)-f(a_i,b_i)[/math] siendo \(f\) la función potencial de [math]\vec u[/math].

Para demostrar esta premisa, vamos a calcular la circulación del campo de las dos formas; mediante la expresión tradicional de la circulación de un campo sobre una curva y mediante el teorema del potencial.

En primer lugar, calculamos la función potencial del campo [math]\vec u[/math]: [math]\vec u=-\frac{x^3}{20}\vec i + \frac{y}{10}\vec j = \frac{\partial f }{\partial x}\vec i+\frac{\partial f}{\partial y}\vec j[/math]

[math] \begin{cases} \frac{\partial f }{\partial x}= -\frac{x^3}{20} \\ \frac{\partial f}{\partial y}= \frac{y}{10} \end{cases}[/math]

[math]f= \int -\frac{x^3}{20}\,dx= -\frac{x^4}{80} + φ(y)[/math]

[math]\frac{\partial f}{\partial y}= φ'(y) = \frac{y}{10}[/math]

[math]f= \int \frac{y}{10}\,dy= \frac{y^2}{20} + C[/math] donde \(C\) es una constante de integración; tal que [math]C ∈ R[/math]

Por lo tanto, la función potencial de [math]\vec u[/math] presenta la siguiente expresión: [math]f(x,y) = -\frac{x^4}{80} + \frac{y^2}{20} + C[/math]

Ahora, escogemos una curva de expresión sencilla para la comprobación del teorema. Escogemos una recta horizontal cuyo punto inicial sea \((x,y)=(2,1)\) y punto final, \((x,y)=(3,1)\). si parametrizamos dicha curva, obtenemos: [math]α(t)= \begin{cases} x= t \\ y= 1\end{cases}[/math]

De tal manera que: [math]\vec r(t)= t \vec i + \vec j[/math] donde [math]t ∈ [2,3][/math]

Por lo tanto: [math]\int_α \vec u\,d\vec r= \int_t \vec u(\vec r(t))\cdot\frac{d\vec r}{dt}\,dt = \int_2^3 (-\frac{t^3}{20})\,dt = [-\frac{t^4}{80}]_2^3 = -\frac{81}{80} + \frac{16}{80} = -\frac{13}{16}[/math]

Aplicando el teorema del potencial: [math]\int_α \vec u\,d\vec r= f(3,1) - f(2,1) = -\frac{81}{16} + \frac{16}{80} = -\frac{13}{16}[/math]

Por lo tanto, dado que los resultados coinciden, podemos afirmar que el campo [math]\vec u[/math] es un campo conservativo; de tal manera que [math]∃ f; \nabla f = \vec u[/math].

A continuación, se muestran las imágenes resultantes del código de Matlab expuesto al final de la sección. En esta imagen, observamos la placa antes de la deformación en [math]R^3[/math] y proyectada sobre el plano horizontal. Esto se debe a que la función del rotacional, al ser nula y analizada con las matrices del mallado inicial de la placa; da una matriz de ceros.

centro

El código de Matlab empleado para obtener los gráficos es el siguiente:

pm=1/10;
x=[-2:pm:2];
y=[0:pm:4];
[X,Y]=meshgrid(x,y);
%Obtención de un vector fila de dos componentes el número de filas y el número de columnas de la matriz X 
[m,n]=size(X);
%Definición de la función del módulo del rotacional
MODROT=inline('0*x','x','y');
%Obtención de valores de MODROT para cada punto del mallado
Z=MODROT(X,Y);
%Representación gráfica del módulo del rotacional
subplot(1,2,1)
surf(X,Y,Z);
subplot(1,2,2)
surf(X,Y,Z);
view(2)
axis([-3 3 0 5])

9 Obtención y visualización de la proyección de las tensiones en las tres direcciones del espacio [math]{\vec i;\vec j;\vec k}[/math]

Definamos [math]\epsilon(\vec u)=(\nabla \vec u + \nabla \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]\sigma_{ij}[/math] a través de la fórmula: [math] \sigma_{ij}=\lambda \nabla \cdot \vec u \delta_{ij} + 2\mu \epsilon_{ij} [/math]

[math] \sigma =\lambda \nabla \cdot \vec u 1 + 2\mu \epsilon [/math] donde [math]\lambda[/math] y [math]\mu[/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]\lambda=\mu=1[/math], vamos a representar las tensiones normales a la dirección que marca el eje [math]\vec i [/math], es decir [math] \vec i \cdot \sigma \cdot \vec i [/math], las tensiones normales en la dirección que marca el eje [math]\vec j[/math], es decir [math] \vec j \cdot \sigma \cdot \vec j [/math] y las correspondientes al eje [math]\vec k[/math], es decir [math] \vec k \cdot \sigma \cdot \vec k [/math].

Antes de calcular la proyección de las tensiones; es necesario calcular el tensor de deformaciones, como la parte simétrica de [math]\nabla\vec u[/math]. Mediante cálculo matricial: [math]\epsilon(\vec u)=\frac{(\nabla \vec u + \nabla \vec u^t)}{2} = \frac{\begin{pmatrix} -\frac{3x^2}{20} & 0 & 0 \\ 0 & \frac{1}{10} & 0 \\ 0 & 0 & 0 \end{pmatrix} + \begin{pmatrix} -\frac{3x^2}{20} & 0 & 0 \\ 0 & \frac{1}{10} & 0 \\ 0 & 0 & 0 \end{pmatrix}}{2} = \begin{pmatrix} -\frac{3x^2}{20} & 0 & 0 \\ 0 & \frac{1}{10} & 0 \\ 0 & 0 & 0 \end{pmatrix}[/math]

El tensor de tensiones viene dado matricialmente mendiante la siguiente expresión: [math]\sigma = \nabla \cdot \vec u 1 + 2\epsilon = \begin{pmatrix} -\frac{3x^2}{20} + \frac{1}{10} & 0 & 0 \\ 0 & -\frac{3x^2}{20} + \frac{1}{10} & 0 \\ 0 & 0 & -\frac{3x^2}{20} + \frac{1}{10} \end{pmatrix} + 2\begin{pmatrix} -\frac{3x^2}{20} & 0 & 0 \\ 0 & \frac{1}{10} & 0 \\ 0 & 0 & 0 \end{pmatrix} = \begin{pmatrix} -\frac{9x^2}{20} + \frac{1}{10} & 0 & 0 \\ 0 & -\frac{3x^2}{20} + \frac{3}{10} & 0 \\ 0 & 0 & -\frac{3x^2}{20} + \frac{1}{10} \end{pmatrix}[/math]

Ahora que disponemos del tensor de tensiones, vamos a calcular la proyección de las tensiones en la dirección [math]{\vec i;\vec j;\vec k}[/math]. Las expresiones de las tensiones a calcular antes mencionada, muestran que las tensiones normales en la dirección [math]{\vec i;\vec j;\vec k}[/math], se corresponden con las componentes de la diagonal principal de la matriz de tensiones.

Por lo tanto, las tensiones normales en la dirección que marca el eje [math]\vec i[/math], vienen dadas por: [math] \vec i \cdot \sigma \cdot \vec i = -\frac{9x^2}{20} + \frac{1}{10}[/math]

Las tensiones normales en la dirección que marca el eje [math]\vec j[/math], vienen dadas por: [math] \vec j \cdot \sigma \cdot \vec j = -\frac{3x^2}{20} + \frac{3}{10}[/math]

Las tensiones normales en la dirección que marca el eje [math]\vec k[/math], vienen dadas por: [math] \vec k \cdot \sigma \cdot \vec k = -\frac{3x^2}{20} + \frac{1}{10}[/math]

Las tensiones son campos escalares; de tal manera que para su representación los tratamos como tal. A continuación, tenemos una imagen que muestra la proyección de las superficies de los campos de tensiones sobre el plano horizontal; dotadas de graduación del color en función del valor de dichas tensiones. En las tres representaciones, observamos que las tensiones son mayores en la sección media de la placa.

centro

El código de Matlab empleado para la representación de las tensiones es el siguiente:

pm=1/10;
x=[-2:pm:2];
y=[0:pm:4];
[X,Y]=meshgrid(x,y);
%Definición de las funciones correspondientes a cada componente de la función de tensiones tomando como incógnitas x e y en todos los casos
TX=inline('(-9*x.^2)./20+1/10','x','y');
TY=inline('(-3*x.^2)./20+3/10','x','y');
TZ=inline('(-3*x.^2)./20+1/10','x','y');
%Obtención de valores de las funciones TX, TY y TZ para cada punto del mallado 
U=TX(X,Y);
V=TY(X,Y);
W=TZ(X,Y);
%Representación gráfica de las tensiones en la dirección del eje i, del eje j y del eje k 
subplot(1,3,1)
pcolor(X,Y,U)
colorbar
axis image
title('Tensiones normales en la dirección que marca el eje i')
subplot(1,3,2)
pcolor(X,Y,V)
colorbar
axis image
title('Tensiones normales en la dirección que marca el eje j')
subplot(1,3,3)
pcolor(X,Y,W)
colorbar
axis image
title('Tensiones normales en la dirección que marca el eje k')

10 Tensiones tangenciales

Además de las tensiones calculadas anteriormente, podemos calcular las tensiones tangenciales respecto al plano ortogonal a \(\vec i\); cuya expresión viene da por:[math]|\sigma \cdot \vec i-(\vec i \cdot \sigma \cdot \vec i) \vec i|[/math]

Procediendo matricialmente y siguiendo la expresión anterior: [math]|\sigma \cdot \vec i-(\vec i \cdot \sigma \cdot \vec i) \vec i| = |\begin{pmatrix} -\frac{9x^2}{20} + \frac{1}{10} & 0 & 0 \\ 0 & -\frac{3x^2}{20} + \frac{3}{10} & 0 \\ 0 & 0 & -\frac{3x^2}{20} + \frac{1}{10} \end{pmatrix} \cdot \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix} - (-\frac{9x^2}{20} + \frac{1}{10}) \cdot \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}| = |(-\frac{9x^2}{20} + \frac{1}{10}; 0 ; 0) - (-\frac{9x^2}{20} + \frac{1}{10}; 0 ; 0)| = 0[/math]

Por lo tanto, afirmamos que no hay tensiones tangenciales respecto al plano ortogonal a \(\vec i\)

11 Tensión de Von Mises y autovalores

La tensión de Von Mises se define por la fórmula: [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). Se trata de una magnitud escalar que se suele usar como indicador para saber cuándo un material inicia un comportamiento plástico (y no elástico puro).

En esta ocasión, no es posible realizar los cálculos analíticamente; ya que es necesario calcular los autovalores y emplear la expresión anterior para cada punto \((x,y)\) del mallado de la placa. Sin embargo, mediante un programa de código Matlaba no es muy difícil; ya que elaborando un bucle, se realizarán los diferentes cálculos para cada iteración del bucle.

A continuación, tenemos la imagen proporcionada por el código de Matlaba expuesto al final de esta sección.

centro

En esta imagen, se puede observar la superficie que representa la tensión de Von Misses, tratado como un campo escalar; de tal manera que para cada punto \((x,y)\) del mallado de la placa le hace corresponder un escalar o altura [math]\sigma_{VM}[/math]. Mediante el bucle antes descrito, se ha ido almacenando el valor de [math]\sigma_{VM}[/math] de cada iteración, en una matriz que presentará las mismas dimensiones que las matrices del mallado de la placa.

Por otro lado, dada la forma de la superficie, ésta se ha proyectado sobre el plano \(XZ\). De esta manera, los puntos máximos que se pueden ver en la proyección representan la proyección de las rectas máximas de la superficie sobre el plano \(XZ\). El valor máximo de [math]\sigma_{VM}[/math] se aporta en forma de mensaje, en el código, mediante el comando fprintf.

El código Matlab elaborado para la representación es el siguiente:

pm=1/10;
x=[-2:pm:2];
y=[0:pm:4];
[X,Y]=meshgrid(x,y);
%Obtención del numero de componentes del vector discretizado 
t=length(x);
TX=inline('(-9*x.^2)./20+1/10','x','y');
TY=inline('(-3*x.^2)./20+3/10','x','y');
TZ=inline('(-3*x.^2)./20+1/10','x','y');
%Obtención de autovalores de la tensión de Von Mises
for i=1:t
    for j=1:t
    %Obtención de valores para cada componente del mallado para las funciones de las componentes de la función de tensión de Von Mises    
    U=TX(X(i,j),Y(i,j));
    V=TY(X(i,j),Y(i,j));
    W=TZ(X(i,j),Y(i,j));
    %Creación de un vector cuyas componentes sean los vectores en los que se han almacenado los valores de las funciones estudiadas
    vec=[U V W]; 
    %Diagonalización de ese vector
    tt=diag(vec);
    %Obtención de los autovalores de la tensión de Von Mises
    auto=eig(tt);
    %Cálculo de la tensión de Von Mises para cada componente del mallado 
    VM=sqrt(((auto(1)-auto(2))^2+((auto(2)-auto(3))^2+((auto(3)-auto(1))^2))*1/2));
    Z(i,j)=VM;
    end
end
%Representación gráfica de la tensión de Von Mises para cada punto del mallado
subplot(1,3,1)
surf(X,Y,Z)
xlabel('Eje x')
ylabel('Eje y')
zlabel('Eje z')
subplot(1,3,2)
pcolor(X,Y,Z)
xlabel('Eje x')
ylabel('Eje y')
colorbar
axis image
%Estudio del máximo valor que alcanza la tensión de Von Mises y en qué punto
subplot(1,3,3)
hold on
m=Z(1,1:length(Z));
plot(x,m,'b');
maxZ=max(m);
for k=1:length(m);
    if m(k)==maxZ
        plot(x(k),maxZ,'xr','markersize',10)
    end
end
axis([-3 3 -0 2]);
xlabel('Eje x')
ylabel('Eje z')
legend('Proyección sobre el plano XOZ','Punto máximo')
hold off
fprintf('El máximo valor de la tensión de Von Mises es %f\n',maxZ)

El máximo valor de la tensión de Von Mises es 1.562050

12 Obtención y visualización del campo de fuerzas

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 = -\nabla \cdot \sigma = -\frac{\partial \sigma_{ji}}{\partial x_{j}} \vec e_{i} [/math].

Si analizamos la expresión anterior, se trata de derivar cada componente de cada fila con respecto a \(x\), \(y\), \(z\); respectivamente. De tal manera que las derivadas de la primera fila se corresponderán con la componente [math]\vec i[/math] del campo de fuerzas; las de la segunda fila, la componente [math]\vec j[/math]; y las de la tercera fila, la componente [math]\vec k[/math].

Realizando los diferentes cálculos: [math]\frac{\partial \sigma_{j1}}{\partial x_{j}} \vec e_{1} = (\frac{\partial (-\frac{9x^2}{20} + \frac{1}{10})}{\partial x} + \frac{\partial (0)}{\partial y} + \frac{\partial (0)}{\partial z})\vec i = -\frac{9x}{10} \vec i[/math]

[math]\frac{\partial \sigma_{j2}}{\partial x_{j}} \vec e_{2} = (\frac{\partial (0)}{\partial x} + \frac{\partial (-\frac{3x^2}{20} + \frac{3}{10})}{\partial y} + \frac{\partial (0)}{\partial z})\vec j = 0[/math]

[math]\frac{\partial \sigma_{j3}}{\partial x_{j}} \vec e_{3} = (\frac{\partial (0)}{\partial x} + \frac{\partial (0)}{\partial y} + \frac{\partial (-\frac{3x^2}{20} + \frac{1}{10})}{\partial z})\vec k = 0[/math]

Por lo tanto: [math] \vec F = -\frac{\partial \sigma_{ji}}{\partial x_{j}} \vec e_{i} = \frac{9x}{10} \vec i[/math]

A continuación, tenemos la imagen de la representación del campo de fuerzas [math] \vec F[/math].

El código Matlab empleado para la representación del campo de fuerzas es el siguiente:

right

pm=1/10;
x=[-2:pm:2];
y=[0:pm:4];
[X,Y]=meshgrid(x,y);
%Definición que la función de la componente i del campo de fuerzas tomando como incógnitas x e y
FX=inline('(9*x)./10','x','y');
%Obtención de los valores de FX para cada punto del mallado
U=FX(X,Y);
%Obtención de las dimensiones del vector en el que se han almacenado los valores que toma FX para cada punto del mallado 
[m,n]=size(U);
%Creación de una matriz nula de las mismas dimensiones que U
V=zeros(m,n);
%Representación gráfica del efecto del campo de fuerzas sobre la placa
hold on
mesh(X,Y,0*X)
quiver(X,Y,U,V,'k')
hold off

13 Cálculo de la masa total de la placa

Supongamos que la densidad de la placa es [math] d(x,y) = 1 + e^{-|x|/(y+1)^2} [/math]. Al disponer de la función densidad de la placa, podemos calcular la masa total aproximando la integral correspondiente numéricamente.

A continuación, tenemos la representación de la función densidad en [math]R^3[/math]; así como también, su proyección sobre el plano \(XY\).

centro

El código Matlab empleado para su representación es el siguiente:

pm=1/10;
x=-2:pm:2;
y=0:pm:4;
[X,Y]=meshgrid(x,y);
d=inline('(1+exp(-abs(x)./(y+1).^2))','x','y');
Z=d(X,Y);
subplot(1,2,1)
surf(X,Y,Z)
shading flat
title('Superficie de la función densidad');
xlabel('Eje x')
ylabel('Eje y')
zlabel('Eje z')
subplot(1,2,2)
pcolor(X,Y,Z)
shading flat
colorbar
axis([-3 3 0 5]);
title('Proyección de la superficie densidad sobre el plano XY');
xlabel('Eje x')
ylabel('Eje y')

Para el cálculo de la masa de la placa, empleamos el código Matlab que se muestra a continuación, aportando el valor de la masa de la placa mediante el comando fprintf:

pm=1/10;
x=-2:pm:2;
y=0:pm:4;
%Obtención de la longitud de x - 1 para tener así el número de
%subintervalos. 
N1=length(x)-1; 
N2=N1; 
%Definición de los límites del intervalo
a=-2; 
b=2; 
c=0; 
d=4;      
%Definición de dos valores que indican la diferencia a tomar entre punto y punto en los intervalos que se definen a continuación como u y v 
h1=(b-a)/N1; 
h2=(d-c)/N2;
X=a:h1:b; 
Y=c:h2:d;
%Creación del mallado a partir de los intervalos u y v
[XX,YY]=meshgrid(X,Y);     
%Definición de la función de densidad en función de las variables XX e YY
D=(1+exp((-abs(XX))./(YY+1).^2));    
P1=ones(N1+1,1);                 
P1(1)=1/2; 
P1(N1+1)=1/2;
P2=ones(N2+1,1);                 
P2(1)=1/2; 
P2(N2+1)=1/2;
%Cálculo de la masa total conociendo h1, h2, P2, D y P1
M=h1*h2*rot90(P2)*D*P1;
fprintf('La masa total de la placa es %f\n', M);

La masa total de la placa es 29.474167 kg