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	Miguel Peña Martín de Prado
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)\).
• Los desplazamientos [math]\vec u(x,y)[/math] producidos por la acción de una fuerza [math]\vec{F}[/math] 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

Se dibuja el mallado que represente los puntos interiores del sólido. Para ello, se toman los ejes en el rectángulo \((x,y) ∈ [−3,3] × [0,5]\) y como paso de muestreo [math]h=\frac{1}{10}[/math] para las variables \(x\) e \(y\).

Mallado de los puntos interiores del sólido

%Paso de muestreo
pm=1/10; 
%Discretización de los 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
mesh(X,Y,0*X)
%Ajuste de la gráfica 
axis([-3,3,0,5]);
view(2);

2 Visualización del campo escalar: Temperatura

La temperatura del sólido viene dada por la función [math]T(x,y)=(x+2)^2+(y+2)^2[/math]. Para representar dicho campo escalar se dibujará tanto la superficie que genera como sus curvas de nivel. En última instancia, se señalará en qué punto del mallado la temperatura es máxima a partir de los gráficos obtenidos.

Campo escalar: Temperatura

Campo escalar: Temperatura

%Paso de muestreo
pm=1/10; 
%Discretización de los parámetros de la superficie
x=[-2:pm:2];          
y=[0:pm:4];
%Creación del mallado
[X,Y]=meshgrid(x,y);
%Definición de la función Temperatura
T=inline('((x+2).^2)+((y+2).^2)','x','y');
%Obtención de la matriz con los valores de T(x,y) sobre los puntos de la retícula anterior
Z=T(X,Y);
%Creación de subventanas y representación del campo escalar y de sus curvas de nivel
subplot(1,2,1);
surfc(X,Y,Z);
shading flat
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

Como puede observarse en ambos gráficos, la temperatura es máxima en el punto \((2,4)\). Esto también es visible si analizamos la fórmula de la temperatura, donde podemos ver que, cuanto mayor sea el valor de \(x\) e \(y\), mayor será el valor de la temperatura.

3 Visualización del campo vectorial: Gradiente de Temperatura

Puesto que el campo del que se desea calcular el gradiente es un campo escalar, el gradiente de dicho campo será un campo vectorial.
En primer lugar, se procede al cálculo del gradiente:

[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]

Por último, se dibuja el \(∇T\) para observar gráficamente que es ortogonal a las curvas de nivel del campo de temperaturas.

Campo vectorial: Gradiente de Temperatura

%Paso de muestreo
pm=1/10;
%Discretización de los parámetros de la superficie
x=[-2:pm:2];          
y=[0:pm:4];
%Creación del mallado
[X,Y]=meshgrid(x,y);
%Definición de la función Temperatura y de las componentes en la dirección de i y de j del gradiente
T=inline('((x+2).^2)+((y+2).^2)','x','y');
TX=inline('2*(x+2)','x','y');
TY=inline('2*(y+2)','x','y');
%Obtención de la matriz con los valores de T(x,y) sobre los puntos de la retícula anterior
Z=T(X,Y);
%Dibujo de las curvas de nivel del campo de temperaturas
contour(X,Y,Z);
colorbar
hold on
%Nueva discretización para representar el campo vectorial Gradiente de Temperatura
xx=[-2:0.25:2];
yy=[0:0.25:4];
%Creación del mallado
[XX,YY]=meshgrid(xx,yy);
%Obtención de la matriz con los valores de TX(x,y) y TY(x,y) sobre los puntos de la retícula anterior
U=TX(XX,YY);
V=TY(XX,YY);
%Representación de los vectores del campo Gradiente de Temperatura
quiver(XX,YY,U,V);
axis equal
hold off

4 Cálculo del campo vectorial: Campo de desplazamientos

Consideramos ahora el campo de desplazamientos [math]\vec{u}[/math] cuya expresión viene dada por:

[math]\vec u(x,y)=-\frac{x^3}{20}\vec i+f(y)\vec j[/math]

Este campo debe calcularse sabiendo:

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

Teniendo esto en cuenta, calculamos el valor \(f(y)\) a partir del campo [math]\vec{u}[/math] y de su divergencia, que nos quedará en función de una constante.
Sabemos que [math]\nabla\cdot\vec{u}=\sum\frac{\partial{u}_{i}}{\partial {x}_{i}}=\frac{\partial{u}_{1}}{\partial{x}}+\frac{\partial{u}_{2}}{\partial{y}}=-\frac{3x^2}{20}+\frac{1}{10}[/math]
Por comparación con la expresión del campo [math]\vec{u}[/math]:
[math]\frac{\partial{u}_{1}}{\partial{x}}=-\frac{3x^2}{20}[/math]
[math]\frac{\partial{u}_{2}}{\partial{y}}=\frac{1}{10}[/math]
Integrando el primer sumando de la divergencia respecto de \(x\), hallaremos la componente en la dirección [math]\vec{i}[/math] del campo de desplazamientos:
[math]\int\frac{-3x^2}{20}dx=-\frac{x^3}{20}+C[/math]
Comparando el resultado con la primera componente del campo de desplazamientos, observamos que \(C=0\)
Repetimos el mismo proceso para hallar la \(f(y)\). En este caso, integramos el segundo sumando de la divergencia respecto de \(y\):
[math]\int\frac{1}{10}dy=\frac{y}{10}+C[/math]
Comparando el resultado con la componente en la dirección de [math]\vec{j}[/math] del campo de desplazamientos, observamos que [math]f(y)=\frac{y}{10}+C[/math]
Haciendo uso de la otra condición dada en el enunciado, sabemos que \(C=0\), para que al introducir \(y=0\) en el campo [math]\vec{u}[/math], los puntos del eje de abscisas no sufran desplazamiento en la dirección [math]\vec{j}[/math]:
[math]\vec{u}=-\frac{x^3}{20}\vec{e}_{1}+f(0)\vec{e}_{2}=-\frac{x^3}{20}\vec{e}_{1}+C\vec{e}_{2}=-\frac{x^3}{20}\vec{e}_{1}[/math]
Por lo tanto, el campo de desplazamientos [math]\vec{u}[/math] tiene la siguiente expresión:

[math]\vec u(x,y)=-\frac{x^3}{20}\vec i+\frac{y}{10}\vec j[/math]

5 Visualización del campo vectorial: Campo de desplazamientos

A partir de la expresión del campo de desplazamientos obtenida en el apartado anterior, se dibuja dicho campo de vectores en los puntos del mallado del sólido

Campo vectorial: Campo de desplazamientos

%Paso de muestreo
pm=1/10;
%Discretización de los parámetros de la superficie
x=-2:pm:2;
y=0:pm:4;
%Creación del mallado
[X,Y]=meshgrid(x,y);
%Definición de las componentes en la dirección de i y de j del campo de desplazamientos 
UX=inline('(-1*x.^3)./20','x','y');
UY=inline('y./10','x','y');
%Obtención de las matrices con los valores de UX(x,y) y UY(x,y) sobre los puntos de la retícula anterior
U=UX(X,Y);
V=UY(X,Y);
hold on
%Dibujo de la malla
mesh(X,Y,0*X)
%Representación de los vectores del campo de desplazamientos
quiver(X,Y,U,V,'r');
axis([-3 3 0 5]);
hold off

6 Visualización del sólido antes y después del desplazamiento

Una vez conocidas tanto la posición inicial como el desplazamiento que sufre cada punto de la malla, puede determinarse su posición final a través de la expresión:

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

Conocida la posición final de cada punto del mallado, se dibuja el sólido antes y después del desplazamiento para poder apreciar el efecto que tiene dicho campo de vectores [math]\vec{u}[/math] sobre los puntos del sólido:

Sólido antes y después del desplazamiento

%Paso de muestreo
pm=1/10;
%Discretización de los parámetros de la superficie
x=-2:pm:2;
y=0:pm:4;
%Creación del mallado
[X,Y]=meshgrid(x,y);
%Definición de las componentes en la dirección de i y de j del vector posición de cada punto de la placa después de la deformación
RX=inline('x-(x.^3)./20','x','y');
RY=inline('y+y./10','x','y');
%Obtención de las matrices con los valores de RX(x,y) y RY(x,y) sobre los puntos de la retícula anterior
POSFX=RX(X,Y);
POSFY=RY(X,Y);
subplot(1,3,1)
%Dibujo del sólido antes del desplazamiento
surf(X,Y,0*X);
view(2)
axis([-3 3 0 5]);
title('Antes del desplazamiento');
subplot(1,3,2)
%Dibujo del sólido después del desplazamiento
surf(POSFX,POSFY,0*POSFX);
view(2)
axis([-3 3 0 5]);
title('Después del desplazamiento');
subplot(1,3,3)
hold on
%Dibujo del sólido antes y después del desplazamiento
surf(X,Y,0*X);
surf(POSFX,POSFY,0*POSFX);
view(2)
axis([-3 3 0 5]);
title ('Comparación')
hold off

Como puede observarse en los gráficos, los puntos que sufren un mayor desplazamiento son aquellos más cercanos a los bordes superior, izquierdo y derecho del sólido, mientras que los puntos más próximos a la parte central del borde inferior apenas sufren desplazamiento.

7 Visualización del campo escalar: Divergencia del campo de desplazamientos

A partir de la expresión de la [math]\nabla\cdot\vec{u}[/math], se desea representar dicho campo escalar:

[math]\nabla\cdot\vec{u}=-\frac{3x^2}{20}+\frac{1}{10}[/math]

Visualización del campo escalar: Divergencia

%Paso de muestreo
pm=1/10; 
%Discretización de los parámetros de la superficie
x=-2:pm:2;          
y=0:pm:4;
%Creación del mallado
[X,Y]=meshgrid(x,y);
%Definición de la función divergencia de u
DIVU=inline('((-3*x.^2)/20)+1/10','x','y');
%Obtención de la matriz con los valores de DIVU(x,y) sobre los puntos de la retícula anterior
Z=DIVU(X,Y);
m=DIVU(x,0);
%Hallamos el valor máximo de la divergencia
divmax=max(m);
%Hallamos el valor mínimo de la divergencia
divmin=min(m);
%Hallamos las raíces del polinomio de la divergencia para saber qué valores toma x cuando la divergencia es nula
t=[-3/20,0,1/10];
raices=roots(t);
subplot(1,3,1)
%Dibujo de la superficie que representa la divergencia
surf(X,Y,Z);
shading flat
xlabel('Eje x')
ylabel('Eje y')
zlabel('Eje z')
subplot(1,3,2)
%Dibujo de la superficie que representa la divergencia y de sus curvas de nivel visto cuando z tiende a infinito
pcolor(X,Y,Z)
axis equal
shading flat
hold on
contour(X,Y,Z,'k')
colorbar
xlabel('Eje x')
ylabel('Eje y')
hold off
subplot(1,3,3)
hold on
%Dibujo de la proyección sobre el plano XOZ de la superficie que representa la divergencia
plot(x,m,'b');
axis([-3 3 -2 2]);
%Bucle para pintar los puntos donde la divergencia es máxima y mínima
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
%Dibujo de los puntos donde la divergencia es nula
    plot(raices(1),0,'xk','markersize',10)
    plot(raices(2),0,'xk','markersize',10)
    xlabel('Eje x')
    ylabel('Eje z')
legend('Proyección de la diverencia sobre el plano XOZ','punto mínimo','punto máximo','punto mínimo','punto nulo','punto nulo');
hold off

Visualización del campo escalar: Divergencia

Proyección de la Divergencia sobre el plano XOZ

Como consecuencia de que la expresión de la [math]\nabla\cdot\vec{u}[/math] dependa únicamente de \(x\), ha sido posible determinar los puntos en los que la divergencia es máxima, mínima y nula gráficamente mediante la proyección de dicho campo sobre el plano \(XOZ\).

7.1 Cálculo de los puntos donde la Divergencia es máxima, mínima y nula

También es posible determinar analíticamente los puntos en los que la divergencia de [math]\vec{u}[/math] es máxima, mínima y nula. Al tratarse de una función que depende únicamente de una variable \(x\), el proceso a seguir es el siguiente:

• Puntos en los que la divergencia es máxima:

[math]\frac{d\nabla\cdot\vec{u}}{dx}=-\frac{6x}{20}=-\frac{3x}{10}[/math]

Igualando esta expresión a 0 obtenemos la coordenada \(x\) de los puntos en los que la divergencia es máxima o mínima:

[math]0=-\frac{3x}{10}[/math]

[math]0=x[/math]

Para [math]x=0[/math] la divergencia toma el valor [math]\frac{1}{10}[/math]. Puesto que para valores de x cercanos a 0, tanto por la izquierda como por la derecha, la divergencia toma valores menores, podemos decir que hemos encontrado un máximo de la función.

• Puntos en los que la divergencia es mínima:
  

Debemos estudiar también los extremos de la función. Al tratarse de una función par: \(f(-x)=f(x)\), estudiaremos el valor de la divergencia tan solo para [math]x=2[/math]:

-Para [math]x=2[/math]:

[math]\nabla\cdot\vec{u}=-\frac{12}{20}+\frac{1}{10}=-\frac{1}{2}[/math]

Por lo tanto, para [math]x=±2[/math] la función alcanza su valor mínimo.

• Puntos en los que la divergencia es nula:
  

Para hallar el punto donde la divergencia es nula procedemos a igualar la expresión de la divergencia a 0:

[math]\nabla\cdot\vec{u}=-\frac{3x^2}{20}+\frac{1}{10}=0[/math]

[math]\frac{3x^2}{20}=\frac{1}{10}[/math]

[math]3x^2=2[/math]

[math]x=±\sqrt{\frac{2}{3}}[/math]

Por lo tanto, para [math]x=±\sqrt{\frac{2}{3}}[/math] la función toma un valor nulo.

7.2 Interpretación de la Divergencia

La divergencia es una medida del cambio de volumen local debido al desplazamiento. Como puede observarse en las gráficas del apartado 6, el sólido sufre una compresión en la dirección horizontal y una tracción en la dirección vertical debidas al campo de desplazamientos [math]\vec{u}[/math].

• En el rectángulo [math][-\sqrt{\frac{2}{3}},\sqrt{\frac{2}{3}}]×[0,4][/math], la [math]\nabla\cdot\vec{u}[/math] es positiva, lo que indica que la compresión es mayor que la tracción, es decir, en ese rectángulo, el espesor del sólido aumenta provocando ese aumento localizado de volumen que aparece representado en el gráfico anterior.
• En el dominio [math][-2,-\sqrt{\frac{2}{3}}]×[0,4][/math] U [math][\sqrt{\frac{2}{3}},2]×[0,4][/math], la [math]\nabla\cdot\vec{u}[/math] es negativa, lo que indica que la tracción es mayor que la compresión, es decir, en ese dominio, el espesor del sólido disminuye provocando esa disminución localizada de volumen que aparece representada en el gráfico anterior.

8 Cálculo del campo vectorial: Rotacional del campo de desplazamientos

8.1 Cálculo del Rotacional

A partir de la expresión del campo de desplazamientos [math]\vec{u}[/math], se desea calcular el módulo del rotacional de dicho campo. Para calcular el [math]|\nabla×\vec{u}|[/math] se procede de la siguiente manera:

[math]|\nabla×\vec{u}|=|\det{\begin{pmatrix} \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{pmatrix}}|[/math]

En primer lugar, calculamos el rotacional y después hallamos su módulo:

[math]\nabla\times\vec {u}=\left|\begin{matrix}\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{matrix}\right|=\left|\begin{matrix}\vec{i} & \vec{j} & \vec{k} \\ \frac{ \partial}{\partial x} & \frac{\partial }{\partial y} & \frac{\partial }{\partial z}\\ -\frac{x^3}{20} & \frac{y}{10} & 0 \end{matrix}\right|=0\vec{i}+0\vec{j}+0\vec{k}=\vec{0}[/math]

Por lo tanto, el módulo del vector nulo es 0:

[math]|\vec{0}|=0[/math]

8.2 Interpretación del Rotacional

Un campo cuyo rotacional es nulo en todos los puntos del espacio se denomina irrotacional, es decir, dicho campo no induce rotación al sólido. Por lo tanto, todos los puntos sufren el mismo rotacional o, lo que es lo mismo, ningún punto sufre rotacional.
Dado que los campos con potencial escalar son irrotacionales, podemos tratar de identificar si el campo [math]\vec{u}[/math] tiene potencial escalar:
Para ello, buscamos \(Ψ(x,y,z)\) tal que [math]\nabla Ψ=\vec{u}[/math]:

[math]\nabla Ψ(x,y,z)=\frac{\partial Ψ}{\partial x}\vec{i}+\frac{\partial Ψ}{\partial y}\vec{j}+\frac{\partial Ψ}{\partial z}\vec{k}=\vec{u}=-\frac{x^{3}}{20}\vec{i}+\frac{y}{10}\vec{j}[/math]

Por lo tanto, hay que resolver este sistema de ecuaciones:

[math]\left \{ \begin{array}{rcrcrcr} \frac{\partial {Ψ}}{\partial {x}}=u_{1}=-\frac{x^{3}}{20}\\ \frac{\partial Ψ}{\partial y}=u_{2}=\frac{y}{10}\\ \frac{\partial Ψ}{\partial z}=u_{3}=0 \end{array}\right .[/math]

• Primero resolvemos la ecuación 1 integrando en x:

[math]Ψ(x,y,z)=-\frac{x^{4}}{80}+f(y,z)[/math]

• Sustituimos en la ecuación 2 y resolvemos:

[math]\frac{\partial Ψ}{\partial y}=\frac{y}{10}[/math]

[math]\frac{\partial}{\partial y}(-\frac{x^{4}}{80}+f(y,z))=\frac{y}{10}[/math]

[math]\frac{\partial f}{\partial y}=\frac{y}{10}[/math]

[math]f(y,z)=\frac{y^{2}}{20}+g(z)[/math]

• Por último, sustituimos en la ecuación 3 y resolvemos:

[math]\frac{\partial Ψ}{\partial z}=0[/math]

[math]\frac{\partial}{\partial z}(-\frac{x^{4}}{80}+\frac{y^{2}}{20}+g(z))=0[/math]

[math]\frac{dg}{dz}=0[/math]

\(g(z)=C\), donde \(C\) es constante.

Por lo tanto, el potencial escalar del campo [math]\vec{u}[/math] es:

[math]Ψ(x,y,z)=-\frac{x^{4}}{80}+\frac{y^{2}}{20}+C[/math]

9 Tensiones normales

9.1 Cálculo del tensor de deformaciones

El tensor de deformaciones [math]e(\vec{u})[/math] representa la parte simétrica del tensor [math]\nabla\vec{u}[/math], cuya expresión es:

[math]e(\vec{u})=\frac{\nabla\vec{u}+\nabla\vec{u}^{t}}{2}[/math]

Para calcularlo es necesario, en primer lugar, calcular el [math]\nabla\vec{u}[/math]:

[math]\nabla\vec{u}=\begin{pmatrix} \frac{\partial u_1}{\partial x} & \frac{\partial u_1}{\partial y} \\ \frac{\partial u_2}{\partial x} & \frac{\partial u_2}{\partial y} \end{pmatrix}=\begin{pmatrix} \frac{{-3x}^{2}}{20} & 0 \\ 0 & \frac{1}{10} \end{pmatrix}[/math]

Como podemos observar, el gradiente de [math]\vec{u}[/math] es un tensor simétrico, por lo tanto:

[math]\nabla\vec{u}+\nabla\vec{u}^{t}=2*\nabla\vec{u}[/math]

Por lo tanto, el tensor de deformaciones sería:

[math]e(\vec{u})=\frac{\nabla\vec{u}+\nabla\vec{u}^{t}}{2}=\frac{2*\nabla\vec{u}}{2}=\nabla\vec{u}=\begin{pmatrix} \frac{{-3x}^{2}}{20} & 0 \\ 0 & \frac{1}{10} \end{pmatrix}[/math]

9.2 Cálculo del tensor de tensiones [math]\sigma_{ij}[/math]

La expresión del tensor de tensiones [math]\sigma_{ij}[/math] es la siguiente:

[math]\sigma_{ij}=λ\nabla\cdot\vec{u}1+2μe(\vec{u})[/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 [math]λ=μ=1[/math], tenemos todo lo necesario para calcular [math]\sigma_{ij}[/math]:

[math]\sigma_{ij}=\nabla\cdot\vec{u}1+2e(\vec{u})=\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]

9.3 Cálculo de las tensiones normales

9.3.1 Tensiones normales en la dirección que marca el eje [math]\vec{i}[/math]

Las tensiones normales en la dirección que marca el eje [math]\vec{i}[/math] se obtienen a partir de la expresión:

[math]\vec{i}\cdot\sigma_{ij}\cdot\vec{i}=\begin{pmatrix} 1, & 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}\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}=\begin{pmatrix} -\frac{{9x}^{2}}{20}+\frac{1}{10}, & 0, & 0 \end{pmatrix}\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}=-\frac{{9x}^{2}}{20}+\frac{1}{10}[/math]

9.3.2 Tensiones normales en la dirección que marca el eje [math]\vec{j}[/math]

Las tensiones normales en la dirección que marca el eje [math]\vec{j}[/math] se obtienen a partir de la expresión:

[math]\vec{j}\cdot\sigma_{ij}\cdot\vec{j}=\begin{pmatrix} 0, & 1, & 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}\begin{pmatrix} 0 \\ 1 \\ 0 \end{pmatrix}=\begin{pmatrix} 0, & -\frac{{3x}^{2}}{20}+\frac{3}{10}, & 0 \end{pmatrix}\begin{pmatrix} 0 \\ 1 \\ 0 \end{pmatrix}=-\frac{{3x}^{2}}{20}+\frac{3}{10}[/math]

9.3.3 Tensiones normales en la dirección que marca el eje [math]\vec{k}[/math]

Las tensiones normales en la dirección que marca el eje [math]\vec{k}[/math] se obtienen a partir de la expresión:

[math]\vec{k}\cdot\sigma_{ij}\cdot\vec{k}=\begin{pmatrix} 0, & 0, & 1 \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}\begin{pmatrix} 0 \\ 0 \\ 1 \end{pmatrix}=\begin{pmatrix} 0, & 0, & -\frac{{3x}^{2}}{20}+\frac{1}{10} \end{pmatrix}\begin{pmatrix} 0 \\ 0 \\ 1 \end{pmatrix}=-\frac{{3x}^{2}}{20}+\frac{1}{10}[/math]

9.4 Visualización de las tensiones normales

A partir de las expresiones de las tensiones normales calculadas en el apartado anterior, es posible visualizar dichos campos escalares:

%Paso de muestreo
pm=1/10;
%Discretización de los parámetros de la superficie
x=-2:pm:2;
y=0:pm:4;
%Creación del mallado
[X,Y]=meshgrid(x,y);
%Definición de los campos escalares que representan las tensiones normales
NI=inline('(-9*x.^2)./20+1/10','x','y');
NJ=inline('(-3*x.^2)./20+3/10','x','y');
NK=inline('(-3*x.^2)./20+1/10','x','y');
%Obtención de las matrices con los valores de NI(x,y), NJ(x,y) y NK(x,y) sobre los puntos de la retícula anterior
U=NI(X,Y);
V=NJ(X,Y);
W=NK(X,Y);
%Visualización de las tensiones normales en la dirección que marca el Eje i
subplot(1,3,1)
pcolor(X,Y,U)
colorbar
xlabel('Eje x')
ylabel('Eje y')
axis image
title('Tensiones normales en la dirección que marca el Eje i')
%Visualización de las tensiones normales en la dirección que marca el Eje j
subplot(1,3,2)
pcolor(X,Y,V)
colorbar
xlabel('Eje x')
ylabel('Eje y')
axis image
title('Tensiones normales en la dirección que marca el Eje j')
%Visualización de las tensiones normales en la dirección que marca el Eje k
subplot(1,3,3)
pcolor(X,Y,W)
colorbar
xlabel('Eje x')
ylabel('Eje y')
axis image
title('Tensiones normales en la dirección que marca el Eje k')

Cabe destacar que al ser las tensiones normales funciones pares que dependen únicamente de la variable \(x\), se obtienen superficies que presentan simetría respecto al plano \(x=0\).

10 Cálculo de las tensiones tangenciales respecto al plano ortogonal a [math]\vec{i}[/math]

Para obtener las tensiones tangenciales respecto al plano ortogonal a [math]\vec{i}[/math] se emplea la siguiente expresión:

[math]|\sigma\cdot\vec{i}-(\vec{i}\cdot\sigma\cdot\vec{i})\cdot\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}\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}-(\begin{pmatrix} 1, & 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}\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix})\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}|=[/math]

[math]=|\begin{pmatrix} -\frac{{9x}^{2}}{20}+\frac{1}{10} \\ 0 \\ 0 \end{pmatrix}-(-\frac{{9x}^{2}}{20}+\frac{1}{10})\begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}|=|\begin{pmatrix} -\frac{{9x}^{2}}{20}+\frac{1}{10} \\ 0 \\ 0 \end{pmatrix}-\begin{pmatrix} -\frac{{9x}^{2}}{20}+\frac{1}{10} \\ 0 \\ 0 \end{pmatrix}|=|\vec{0}|=0[/math]

11 La tensión de Von Mises

11.1 Interpretación de la tensión de Von Mises

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}, \sigma_{2}[/math] y [math]\sigma_{3}[/math] son los autovalores de [math]\sigma[/math] (también conocidos como tensiones principales) cuya expresión es:

[math]\sigma=\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]

La tensión de Von Mises es una magnitud escalar que se suele usar como indicador para saber cuando un material inicia un comportamiento plástico (y no elástico puro).

11.2 Visualización de la tensión de Von Mises

Se desea representar dicho campo escalar en todos los puntos del mallado, además de determinarse los puntos donde el campo alcanza su máximo valor.

Tensión de Von Mises

%Paso de muestreo
pm=1/10;
%Discretización de los parámetros de la superficie
x=-2:pm:2;
y=0:pm:4;
%Creación del mallado
[X,Y]=meshgrid(x,y);
t=length(x);
%Definición de las componentes de la diagonal del tensor de tensiones sigma
COMP1=inline('(-9*x.^2)./20+1/10','x','y');
COMP2=inline('(-3*x.^2)./20+3/10','x','y');
COMP3=inline('(-3*x.^2)./20+1/10','x','y');
%Bucle que haya los autovalores para cada punto del mallado y calcula la tensión de Von Mises a partir de ellos
for i=1:t
    for j=1:t
    U=COMP1(X(i,j),Y(i,j));
    V=COMP2(X(i,j),Y(i,j));
    W=COMP3(X(i,j),Y(i,j));
    vec=[U V W]; 
    tt=diag(vec);
    auto=eig(tt);
    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
%Dibujo de la superficie que genera el campo escalar: Tensión de Von Mises
subplot(1,3,1)
surf(X,Y,Z)
xlabel('Eje x')
ylabel('Eje y')
zlabel('Eje z')
%Dibujo de la superficie que genera el campo escalar: Tensión de Von Mises, con vista cuando z tiende a infinito
subplot(1,3,2)
pcolor(X,Y,Z)
xlabel('Eje x')
ylabel('Eje y')
colorbar
axis image
%Estudio del máximo valor de la tensión de Von Mises
subplot(1,3,3)
hold on
m=Z(1,1:length(Z));
%Dibujo de la proyección sobre el plano XOZ de la superficie que genera la tensión de Von Mises
plot(x,m,'b');
maxZ=max(m);
%Bucle que pinta los puntos donde la tensión de Von Mises es máxima sobre la proyección
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)

Como consecuencia de que la expresión del tensor de tensiones [math]\sigma[/math] dependa únicamente de \(x\), es posible determinar los puntos donde la tensión de Von Mises alcanza los valores máximos gráficamente mediante la proyección de dicho campo sobre el plano \(XOZ\).
Como es posible observar en el gráfico superior derecho, el máximo valor de la tensión de Von Mises se alcanza para [math]x=±2[/math] y, como bien muestra el programa por pantalla, el valor de la tensión de Von Mises para dicho punto es de \(1,562050\)

12 Campo de fuerzas [math]\vec{F}[/math] que actúan sobre la placa

12.1 Cálculo del campo de fuerzas [math]\vec{F}[/math]

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]

Calculamos las tres componentes del campo de vectores:

• [math]F_1=-\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]F_2=-\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]F_3=-\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, la expresión del campo de fuerzas [math]\vec{F}[/math] es:

[math]\vec{F}=\frac{9x}{10} \vec i[/math]

12.2 Visualización del campo de fuerzas [math]\vec{F}[/math]

Visualización del campo vectorial: Campo de fuerzas

%Paso de muestreo
pm=1/10;
%Discretización de los parámetros de la superficie
x=-2:pm:2;
y=0:pm:4;
%Creación del mallado
[X,Y]=meshgrid(x,y);
%Definición de la componente en la dirección de i del campo de fuerzas por tratarse de la única componente no nula
FX=inline('(9*x)./10','x','y');
%Obtención de la matriz con las componentes en la dirección de i de los vectores del campo
U=FX(X,Y);
%Obtención de la matriz con las componentes en la dirección de j de los vectores del campo
[m,n]=size(U);
V=zeros(m,n);
hold on
%Dibujo del mallado
mesh(X,Y,0*X)
%Dibujo del campo vectorial
quiver(X,Y,U,V,'k')
hold off

12.3 Interpretación del campo de fuerzas [math]\vec{F}[/math]

Como puede verse en el gráfico, sobre los puntos del eje de ordenadas no actúa ninguna fuerza, sin embargo, a medida que nos separamos del eje de ordenadas, el módulo de la fuerza aumenta. Además, el campo de fuerzas [math]\vec{F}[/math] solo depende de \(x\) y, por lo tanto, todos los puntos del mallado con la misma coordenada \(x\), están sometidos a la misma fuerza [math]\vec{F}[/math]. Por último, cabe destacar que los vectores actúan únicamente en la dirección de [math]\vec{i}[/math] dado que las componentes del campo [math]\vec{F}[/math] en la dirección de [math]\vec{j}[/math] y de [math]\vec{k}[/math] son nulas.

13 Cálculo de la masa total de la placa

Suponiendo que la placa tiene una densidad que responde a la ecuación:

[math]d(x,y)=1+e^{\frac{-|x|}{(y+1)^{2}}}[/math]

La masa total de la placa se calcularía de la siguiente forma:

[math]\int_{A} d(x,y) dA=\int_{-2}^{2}\int_{0}^{4} d(x,y) dydx=\int_{-2}^{2}\int_{0}^{4} 1+e^{\frac{-|x|}{(y+1)^{2}}} dydx=29,4754[/math]

Sin embargo, puesto que, en esta ocasión, la primitiva es difícil de calcular, se recurre a la utilización de métodos numéricos. Por esta razón, empleamos la fórmula del trapecio para la aproximación de integrales:

%Paso de muestreo
pm=1/10;
x=-2:pm:2;
y=0:pm:4;
%length(x) nos da el número de puntos, pero le restamos 1 para tener el número de subintervalos
N1=length(x)-1; N2=N1;
%Extremos del intervalo
a=-2; b=2; c=0; d=4;             
h1=(b-a)/N1; h2=(d-c)/N2;
%Coordenadas de la partición
u=a:h1:b; v=c:h2:d;              
%Creación del mallado
[uu,vv]=meshgrid(u,v);   
%Definición de la función densidad
f=1+exp((-abs(uu))./(vv+1).^2);                
w1=ones(N1+1,1);                 
w1(1)=1/2; w1(N1+1)=1/2;
w2=ones(N2+1,1);                
w2(1)=1/2; w2(N2+1)=1/2;
result=h1*h2*rot90(w2)*f*w1

El resultado que el programa muestra por pantalla es \(29,4742\), que es un valor muy aproximado al que se obtiene analíticamente.