Visualización de campos escalares y vectoriales en elasticidad - Grupo A1

Trabajo realizado por estudiantes
Título	Visualización de campos escalares y vectoriales en elasticidad - Grupo A1
Asignatura	Teoría de Campos
Curso	2021-22
Autores	Héctor Alcañiz Poyatos Iván Calle Rodríguez Facundo Pérez Cristina Ye Wang
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

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

[math]T(x, y)=(x-3)^2+(10y)^2[/math]

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

[math]\vec{r}(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{yf(x)}{20}\vec{i}-\frac{y^2}{40}\vec{j}[/math]

donde [math]f(x)[/math] es una cierta función que no conocemos.

1 Representación de la placa

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

Con las 3 primeras líneas de código (clc, clear y close all) realizamos una limpieza de programas anteriores para evitar que el programa o los gráficos no se ejecuten correctamente.

Tras crear el paso de muestreo [math]h[/math]. y las variables [math]x[/math] e [math]y[/math] con los valores que se nos indican, utilizamos el comando [math]meshgrid ()[/math], que nos devuelve las matrices [math]X[/math] e [math]Y[/math].

A continuación, utilizamos el comando [math]mesh ()[/math] con las matrices obtenidas. Sin embargo, al ser un comando que requiere de 3 matrices ([math]X[/math], [math]Y[/math] y [math]Z[/math]) y solo neesitamos 2, aquella que se corresponde con la tercera matriz será 0, por lo que multiplicaremos por 0 una de las dos matrices existentes, por ejemplo la [math]Y[/math].

Por último, ajustaremos los ejes según los límites que nos han indicado y escribiremos el título del gráfico y de los ejes.

clc
clear
close all
% Definimos el paso de muestreo h para las variables x e y
h=1/10; 
% Definimos los parámetros que definen la placa sólida: x e y
x=[0:h:10];          
y=[-1:h:1];
% Creamos el mallado con las matrices X e Y
[X,Y]=meshgrid(x,y); 
% Dibujamos la malla de puntos. mesh () es un comando que requiere de 3 elementos de entrada, por lo que la matriz Z será nula: 0 veces la matriz Y, por ejemplo
mesh(X,Y,0*Y)
% Utilizamos el comando axis (), pues queremos que el eje de abscisas vaya desde el -0.5 hasta el 10.5 y el eje Y desde el -1.5 hasta el 1.5
% Establecemos los límites de los ejes
axis([-0.5,10.5,-1.5,1.5]);
% Escribimos el título del gráfico y los nombres de los ejes 
title('Mallado');
xlabel('Eje X');
ylabel('Eje Y');
% Con el comando view (2), visuaizamos el mallado en 2 dimensiones
view(2);

2 Curvas de nivel de la temperatura

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.

A partir de las variables [math]h[/math]., [math]x[/math] e [math]y[/math] creadas anteriormente, volvemos a utilizar el comando [math]meshgrid ()[/math] para obtener las matrcies [math]X[/math] e [math]Y[/math].

A continuación, definimos el campo de temperatura T en función de las matrices X e Y. Utilizaremos [math].[/math] cuando realicemos operaciones como elevar al cuadrado, pues el punto permite realizar esa operación a cada elemento de la matriz (la operación se realiza elemento a elemento gracias al punto).

Utilizaremos el comando [math]mesh()[/math] para representar la función de temperatura utilizando las matrices [math]X[/math], [math]Y[/math] y [math]T[/math]. Para representar las curvas de nivel, utilizaremos el comando [math]contour ()[/math], mediante el que dibujaremos 50 curvas de nivel. Como obtendremos dos gráficos, los representaremos en la misma ventana con el comando [math]subplot ()[/math].

Por último, ajustaremos los ejes según los límites indicados, escribiremos el título del gráfico y de los ejes y nos ayudaremos del comando [math]colorbar[/math] para incluir una barra de colores que permita identificar los valores del campo de temperaturas. De esta manera, los colores fríos indican un menor valor de la temperatura, mientras que los colores cálidos se corresponden con valores mayores.

Como podemos observar al ampliar el gráfico de la temperatura, el valor máximo que esta alcanza es de 149ºC y se alcanza en dos puntos: en el (0,10) y en el (10,1).

Curvas de nivel de la temperatura

clear
 close all
  % Paso de muestreo h para las variables x e y
  h = 2/10;
  % Región que ocupa la placa rectangular
  x = [0:h:10];
  y = [0:h:2];
  % Creación del mallado
 [X,Y]= meshgrid(x,y);
 % Función temperatura: T(x,y) = (x-6)^2 + (10(Y-3/2))^2
 T = (X-6).^2 + (10*(Y-3/2)).^2;
 % Escribimos el título del gráfico
 title('Curvas de nivel de la temperatura');
 % Escribimos los nombres de los ejes
 xlabel('Eje X');
 ylabel('Eje Y');
% Establecemos el límite de los ejes
 axis([-0.5,10.5,-0.5,2.5]);
 % Representación de la temperatura y las curvas de nivel
 subplot(1,2,1);
 mesh(X,Y,T);
 subplot(1,2,2);
 contour(X,Y,T,50);
 colorbar

3 Gradiente de la temperatura

Calcular [math]∇T[/math] y dibujarlo como campo vectorial. Observar gráficamente que [math]∇T[/math] es ortogonal a dichas curvas. (Nota: es muy importante que la escala de los ejes sea la misma para apreciar el que forman las curvas de nivel y el gradiente).

Nos ayudaremos de gran parte del código del programa del ejercicio anterior, por lo que mantendremos las variables [math]h[/math], [math]x[/math] e [math]y[/math], así como las matrices [math]X[/math] e [math]Y[/math] obtenidas mediante el comando [math]meshgrid ()[/math] y la matriz [math]T[/math] que representa la temperatura.

En segundo lugar, calcularemos el gradiente de la temperatura, que es un vector cuya primera componente es la derivada del campo escalar de la temperatura respecto de [math]x[/math], multiplicada por el vector [math]\vec{i}[/math] y como segunda componente la derivada del campo escalar respecto de [math]y[/math], multiplicada por el vector [math]\vec{j}[/math].

[math]∇T = \frac{\partial T}{\partial x}\vec{i} + \frac{\partial T}{\partial y}\vec{j} [/math]

Una vez calculado el gradiente, crearemos las variables [math]dx[/math] y [math]dy[/math] en función de las matrices [math]X[/math] e [math]Y[/math]. Por último, utilizaremos el comando [math]quiver ()[/math] para representar el gradiente de la temperatura.

Al tener que comprobar que las curvas de nivel del campo de temperatura y el gradiente son perprendiculares, nos ayudaremos del comando [math]hold on[/math] para poder representar las curvas de nivel y el gradiente en el mismo gráfico.

Gradiente

%pregunta3
% 3. Calcular ∇T y dibujarlo como campo vectorial. Observar gr´aficamente que ∇T es ortogonal a
% dichas curvas. (Nota: es muy importante que la escala de los ejes sea la misma para apreciar el
% ´angulo que forman las curvas de nivel y el gradiente)
clc
clear
close all
% Paso de muestreo h para las variables x e y: h = 1/10
h = 2/10;
% Región que ocupa la placa rectangular
x = [0:h:10];
y = [0:h:2];
% Creación del mallado
[X,Y]= meshgrid(x,y);
% Función temperatura: T(x,y) = (x-3)^2 + (10Y)^2
T = (X-6).^2 + (10*(Y-3/2)).^2;
contour(X,Y,T,30);
% Gradiente de T
dx = 2.*(X-6);
dy = 200*Y-300;
% Título
title('Gradiente de temperatura');
% Escribimos los nombres de los ejes
xlabel('Eje X');
ylabel('Eje Y');
% Representación de la temperatura y las curvas de nivel
hold on
quiver(x,y,dx,dy);
axis equal
colorbar

4 Gradiente del campo de desplazamientos

Consideramos ahora el campo de desplazamientos [math]\vec{u}[/math]. Calcularlo sabiendo:

• Los puntos situados en el eje [math]x = 0[/math] no sufren desplazamiento en la dirección [math]\vec{i}[/math].
• [math]∇ · \vec{u} = 0[/math].

El gradiente de un campo vectorial es una matriz 3x3 en cuyas columnas se encuentran las derivadas del campo respecto a cada componente.

[math]∇\vec{u} = \begin{pmatrix} \frac{\partial u_{1} }{\partial x} & \frac{\partial u_{1} }{\partial y} & \frac{\partial u_{1} }{\partial z} \\ \frac{\partial u_{2} }{\partial x} & \frac{\partial u_{2} }{\partial y} & \frac{\partial u_{2} }{\partial z} \\ \frac{\partial u_{3} }{\partial x} & \frac{\partial u_{3} }{\partial y} & \frac{\partial u_{3} }{\partial z} \end{pmatrix} [/math]

[math]\vec{u}(x, y)=\frac{yf(x)}{20}\vec{i}-\frac{y^2}{40}\vec{j}[/math]

Al no aparecer la componente [math]z[/math], la tercera columna de la matriz será una columna de ceros.

[math]∇\vec{u} = \begin{pmatrix} \frac{\partial u_{1} }{\partial x} & \frac{\partial u_{1} }{\partial y} & 0\\ \frac{\partial u_{2} }{\partial x} & \frac{\partial u_{2} }{\partial y} & 0 \\ \frac{\partial u_{3} }{\partial x} & \frac{\partial u_{3} }{\partial y} & 0 \end{pmatrix} [/math]

[math]\frac{\partial u_{1} }{\partial x} = \frac{y f^\prime(x)}{20} [/math] , [math]\frac{\partial u_{1} }{\partial y} = \frac{f(x)}{20} [/math]

[math]\frac{\partial u_{2} }{\partial y} = \frac{-2y}{40} = \frac{-y}{20} [/math] , [math]\frac{\partial u_{2} }{\partial x} = 0 [/math]

[math]\frac{\partial u_{3} }{\partial x} =\frac{\partial u_{3} }{\partial y} = \frac{\partial u_{3} }{\partial z} = 0 [/math]

Como [math]∇ · \vec{u} = \frac{\partial u_{1} }{\partial x} + \frac{\partial u_{2} }{\partial y} + \frac{\partial u_{3} }{\partial z} = 0 = \frac{y f^\prime(x)}{20} - \frac{y}{20} + 0 [/math]

Despejando: [math]\frac{y f^\prime(x)}{20} = \frac{y}{20} [/math]. Eliminando los denominadores: [math] y f^\prime(x) = y [/math], por lo que [math] f^\prime(x) = \frac{y}{y} = 1 [/math]

Integramos a ambos lados de la igualdad: [math] \int f^\prime(x) \, \text{d}x \,\! = f(x) = \int 1\, \text{d}x \,\! = x [/math]

Por tanto, [math] f(x) = x [/math]

[math]∇\vec{u} = \begin{pmatrix} \frac{yf^\prime(x)}{20} & \frac{f(x)}{20} & 0\\ \frac{\partial u_{2} }{\partial x} & \frac{f(x)}{20} & 0 \\ \frac{\partial u_{3} }{\partial x} & \frac{\partial u_{3} }{\partial y} & 0 \end{pmatrix} = \begin{pmatrix} \frac{y}{20} & \frac{x}{20} & 0\\ \frac{\partial u_{2} }{\partial x} & \frac{x}{20} & 0 \\ \frac{\partial u_{3} }{\partial x} & \frac{\partial u_{3} }{\partial y} & 0 \end{pmatrix} = \begin{pmatrix} \frac{y}{20} & \frac{x}{20} & 0\\ 0 & \frac{x}{20} & 0 \\ 0 & 0 & 0 \end{pmatrix}[/math]

[math]∇\vec{u} = \begin{pmatrix} \frac{y}{20} & \frac{x}{20} & 0\\ 0 & \frac{x}{20} & 0 \\ 0 & 0 & 0 \end{pmatrix}[/math]

5 Representación del campo de desplazamientos

Dibujar el campo de vectores en los puntos del mallado del sólido.

Al igual que en los ejercicios anteriores, representaremos las variables [math]h[/math], [math]x[/math] e [math]y[/math] y utilizaremos el comando [math]meshgrid ()[/math] para obtener las matrices [math]X[/math] e [math]Y[/math].

Definiremos el campo de desplazamientos mediante las variables [math]ux[/math] y [math]uy[/math] en función de las matrices anteriores, y mediante el comando [math]quiver ()[/math] representaremos las variables [math]x[/math], [math]y[/math], [math]ux[/math] y [math]uy[/math].

Por último, ajustaremos los ejes según los límites que nos han indicado y escribiremos el título del gráfico y de los ejes.

Campo de desplazamientos

clc
clear
close all
% Paso de muestreo h para las variables x e y
h=1/10; 
% Región que ocupa la placa rectangular
x=[0:h:10];          
y=[-1:h:1];
% Creación del mallado
[X,Y]=meshgrid(x,y);
% Definimos el campo de desplazamientos
ux=X.*Y./20;
uy=-(Y.^2)./40;
% Dibujamos el campo de vectores en los puntos del mallado del sóldio
quiver(x,y,ux,uy)
% Establecemos el límite de los ejes
axis([-0.5,10.5,-1.5,1.5]);
axis equal
% Escribimos el título del gráfico
title('Campo de desplazamientos');
% Escribimos los nombres de los ejes
xlabel('Eje X');
ylabel('Eje Y');
view(2);

6 Sólido antes y después del desplazamiento

Dibujar el sólido antes y después del desplazamiento dado por el campo de vectores [math]u[/math]. Dibujar ambos en la misma figura usando el comando subplot.

Situación antes del desplazamiento

Situación después del desplazamiento

Comparación entre ambas situaciones

%Limpieza de programas anteriores
clc
clear
close all

%Definición de regiones
h=1/10; 
x=0:h:10;          
y=-1:h:1;

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

%Vector desplazamiento
Ux=Mx.*My./20;
Uy=-(My.^2)./40;

figure

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

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

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

7 Divergencia del campo vectorial sobre la placa

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?

En el apartado 4 se indica que la divergencia del desplazamiento es nula, por lo que ésta será nula en todos los puntos. En la representación la placa tendrá un color uniforme y no habrá cambio de volumen debido al desplazamiento en ningún punto.

Divergencia del desplazamiento

%Limpieza de programas anteriores
clc
clear
close all

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

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

%Según lo indicado en el apartado 4, la divergencia es nula
div=Mx.*0+My.*0;

%Se puede representar la divergencia, pero la placa tendrá un color uniforme porque es nula
surf(Mx,My,div)
shading flat
axis([-0.5,10.5,-1.5,1.5]);
axis equal
view(2);
title('Divergencia del vector desplazamiento (nula)');
xlabel('Eje X');
ylabel('Eje Y');
colorbar

8 Rotacional del sólido

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

8.1 Cálculo del rotacional

[math]∇ × \vec{u}= \begin{vmatrix}\vec{i} & \vec{j}& \vec{k}\\ \frac{\partial }{\partial x} & \frac{\partial }{\partial y}& \frac{\partial }{\partial z}\\ u_{1}& u_{2}& u_{3} \end{vmatrix} = \begin{vmatrix} \vec{i} & \vec{j} & \vec{k}\\ \frac{\partial }{\partial x} & \frac{\partial }{\partial y} &\frac{\partial }{\partial z} \\ \frac{yx}{20} & -\frac{y^2}{40} & 0\end{vmatrix} = -\frac{x}{20}\vec{k} [/math].

8.2 Cálculo del módulo del rotacional

[math]|∇ × \vec{u}|=\frac{x}{20}[/math]

8.3 Representación del módulo del rotacional

Mientras que el rotacional se representa como un campo escalar, el módulo del rotacional se trata de un campo escalar. Como puede verse en la imagen adjunta, los puntos incrementan su módulo del rotacional según avanzan positivamente en el eje x. Así, los puntos de la placa con mayor módulo del rotacional serán los situados en [math]x = 10[/math].

Módulo del rotacional (campo escalar)

%Limpieza de programas anteriores
clc
clear
close all

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

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

%Cálculo del módulo del rotacional
rot=Mx./20;

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

9 Tensor de tensiones

Definamos [math]Ԑ(\vec{u}) = (∇\vec{u} + ∇\vec{u}t)/2[/math], la parte simétrica del tensor gradiente de [math]\vec{u}[/math] conocido como tensor de deformaciones. En un medio elástico lineal, isótropo y homogéneo los desplazamientos permiten escribir el tensor de tensiones [math]σ_{ij}[/math] a través de la fórmula [math]σ = λ∇ · \vec{u} 1 + 2 μԐ[/math], donde λ 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], dibujar las tensiones normales en la dirección que marca el eje [math]\vec{i}[/math], es decir [math]\vec{i}·σ·\vec{i}[/math], las tensiones normales en la dirección que marca el eje [math]\vec{j}[/math], es decir [math]\vec{j}·σ·\vec{j}[/math] y las correspondientes al eje [math]\vec{k}[/math], es decir [math]\vec{k}·σ·\vec{k}[/math].

Para el cálculo de la parte simétrica del tensor gradiente del campo vectorial [math]U[/math], utilizaremos la fórmula correspondiente, con la cual construiremos una matriz cuadrada como se vio anteriormente. La tercera columna es toda 0 por no tener el campo componente en [math]\vec{k}[/math]. Ahora, sumando este resultado a su traspuesto y dividiendo todo entre dos, obtenemos el tensor deformación, el cual es una matriz con las siguientes filas:

[math]Ԑ(\vec{u})=[/math] \begin{pmatrix} \frac{y}{20} & \frac{x}{40} & 0 \\ \frac{x}{40}& \frac{-y}{20} & 0 \\ 0 & 0 & 0 \end{pmatrix}

Al ser la divergencia del campo igual a 0 y el parámetro necesario igual a 1, al aplicar la fórmula proporcionada por el enunciado para el cálculo del tensor de tensiones, obtenemos un tensor de tensiones que no es sino el resultado de multiplicar por dos el anterior.

Aplicando las fórmulas correspondientes y proporcionadas por el enunciado, se obtienen unas tensiones normales en las direcciones de los ejes [math]\vec{i}[/math], [math]\vec{j}[/math] y [math]\vec{k}[/math], respectivamente de:

[math]\vec{i}·σ·\vec{i}=\frac {y}{10}[/math]

[math]\vec{j}·σ·\vec{j}=\frac {-y}{10}[/math]

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

Representación de las tensiones normales en la dirección [math]\vec{i}[/math]

Representación de las tensiones normales en la dirección [math]\vec{j}[/math]

Representación de las tensiones normales en la dirección [math]\vec{k}[/math]

%limpieza de todo lo anterior
clc
clear
close all
%definición de las variables
h=1/10; 
x=(0:h:10);          
y=(-1:h:1);

[X,Y]=meshgrid(x,y);
tensioni=Y./10;
tensionj=-Y./10;
tensionk=X-X;
%tensiones normales en la dirección normal del eje i
figure

surf(X,Y,tensioni)
shading flat
axis equal
colorbar()
axis([-0.5,10.5,-1.5,1.5]);
title('Tensión normal en dirección i');
xlabel('Eje X');
ylabel('Eje Y');
view(2);
%tensiones normales en la dirección normal del eje j
figure

surf(X,Y,tensionj)
shading flat
axis equal
colorbar()
axis([-0.5,10.5,-1.5,1.5]);
title('Tensión normal en dirección j');
xlabel('Eje X');
ylabel('Eje Y');
view(2);
%tensiones normales en la dirección normal del eje k (la variable
%tensionesk se definió como x-x para hacer ver que siempre vale 0)
figure

surf(X,Y,tensionk)
shading flat
axis equal
colorbar()
axis([-0.5,10.5,-1.5,1.5]);
title('Tensión normal en dirección k');
xlabel('Eje X');
ylabel('Eje Y');
view(2);

10 Tensiones tangenciales

Calcular las tensiones tangenciales respecto al plano ortogonal a [math]\vec{i}[/math], es decir [math]|σ·\vec{i} − (\vec{i}·σ·\vec{i})\vec{i}|[/math].

La fórmula proporcionada por el enunciado permite calcular las componentes de la tensión que experimenta el sólido en cada punto pertenecientes a un plano perpendicular a [math]\vec{i}[/math]. El cálculo resulta muy sencillo al cancelarse entre sí dos de los tres términos. Recordemos que para calcular el vector resultado de aplicar un tensor a otro vector se obtiene simplemente haciendo un producto matricial de la matriz de componentes del tensor en cuestión y la matriz columna de las componentes del vector en la misma base ortonormal.

[math]|σ·\vec{i} − (\vec{i}·σ·\vec{i})\vec{i}|=|\frac {y}{10}\vec{i}+\frac {x}{20}\vec{j}-\frac {y}{10}\vec{i})|=\frac {x}{20}[/math]

11 Tensión de Von Mises

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

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

donde [math]σ_{1}[/math], [math]σ_{2}[/math] y [math]σ_{3}[/math] son los autovalores de [math]σ[/math] (también conocidos como tensiones principales). Se trata de una magnitud escalar que se suele usar como indicador para saber cuando un material inicia un comportamiento plástico (y no elástico puro).

Sus valores aumentan una cantidad considerable al alejarse del eje Y y del eje [math]X[/math]. Dado esto sus valores máximos son los que se dan en las esquinas de la derecha de la placa. Ambas tienen el mismo valor de [math]0.88318[/math]. A continuacion se muestra su calculo en Octave/MatLab.

Tension de Von Mises

%Limpieza de programas anteriores
clear
clc
close all

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

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

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

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

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

12 Campo de fuerzas que actúa sobre la placa

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

[math]\vec{F}=-∇· σ[/math]

donde [math]∇ · σ[/math] es el campo vectorial que se obtiene al hacer la divergencia de los vectores cuyas componentes son las filas de la matriz [math]σ[/math].

12.1 Cálculo del campo de fuerzas

[math] \vec F = -\frac{∂σji}{∂xj}\vec e = \frac{\frac{∂y}{10}}{∂x}\vec i+ \frac{\frac{∂x}{20}}{∂y}\vec i+ \frac{∂0}{∂z}\vec i+ \frac{\frac{∂x}{20}}{∂x}\vec j+ \frac{\frac{∂-y}{10}}{∂y}\vec j+ \frac{∂0}{∂z}\vec j + \frac{∂0}{∂x}\vec k + \frac{∂0}{∂y}\vec k + \frac{∂0}{∂z}\vec k = - \frac{1}{20}\vec j + \frac{1}{10}\vec j[/math]

[math]\vec F = \frac{1}{20}\vec j [/math]

En este caso el campo de fuerzas es uniforme en toda la placa, tiene un valor de [math]\frac {1}{20}[/math] en la direccion de [math]Y[/math].

12.2 Representación del campo de fuerzas

Campo de fuerzas

%Limpieza de programas anteriores
clear
clc
close all

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

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

%Valor del campo de fuerzas para todos los puntos
F=1/20+Y.*0;

%Graficamos
quiver(x,y,0.*X,F)
axis([-0.5,10.5,-1.5,1.5]);
axis equal
title('Campo de Fuerzas');
xlabel('Eje X');
ylabel('Eje Y');
view(2);