MateWiki - Contribuciones del usuario [es]

Visualización de campos escalares y vectoriales en elasticidad. Grupo 10-A

2016-12-04T16:41:05Z

Alejandro de Benito Andrés:

{{ TrabajoED | Visualizaciones de campos escalares y vectoriales. Grupo 10-A | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC16/17|2016-17]] | Alejandro de Benito, Lucía Rico , Adrián Piñeiro, Daniel Apellániz, Guillermo Alfaya. }}

Dada una superficie plana en forma de anillo circular, concretamente la cuarta parte situada en el primer cuadrante, con centro en el origen de radio exterior 3 e interior 1. Se han utilizado como ejes el cuadrado [-1,4]x[-1,4] para una mejor visualización de la placa, a la que se aplicará una temperatura T(x,y)=(y+2)·x2, en la que se toman como variables independientes "x, y" y el campo de desplazamientos u(ρ,θ)=θf(ρ)gρ, conclusión de un agente exterior desconocido. En apartados posteriores veremos como afectan la temperatura o el campo de desplazamientos a la placa en cuestión con las tensiones ejercidas. Por comodidad han sido elegidas las coordenadas cilíndricas.

== Mallado ==
Representaremos una malla de los puntos que componen la placa sólida con forma de corona semicircular.
El paso de muestreo dado es 0,1 tanto para la variable X como para la variable Y.
[[Archivo:Malladogrupo10A.png|thumb|400px|derecha]]

{{matlab|codigo=
%Definimos las variables con el paso de muestreo "h"
h=0.1
rho=[1:h:3];
theta=[0:h:pi/2];

%Creamos el cojunto de puntos que forman la malla y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis([-1,4,-1,4]);

view(2);
title('Mallado')
}}

== Temperatura ==

En la corona la temperatura varía según la función dada, "T(x,y)=(y+2)·x2". Para tener una idea más realista de la temperatura la hemos dibujado tanto como en 2D como curvas de nivel. El punto que ha resultado tener mayor temperatura T(x,y)=24.1902 es el punto x=1.1683,y=2.7632.

[[File:temperaturaGrupo10A2.png|thumb|800px|derecha]]

{{matlab|codigo=
%Definimos las variables con el paso de muestreo "h"
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H); %Así se ajusta mejor a theta=pi/2

%Creamos el cojunto de puntos que forman la malla y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis equal;

view(2);

%CAMPO ESCALAR DE TEMPERATURAS
Temperatura=(Y+2).*X.^2
Maximo=max(max(Temperatura))

%Líneas de programa usadas para obtener posición del máximo
for i=1:16
for j=1:21
Temp=(Y(i,j)+2).*X(i,j).^2;
if Temp==Maximo
disp(X(i,j))
disp(Y(i,j))
end
end
end

%Representación de la malla en 2D en la subventana primera
subplot(1,2,1)

contour(X,Y,Temperatura,20);

view(2)
axis([-1,4,-1,4]);
title('Curvas de nivel de la Temperatura')

%Representación de la malla en 3D en la subventana segunda
subplot(1,2,2);

surf(X,Y,Temperatura);

view(2)
axis([-1,4,-1,4]);
title('Temperatura en 2D')
}}

== Gradiente de la temperatura ==

Para conocer la variación de temperatura por la placa usaremos el gradiente del campo de temperaturas(∇T), el módulo de los vectores gradiente representa la rapidez con la que varía la temperatura, es decir apunta hacia la dirección en la que la temperatura tiene un aumento más notable, donde la derivada direccional es máxima. Se aprecia en la figura que el vector gradiente en perpendicular a la superficie de sólido.

[[Archivo:gradienteDeLaTemperaturaGrupo10A1.png|thumb|800px|Izquierda: Gráfica de la temperatura en 2D y vectores gradientes. Derecha: Vectores del gradiente|derecha]]
{{matlab|codigo=
%Definimos las variables con el paso de muestreo "h"
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el cojunto de puntos que forman la malla y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis([-1,4,-1,4]);
axis square

view(2);

%GRADIENTE
GradX=2*X.*Y+4*X;
GradY=X.^2;

%CAMPO ESCALAR DE TEMPERATURAS
Temperatura=(Y+2).*X.^2;

%Creación de la subventana primera: gradiente en el campo de temperaturas

subplot(1,2,1)
hold on
surf(X,Y,Temperatura);
axis([-1,4,-1,4]);
axis square

%Dibujo del gradiente
quiver3(X,Y,Temperatura,GradX,GradY,0*Temperatura);
axis([-1,4,-1,4]);
axis square

hold off

%Creación de la segunda subventana: campo vectorial del gradiente

subplot(1,2,2)
quiver(X,Y,GradX,GradY);
axis([-1,4,-1,4]);
axis equal
colorbar
}}

== Campo de desplazamientos ==

En este apartado procederemos a dibujar el campo de desplazamientos "u(ρ,θ)=θf(ρ)gρ". Las condiciones que conocemos son: los puntos de ρ=1 no sufren desplazamiento y div(u)=θ(2-1/ρ)/5, lo que nos da la ecuación diferencial f'(ρ)+(1/ρ)*f(ρ)=(2-1/ρ)/5 con la condición de contorno f(1)=0. Resultando f(ρ)=(ρ-1)/5. Lo metemos en el campo de desplazamiento dando u(ρ,θ)=θ*((ρ-1)/5)gρ .Cogiendo el campo en función de los dos componente i,j. Definiremos pues FX y FY. Sabiendo que gρ= cosθi+senθj.

[[Archivo:CampoDesplazamientosGrupo10A1.png|marco|derecha]]

{{matlab|codigo=
%Definimos las variables con el paso de muestreo "h"
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el cojunto de puntos que forman la malla y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis([-1,4,-1,4]);

view(2);

%CAMPO VECTORIAL DE DESPLAZAMIENTOS
%Componente x(X,Y) del campo vectorial
FX=THETA.*((RHO-1)/5).*(cos(THETA));

%Componente y(X,Y) del campo vectorial
FY=THETA.*((RHO-1)/5).*(sin(THETA));

%Representamos el campo vectorial de desplazamientos
quiver(X,Y,FX,FY)

axis([-1,4,-1,4])
view(2)

xlabel('Eje x')
ylabel('Eje y')
title('Campo de desplazamientos sobre la placa')
}}

== Movimiento del sólido ==

En este caso se muestran dos gráficas una con la placa en reposo y otra después de aplicarle el campo de desplazamientos a la superficie, de esta forma es perfectamente apreciable la variación que sufre. El desplazamiento viene dado por dos componentes, el vector desplazamiento y los puntos de la placa respectivos a cada coordenada. Cuando ambas se suman, dan lugar a la figura final.

[[Archivo:movimientoDelSolidoGrupo10A1.png|thumb|800px|derecha]]

{{matlab|codigo=
%Inicializamos las variables
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis([-1,4,-1,4]);

view(2);

%CAMPO VECTORIAL DE DESPLAZAMIENTOS
%Componente x(X,Y) del campo vectorial
FX=THETA.*((RHO-1)/5).*(cos(THETA));

%Componente y(X,Y) del campo vectorial
FY=THETA.*((RHO-1)/5).*(sin(THETA));

%Dibujo de la placa sin movimiento
subplot(1,2,1)

surf(X,Y,0*X);

axis([-1,4,-1,4]);
view(2);
title('Reposo')

%Dibujo de la placa en movimiento
%Movimiento de la placa en el eje X
MovimientoX=X+FX;

%Movimiento de la placa en el eje Y
MovimientoY=Y+FY;

subplot(1,2,2)

surf(MovimientoX,MovimientoY,Y.*0)

axis([-1,4,-1,4]);
view(2);
title('Aplicado el campo de movimientos')
}}

== Divergencia ==

Dada la divergencia usada para definir el campo "∇·(u)=θ(2-1/ρ)/5", procederemos a evaluar sus máximos, mínimos y nulos con el gradiente<big>∇(∇·u)=θ/(5*ρ²)gρ+(2-1/ρ)/5*(1/ρ²)gθ+0gz</big>.
Al tener un dominio de ρ contenido en (1,3) y de θ en (0,pi/2), la derivada direccional gρ y gθ son positivas, también podemos sacar en claro que cuanto mayor sean ρ y θ mayor será su valor, llegando a su punto máximo en (3,pi/2) con: ∇·u(ρ,θ)=0.5236. A su vez el mínimo se encontrará donde los θ sean menores en nuestro dominio θ=0 dando lugar a: ∇·(u)=0 por estar ∇·(u) en función de θ.

[[Archivo:DivergenciaGrupo10A.png|marco|derecha]]

{{matlab|codigo=
%Inicializamos las variables
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis([-1,4,-1,4]);

view(2);

%DIVERGENCIA
DIVERGENCIA=THETA.*((2-1./RHO)./5);

MaxDivergencia=max(max(DIVERGENCIA))

%Representación de la divergencia
surf(X,Y,DIVERGENCIA)

view(2)
colorbar
title('Divergencia')
}}

== Rotacional ==

Se ha calculado el rotacional del campo de desplazamientos u en todos los puntos del sólido y lo representamos. Hemos obtenido que el rotacional es: Rot [u(rho,theta)] = (rho-1)/(5*rho)gz
Observando la gráfica se obtiene cuáles son los puntos que sufren un mayor rotacional: Los puntos con mayor rotacional son aquellos con theta=PI/4 y '''su valor es 0.1333'''

El rotacional del campo de desplazamientos se ha calculado y aplicado a todos los puntos de la placa y han sido representados con colores. El ca

[[Archivo:RotacionalGrupo10A.png|marco|derecha]]

{{matlab|codigo=
%Inicializamos las variables
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis([-1,4,-1,4]);

view(2);

%ROTACIONAL
ROTACIONAL=(RHO-1)./(5.*RHO);

%Representación del rotacional
surf(X,Y,ROTACIONAL)

colorbar;
view(2)
title('Rotacional')

%Puntos con mayor rotacional
Maximo_Rotacional=max(max(ROTACIONAL))

%El máximo del rotacional corresponde a puntos con THETA = (PI/4)
title('Rotacional')

%Puntos con mayor rotacional
Maximo_Rotacional=max(max(ROTACIONAL))

}}

== Tensor de tensiones ==
Si consideramos que estamos en un medio elástico lineal, isótropo y homogéneo; en mecánica de medios continuos podemos definir el tensor de tensiones "σ"como:

[[Archivo:TensorDeTensionesGrupo10A.png|marco|centro]]

siendo λ y μ los coeficientes de Lamé, que dependen de las propiedades elásticas de cada material, considerando en nuestro caso que λ y μ son iguales a 1 y que ε es la parte simétrica del vector gradiente u.

Para hacer este apartado se requiere utilizar bases generales, en lo cual hemos llegado a la conclusión de que las componentes de ∇ū son:

[[Archivo:GradienteDeUGrupo10A.png|marco|centro]]

Siendo ε =(∇u + ∇ut)/2, con todo ello, calculamos la matriz σ:

[[Archivo:MatrizDeTensionesGrupo10A.png|marco|centro]]

Se van a representar las tensiones en las tres direcciones del espacio. En el caso de estar utilzando coordenadas cilíndricas, se representan las tensiones en las direcciones de la base cilíndrica ḡρ, ḡθ, ḡz. '''Para hallar las tensiones en estas direcciones debemos multiplicar escalarmente la matriz σ por cada uno de los vectores ḡi por ambos lados: Tensión normal a la dirección ḡi = ḡi*σ*ḡi'''

[[Archivo:DireccionGrupo10A.png|marco|centro]]

{{matlab|codigo=
%Inicializamos las variables
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%En 2D

%Dibujo de las tensiones respecto a g sub rho
subplot(2,3,1);

%Tomamos el valor (1,1) de la matriz sigma y lo dibujamos
r=2.*THETA./5+THETA.*(2-1./RHO)./5;

surf(X,Y,r)

axis equal
view(2)
title('Tensiones respecto a gRho')
colorbar

%Dibujo de las tensiones respecto a g sub tetha
subplot(2,3,2);

%Tomamos el valor (2,2) de la matriz sigma y lo dibujamos
t=2.*THETA.*(RHO-1)./(5.*RHO.^5)+THETA.*(2-1./RHO)./(5.*RHO.^2);

surf(X,Y,t)

axis equal
title('Tensiones respecto a gTheta')
view(2)

colorbar

%Dibujo de las tensiones respecto a g sub z
subplot(2,3,3)

%Tomamos el valor (3,3) de la matriz sigma y lo dibujamos
z=THETA.*(2-1./RHO)./5;

surf(X,Y,z)

axis equal
view(2)
title('Tensiones respecto a gz')
colorbar

%Cálculo en 3D:

%Dibujo de las tensiones respecto a g sub rho 3D
subplot(2,3,4);

%Tomamos el valor (1,1) de la matriz sigma y lo dibujamos

surf(X,Y,r)

axis equal
title('Tensiones respecto a gRho 3D')

%Dibujo de las tensiones respecto a g sub tetha 3D
subplot(2,3,5);

%Tomamos el valor (2,2) de la matriz sigma y lo dibujamos

surf(X,Y,t)

axis equal
title('Tensiones respecto a gTheta 3D')

%Dibujo de las tensiones respecto a g sub z 3D
subplot(2,3,6)

%Tomamos el valor (3,3) de la matriz sigma y lo dibujamos

surf(X,Y,z)

axis equal
title('Tensiones respecto a gz 3D')
}}

[[Archivo:tensionesNormalesgrupo10A.png|centro|thumb|1300px]]

== Tensiones tangenciales ==

'''Ambas tensiones tangenciales son iguales, pues la matriz σ es simétrica.
Las tensiones del apartado nos dan: σ(2,1)= σ(2,1)= 1/ρ*(1/5*sen(2θ)-1)'''

== Tensiones tangenciales respecto al plano ortogonal a ḡρ ==

[[Archivo:Tensionestangencialesgrupo10A.png||thumb|390px|derecha]]

{{matlab|codigo=
%Inicializamos las variables
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos las tensiones tangentes al plano g sub rho
Tang=(RHO-1)./(5.*RHO.^2);

%Visualización en 3D
subplot(2,1,1)

surf(X,Y,Tang)

colorbar
title('Tensiones tangenciales respecto al plano ortogonal a gRho en 3D')

%Visualización en 2D
subplot(2,1,2)

surf(X,Y,Tang)

colorbar
view(2)
title('Tensiones tangenciales respecto al plano ortogonal a gRho en 2D')
}}

== Tensión de Von Mises ==

La tensión de Von Mises se define mediante la expresión:

[[Archivo:FormulaVonMises.png|marco|centro]]

La tensión de Von Mises es una magnitud escalar que se suele utilizar como indicador para saber cuando un material inicia un comportamiento plástico y no elástico puro. En nuestro caso, el valor máximo de la Tensión de Von Mises es igual a 0.6283.

[[Archivo:tensionVonMisesgrupo10A.png|marco|derecha]]

{{matlab|codigo=
%Inicializamos las variables
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos la parametrización
[RHO,THETA]=meshgrid(rho,theta);

%Creamos la matriz sigma
Sigma=zeros(3,3);

%Creamos la matriz de Von Mises
MatrizVonMises=zeros(16,21);

%Obtenemos los valores de la tensión de Von Mises a partir de los autovalores
for i=1:16*21

r=RHO(i)';

t=THETA(i)';

Sigma(1,1)=2*t/5+t*(2-1/r)/5;

Sigma(1,2)=(r-1)/(5*r^2);

Sigma(2,1)= (r-1)/(5*r^2);

Sigma(2,2)=2*t*(r-1)/(5*r^3)+t*(2-1/r)/(5);
Sigma(3,3)=t*(2-1/r)/5;

%Obtención de los autovalores de la matriz sigma
[v,d]=eig(Sigma);

%Fórmula de Von Mises
MatrizVonMises(i)=sqrt(((d(1,1)-d(2,2))^2+(d(2,2)-d(3,3))^2+(d(3,3)-d(1,1))^2)/2);

end

%Cambiamos a coordenadas cartesianas
X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos la tensión de Von Mises
surf(X,Y,MatrizVonMises)

colorbar
view(2)
title('Tensión de Von Mises')

%Valor Máximo
[i,j]=find(MatrizVonMises==max(max((MatrizVonMises))))
z=max(max(MatrizVonMises))

%El valor máximo de la Tensión de Von Mises es igual a 0.6283
}}

Visualización de campos escalares y vectoriales en elasticidad. Grupo 10-A

2016-12-04T16:31:08Z

Alejandro de Benito Andrés:

{{ TrabajoED | Visualizaciones de campos escalares y vectoriales. Grupo 10-A | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC16/17|2016-17]] | Alejandro de Benito, Lucía Rico , Adrián Piñeiro, Daniel Apellániz, Guillermo Alfaya. }}

Dada una superficie plana en forma de anillo circular, concretamente la cuarta parte situada en el primer cuadrante, con centro en el origen de radio exterior 3 e interior 1. Se han utilizado como ejes el cuadrado [-1,4]x[-1,4] para una mejor visualización de la placa. Por comodidad han sido elegidas las coordenadas cilíndricas, a la que se aplicará una temperatura T(x,y)=(y+2)·x2, en la que se toman como variables independientes "x, y" y el campo de desplazamientos u(ρ,θ)=θf(ρ)gρ, conclusión de un agente exterior desconocido. En apartados posteriores veremos como afectan la temperatura o el campo de desplazamientos a la placa en cuestión con las tensiones ejercidas.

== Mallado ==
Representaremos una malla de los puntos que componen la placa sólida con forma de corona semicircular.
El paso de muestreo dado es 0,1 tanto para la variable X como para la variable Y.
[[Archivo:Malladogrupo10A.png|thumb|400px|derecha]]

{{matlab|codigo=
%Definimos las variables con el paso de muestreo "h"
h=0.1
rho=[1:h:3];
theta=[0:h:pi/2];

%Creamos el cojunto de puntos que forman la malla y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis([-1,4,-1,4]);

view(2);
title('Mallado')
}}

== Temperatura ==

En la corona la temperatura varía según la función dada, "T(x,y)=(y+2)·x2". Para tener una idea más realista de la temperatura la hemos dibujado tanto como en 2D como curvas de nivel. El punto que ha resultado tener mayor temperatura T(x,y)=24.1902 es el punto x=1.1683,y=2.7632.

[[File:temperaturaGrupo10A2.png|thumb|800px|derecha]]

{{matlab|codigo=
%Definimos las variables con el paso de muestreo "h"
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H); %Así se ajusta mejor a theta=pi/2

%Creamos el cojunto de puntos que forman la malla y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis equal;

view(2);

%CAMPO ESCALAR DE TEMPERATURAS
Temperatura=(Y+2).*X.^2
Maximo=max(max(Temperatura))

%Líneas de programa usadas para obtener posición del máximo
for i=1:16
for j=1:21
Temp=(Y(i,j)+2).*X(i,j).^2;
if Temp==Maximo
disp(X(i,j))
disp(Y(i,j))
end
end
end

%Representación de la malla en 2D en la subventana primera
subplot(1,2,1)

contour(X,Y,Temperatura,20);

view(2)
axis([-1,4,-1,4]);
title('Curvas de nivel de la Temperatura')

%Representación de la malla en 3D en la subventana segunda
subplot(1,2,2);

surf(X,Y,Temperatura);

view(2)
axis([-1,4,-1,4]);
title('Temperatura en 2D')
}}

== Gradiente de la temperatura ==

Para conocer la variación de temperatura por la placa usaremos el gradiente del campo de temperaturas(∇T), el módulo de los vectores gradiente representa la rapidez con la que varía la temperatura, es decir apunta hacia la dirección en la que la temperatura tiene un aumento más notable, donde la derivada direccional es máxima. Se aprecia en la figura que el vector gradiente en perpendicular a la superficie de sólido.

[[Archivo:gradienteDeLaTemperaturaGrupo10A1.png|thumb|800px|Izquierda: Gráfica de la temperatura en 2D y vectores gradientes. Derecha: Vectores del gradiente|derecha]]
{{matlab|codigo=
%Definimos las variables con el paso de muestreo "h"
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el cojunto de puntos que forman la malla y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis([-1,4,-1,4]);
axis square

view(2);

%GRADIENTE
GradX=2*X.*Y+4*X;
GradY=X.^2;

%CAMPO ESCALAR DE TEMPERATURAS
Temperatura=(Y+2).*X.^2;

%Creación de la subventana primera: gradiente en el campo de temperaturas

subplot(1,2,1)
hold on
surf(X,Y,Temperatura);
axis([-1,4,-1,4]);
axis square

%Dibujo del gradiente
quiver3(X,Y,Temperatura,GradX,GradY,0*Temperatura);
axis([-1,4,-1,4]);
axis square

hold off

%Creación de la segunda subventana: campo vectorial del gradiente

subplot(1,2,2)
quiver(X,Y,GradX,GradY);
axis([-1,4,-1,4]);
axis equal
colorbar
}}

== Campo de desplazamientos ==

En este apartado procederemos a dibujar el campo de desplazamientos "u(ρ,θ)=θf(ρ)gρ". Las condiciones que conocemos son: los puntos de ρ=1 no sufren desplazamiento y div(u)=θ(2-1/ρ)/5, lo que nos da la ecuación diferencial f'(ρ)+(1/ρ)*f(ρ)=(2-1/ρ)/5 con la condición de contorno f(1)=0. Resultando f(ρ)=(ρ-1)/5. Lo metemos en el campo de desplazamiento dando u(ρ,θ)=θ*((ρ-1)/5)gρ .Cogiendo el campo en función de los dos componente i,j. Definiremos pues FX y FY. Sabiendo que gρ= cosθi+senθj.

[[Archivo:CampoDesplazamientosGrupo10A1.png|marco|derecha]]

{{matlab|codigo=
%Definimos las variables con el paso de muestreo "h"
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el cojunto de puntos que forman la malla y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis([-1,4,-1,4]);

view(2);

%CAMPO VECTORIAL DE DESPLAZAMIENTOS
%Componente x(X,Y) del campo vectorial
FX=THETA.*((RHO-1)/5).*(cos(THETA));

%Componente y(X,Y) del campo vectorial
FY=THETA.*((RHO-1)/5).*(sin(THETA));

%Representamos el campo vectorial de desplazamientos
quiver(X,Y,FX,FY)

axis([-1,4,-1,4])
view(2)

xlabel('Eje x')
ylabel('Eje y')
title('Campo de desplazamientos sobre la placa')
}}

== Movimiento del sólido ==

En este caso se muestran dos gráficas una con la placa en reposo y otra después de aplicarle el campo de desplazamientos a la superficie, de esta forma es perfectamente apreciable la variación que sufre. El desplazamiento viene dado por dos componentes, el vector desplazamiento y los puntos de la placa respectivos a cada coordenada. Cuando ambas se suman, dan lugar a la figura final.

[[Archivo:movimientoDelSolidoGrupo10A1.png|thumb|800px|derecha]]

{{matlab|codigo=
%Inicializamos las variables
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis([-1,4,-1,4]);

view(2);

%CAMPO VECTORIAL DE DESPLAZAMIENTOS
%Componente x(X,Y) del campo vectorial
FX=THETA.*((RHO-1)/5).*(cos(THETA));

%Componente y(X,Y) del campo vectorial
FY=THETA.*((RHO-1)/5).*(sin(THETA));

%Dibujo de la placa sin movimiento
subplot(1,2,1)

surf(X,Y,0*X);

axis([-1,4,-1,4]);
view(2);
title('Reposo')

%Dibujo de la placa en movimiento
%Movimiento de la placa en el eje X
MovimientoX=X+FX;

%Movimiento de la placa en el eje Y
MovimientoY=Y+FY;

subplot(1,2,2)

surf(MovimientoX,MovimientoY,Y.*0)

axis([-1,4,-1,4]);
view(2);
title('Aplicado el campo de movimientos')
}}

== Divergencia ==

Dada la divergencia usada para definir el campo "∇·(u)=θ(2-1/ρ)/5", procederemos a evaluar sus máximos, mínimos y nulos con el gradiente<big>∇(∇·u)=θ/(5*ρ²)gρ+(2-1/ρ)/5*(1/ρ²)gθ+0gz</big>.
Al tener un dominio de ρ contenido en (1,3) y de θ en (0,pi/2), la derivada direccional gρ y gθ son positivas, también podemos sacar en claro que cuanto mayor sean ρ y θ mayor será su valor, llegando a su punto máximo en (3,pi/2) con: ∇·u(ρ,θ)=0.5236. A su vez el mínimo se encontrará donde los θ sean menores en nuestro dominio θ=0 dando lugar a: ∇·(u)=0 por estar ∇·(u) en función de θ.

[[Archivo:DivergenciaGrupo10A.png|marco|derecha]]

{{matlab|codigo=
%Inicializamos las variables
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis([-1,4,-1,4]);

view(2);

%DIVERGENCIA
DIVERGENCIA=THETA.*((2-1./RHO)./5);

MaxDivergencia=max(max(DIVERGENCIA))

%Representación de la divergencia
surf(X,Y,DIVERGENCIA)

view(2)
colorbar
title('Divergencia')
}}

== Rotacional ==

Se ha calculado el rotacional del campo de desplazamientos u en todos los puntos del sólido y lo representamos. Hemos obtenido que el rotacional es: Rot [u(rho,theta)] = (rho-1)/(5*rho)gz
Observando la gráfica se obtiene cuáles son los puntos que sufren un mayor rotacional: Los puntos con mayor rotacional son aquellos con theta=PI/4 y '''su valor es 0.1333'''

El rotacional del campo de desplazamientos se ha calculado y aplicado a todos los puntos de la placa y han sido representados con colores. El ca

[[Archivo:RotacionalGrupo10A.png|marco|derecha]]

{{matlab|codigo=
%Inicializamos las variables
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis([-1,4,-1,4]);

view(2);

%ROTACIONAL
ROTACIONAL=(RHO-1)./(5.*RHO);

%Representación del rotacional
surf(X,Y,ROTACIONAL)

colorbar;
view(2)
title('Rotacional')

%Puntos con mayor rotacional
Maximo_Rotacional=max(max(ROTACIONAL))

%El máximo del rotacional corresponde a puntos con THETA = (PI/4)
title('Rotacional')

%Puntos con mayor rotacional
Maximo_Rotacional=max(max(ROTACIONAL))

}}

== Tensor de tensiones ==
Si consideramos que estamos en un medio elástico lineal, isótropo y homogéneo; en mecánica de medios continuos podemos definir el tensor de tensiones "σ"como:

[[Archivo:TensorDeTensionesGrupo10A.png|marco|centro]]

siendo λ y μ los coeficientes de Lamé, que dependen de las propiedades elásticas de cada material, considerando en nuestro caso que λ y μ son iguales a 1 y que ε es la parte simétrica del vector gradiente u.

Para hacer este apartado se requiere utilizar bases generales, en lo cual hemos llegado a la conclusión de que las componentes de ∇ū son:

[[Archivo:GradienteDeUGrupo10A.png|marco|centro]]

Siendo ε =(∇u + ∇ut)/2, con todo ello, calculamos la matriz σ:

[[Archivo:MatrizDeTensionesGrupo10A.png|marco|centro]]

Se van a representar las tensiones en las tres direcciones del espacio. En el caso de estar utilzando coordenadas cilíndricas, se representan las tensiones en las direcciones de la base cilíndrica ḡρ, ḡθ, ḡz. '''Para hallar las tensiones en estas direcciones debemos multiplicar escalarmente la matriz σ por cada uno de los vectores ḡi por ambos lados: Tensión normal a la dirección ḡi = ḡi*σ*ḡi'''

[[Archivo:DireccionGrupo10A.png|marco|centro]]

{{matlab|codigo=
%Inicializamos las variables
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%En 2D

%Dibujo de las tensiones respecto a g sub rho
subplot(2,3,1);

%Tomamos el valor (1,1) de la matriz sigma y lo dibujamos
r=2.*THETA./5+THETA.*(2-1./RHO)./5;

surf(X,Y,r)

axis equal
view(2)
title('Tensiones respecto a gRho')
colorbar

%Dibujo de las tensiones respecto a g sub tetha
subplot(2,3,2);

%Tomamos el valor (2,2) de la matriz sigma y lo dibujamos
t=2.*THETA.*(RHO-1)./(5.*RHO.^5)+THETA.*(2-1./RHO)./(5.*RHO.^2);

surf(X,Y,t)

axis equal
title('Tensiones respecto a gTheta')
view(2)

colorbar

%Dibujo de las tensiones respecto a g sub z
subplot(2,3,3)

%Tomamos el valor (3,3) de la matriz sigma y lo dibujamos
z=THETA.*(2-1./RHO)./5;

surf(X,Y,z)

axis equal
view(2)
title('Tensiones respecto a gz')
colorbar

%Cálculo en 3D:

%Dibujo de las tensiones respecto a g sub rho 3D
subplot(2,3,4);

%Tomamos el valor (1,1) de la matriz sigma y lo dibujamos

surf(X,Y,r)

axis equal
title('Tensiones respecto a gRho 3D')

%Dibujo de las tensiones respecto a g sub tetha 3D
subplot(2,3,5);

%Tomamos el valor (2,2) de la matriz sigma y lo dibujamos

surf(X,Y,t)

axis equal
title('Tensiones respecto a gTheta 3D')

%Dibujo de las tensiones respecto a g sub z 3D
subplot(2,3,6)

%Tomamos el valor (3,3) de la matriz sigma y lo dibujamos

surf(X,Y,z)

axis equal
title('Tensiones respecto a gz 3D')
}}

[[Archivo:tensionesNormalesgrupo10A.png|centro|thumb|1300px]]

== Tensiones tangenciales ==

'''Ambas tensiones tangenciales son iguales, pues la matriz σ es simétrica.
Las tensiones del apartado nos dan: σ(2,1)= σ(2,1)= 1/ρ*(1/5*sen(2θ)-1)'''

== Tensiones tangenciales respecto al plano ortogonal a ḡρ ==

[[Archivo:Tensionestangencialesgrupo10A.png||thumb|390px|derecha]]

{{matlab|codigo=
%Inicializamos las variables
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos las tensiones tangentes al plano g sub rho
Tang=(RHO-1)./(5.*RHO.^2);

%Visualización en 3D
subplot(2,1,1)

surf(X,Y,Tang)

colorbar
title('Tensiones tangenciales respecto al plano ortogonal a gRho en 3D')

%Visualización en 2D
subplot(2,1,2)

surf(X,Y,Tang)

colorbar
view(2)
title('Tensiones tangenciales respecto al plano ortogonal a gRho en 2D')
}}

== Tensión de Von Mises ==

La tensión de Von Mises se define mediante la expresión:

[[Archivo:FormulaVonMises.png|marco|centro]]

La tensión de Von Mises es una magnitud escalar que se suele utilizar como indicador para saber cuando un material inicia un comportamiento plástico y no elástico puro. En nuestro caso, el valor máximo de la Tensión de Von Mises es igual a 0.6283.

[[Archivo:tensionVonMisesgrupo10A.png|marco|derecha]]

{{matlab|codigo=
%Inicializamos las variables
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos la parametrización
[RHO,THETA]=meshgrid(rho,theta);

%Creamos la matriz sigma
Sigma=zeros(3,3);

%Creamos la matriz de Von Mises
MatrizVonMises=zeros(16,21);

%Obtenemos los valores de la tensión de Von Mises a partir de los autovalores
for i=1:16*21

r=RHO(i)';

t=THETA(i)';

Sigma(1,1)=2*t/5+t*(2-1/r)/5;

Sigma(1,2)=(r-1)/(5*r^2);

Sigma(2,1)= (r-1)/(5*r^2);

Sigma(2,2)=2*t*(r-1)/(5*r^3)+t*(2-1/r)/(5);
Sigma(3,3)=t*(2-1/r)/5;

%Obtención de los autovalores de la matriz sigma
[v,d]=eig(Sigma);

%Fórmula de Von Mises
MatrizVonMises(i)=sqrt(((d(1,1)-d(2,2))^2+(d(2,2)-d(3,3))^2+(d(3,3)-d(1,1))^2)/2);

end

%Cambiamos a coordenadas cartesianas
X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos la tensión de Von Mises
surf(X,Y,MatrizVonMises)

colorbar
view(2)
title('Tensión de Von Mises')

%Valor Máximo
[i,j]=find(MatrizVonMises==max(max((MatrizVonMises))))
z=max(max(MatrizVonMises))

%El valor máximo de la Tensión de Von Mises es igual a 0.6283
}}

Visualización de campos escalares y vectoriales en elasticidad. Grupo 10-A

2016-12-01T17:00:50Z

Alejandro de Benito Andrés:

{{ TrabajoED | Visualizaciones de campos escalares y vectoriales. Grupo 10-A | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC16/17|2016-17]] | Alejandro de Benito, Lucía Rico , Adrián Piñeiro, Daniel Apellániz, Guillermo Alfaya. }}

EDITAR Consideramos una placa plana que ocupa el anillo circular centrado en el origen que ocupa el cuarto
anillo circular centrado en el origen y comprendido entre los radios 1 y 3, en el primer cuadrante
x, y ≥ 0. Dibujar con los ejes en el cuadrado [−1, 4]×[−1, 4]. Trabajaremos en coordenadas cilíndricas, en la cual definimos la temperatura T(x,y)=(y+2)·x2, dependiente de las variables (x,y), y el campo de desplazamientos u(ρ,θ)=θf(ρ)gρ. producido por la acción de una fuerza. De ellos, estudiaremos el efecto de la temperatura dada por \(T(x,y)\) y el estudio del movimiento de la placa con el campo \(\vec u(\rho,\theta)\), con las consiguientes tensiones que producirá.

'''Dada una superficie plana en forma de anillo circular, concretamente la cuarta parte situada en el primer cuadrante, con centro en el origen de radio exterior 3 e interior 1. Se han utilizado como ejes el cuadrado [-1,4]x[-1,4] para una mejor visualización de la placa. Por comodidad han sido elegidas las coordenadas cilíndricas, a la que se aplicará una temperatura T(x,y)=(y+2)·x2, en la que se toman como variables independientes "x, y" y el campo de desplazamientos u(ρ,θ)=θf(ρ)gρ, conclusión de un agente exterior desconocido. En apartados posteriores veremos como afectan la temperatura o el campo de desplazamientos a la placa en cuestión con las tensiones ejercidas.
'''
== Mallado ==
Representaremos una malla de los puntos que componen la placa sólida con forma de corona semicircular.
El paso de muestreo dado es 0,1 tanto para la variable X como para la variable Y.
[[Archivo:Malladogrupo10A.png|thumb|400px|derecha]]

{{matlab|codigo=
%Definimos las variables con el paso de muestreo "h"
h=0.1
rho=[1:h:3];
theta=[0:h:pi/2];

%Creamos el cojunto de puntos que forman la malla y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis([-1,4,-1,4]);

view(2);
title('Mallado')
}}

== Temperatura ==

La temperatura de la corona semicircular cambia según la función T(x,y)=(y+2)·x2: La función de únicamente depende de la variable Y y se mantendrá constante para todo valor de X.
Se ha dibujado el campo de temperaturas tanto en 2D como en 3D.
El punto que ha resultado tener mayor temperatura T(x,y)=24.1902 es el punto x=1.1683,y=2.7632.

'''En la corona la temperatura varía según la función dada, "T(x,y)=(y+2)·x2". Para tener una idea más realista de la temperatura la hemos dibujado tanto como en 2D como curvas de nivel. El punto que ha resultado tener mayor temperatura T(x,y)=24.1902 es el punto x=1.1683,y=2.7632.

'''
[[File:temperaturaGrupo10A2.png|thumb|800px|derecha]]

{{matlab|codigo=
%Definimos las variables con el paso de muestreo "h"
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H); %Así se ajusta mejor a theta=pi/2

%Creamos el cojunto de puntos que forman la malla y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis equal;

view(2);

%CAMPO ESCALAR DE TEMPERATURAS
Temperatura=(Y+2).*X.^2
Maximo=max(max(Temperatura))

%Líneas de programa usadas para obtener posición del máximo
for i=1:16
for j=1:21
Temp=(Y(i,j)+2).*X(i,j).^2;
if Temp==Maximo
disp(X(i,j))
disp(Y(i,j))
end
end
end

%Representación de la malla en 2D en la subventana primera
subplot(1,2,1)

contour(X,Y,Temperatura,20);

view(2)
axis([-1,4,-1,4]);
title('Curvas de nivel de la Temperatura')

%Representación de la malla en 3D en la subventana segunda
subplot(1,2,2);

surf(X,Y,Temperatura);

view(2)
axis([-1,4,-1,4]);
title('Temperatura en 2D')
}}

== Gradiente de la temperatura ==

El gradiente del campo de temperaturas (∇T) muestra la variación de la temperatura al movernos por el plano, siendo el módulo del gradiente lo que nos indica la rapidez con la que aumenta la temperatura.
El gradiente apunta en la dirección en la cual la temperatura aumenta más rápido: donde la derivada direccional será máxima.
El vector gradiente es perpendicular a la superficie del sólido.
Cuanto más aumenta la Y, mayor es la temperatura: En el dibujo podemos observar el aumento de la temperatura a medida que los puntos del sólido avanzan en el eje Y.

'''Para conocer la variación de temperatura por la placa usaremos el gradiente del campo de temperaturas(∇T), el módulo de los vectores gradiente representa la rapidez con la que varía la temperatura, es decir apunta hacia la dirección en la que la temperatura tiene un aumento más notable, donde la derivada direccional es máxima. Se aprecia en la figura que el vector gradiente en perpendicular a la superficie de sólido.'''

[[Archivo:gradienteDeLaTemperaturaGrupo10A1.png|thumb|800px|Izquierda: Gráfica de la temperatura en 2D y vectores gradientes. Derecha: Vectores del gradiente|derecha]]
{{matlab|codigo=
%Definimos las variables con el paso de muestreo "h"
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el cojunto de puntos que forman la malla y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis([-1,4,-1,4]);
axis square

view(2);

%GRADIENTE
GradX=2*X.*Y+4*X;
GradY=X.^2;

%CAMPO ESCALAR DE TEMPERATURAS
Temperatura=(Y+2).*X.^2;

%Creación de la subventana primera: gradiente en el campo de temperaturas

subplot(1,2,1)
hold on
surf(X,Y,Temperatura);
axis([-1,4,-1,4]);
axis square

%Dibujo del gradiente
quiver3(X,Y,Temperatura,GradX,GradY,0*Temperatura);
axis([-1,4,-1,4]);
axis square

hold off

%Creación de la segunda subventana: campo vectorial del gradiente

subplot(1,2,2)
quiver(X,Y,GradX,GradY);
axis([-1,4,-1,4]);
axis equal
colorbar
}}

== Campo de desplazamientos ==

El campo de desplazamientos dado es: u(ρ,θ)=θf(ρ)gρ. Vamos a dibujar el campo de vectorial, sabiendo que: Los puntos situados en ρ=1 no sufren desplazamiento y div(u)=θ(2-1/ρ)/5 ,
lo que nos da la ecuación diferencial f'(ρ)+(1/ρ)*f(ρ)=(2-1/ρ)/5 con la condición de contorno f(1)=0. Resultando f(ρ)=(ρ-1)/5.
Lo metemos en el campo de desplazamiento dando u(ρ,θ)=θ*((ρ-1)/5)gρ .Cogiendo el campo en función de los dos componente i,j. Definiremos pues FX y FY. Sabiendo que gρ= cosθi+senθj.

'''En este apartado procederemos a dibujar el campo de desplazamientos "u(ρ,θ)=θf(ρ)gρ". Las condiciones que conocemos son: los puntos de ρ=1 no sufren desplazamiento y div(u)=θ(2-1/ρ)/5, lo que nos da la ecuación diferencial f'(ρ)+(1/ρ)*f(ρ)=(2-1/ρ)/5 con la condición de contorno f(1)=0. Resultando f(ρ)=(ρ-1)/5. Lo metemos en el campo de desplazamiento dando u(ρ,θ)=θ*((ρ-1)/5)gρ .Cogiendo el campo en función de los dos componente i,j. Definiremos pues FX y FY. Sabiendo que gρ= cosθi+senθj.'''

[[Archivo:CampoDesplazamientosGrupo10A1.png|marco|derecha]]

{{matlab|codigo=
%Definimos las variables con el paso de muestreo "h"
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el cojunto de puntos que forman la malla y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis([-1,4,-1,4]);

view(2);

%CAMPO VECTORIAL DE DESPLAZAMIENTOS
%Componente x(X,Y) del campo vectorial
FX=THETA.*((RHO-1)/5).*(cos(THETA));

%Componente y(X,Y) del campo vectorial
FY=THETA.*((RHO-1)/5).*(sin(THETA));

%Representamos el campo vectorial de desplazamientos
quiver(X,Y,FX,FY)

axis([-1,4,-1,4])
view(2)

xlabel('Eje x')
ylabel('Eje y')
title('Campo de desplazamientos sobre la placa')
}}

== Movimiento del sólido ==

Se ha representado el sólido antes y después del desplazamiento dado por el campo vectorial para así poder compararlos. Definiremos las funciones de desplazamiento en función de ambas componentes, sumando las componentes del vector de desplazamientos a las coordenadas respectivas de posición (Las de la posición incial).
Se puede ver que el desplazamiento es mínimo.

'''En este caso se muestran dos gráficas una con la placa en reposo y otra después de aplicarle el campo de desplazamientos a la superficie, de esta forma es perfectamente apreciable la variación que sufre. El desplazamiento viene dado por dos componentes, el vector desplazamiento y los puntos de la placa respectivos a cada coordenada. Cuando ambas se suman, dan lugar a la figura final.
'''
[[Archivo:movimientoDelSolidoGrupo10A1.png|thumb|800px|derecha]]

{{matlab|codigo=
%Inicializamos las variables
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis([-1,4,-1,4]);

view(2);

%CAMPO VECTORIAL DE DESPLAZAMIENTOS
%Componente x(X,Y) del campo vectorial
FX=THETA.*((RHO-1)/5).*(cos(THETA));

%Componente y(X,Y) del campo vectorial
FY=THETA.*((RHO-1)/5).*(sin(THETA));

%Dibujo de la placa sin movimiento
subplot(1,2,1)

surf(X,Y,0*X);

axis([-1,4,-1,4]);
view(2);
title('Reposo')

%Dibujo de la placa en movimiento
%Movimiento de la placa en el eje X
MovimientoX=X+FX;

%Movimiento de la placa en el eje Y
MovimientoY=Y+FY;

subplot(1,2,2)

surf(MovimientoX,MovimientoY,Y.*0)

axis([-1,4,-1,4]);
view(2);
title('Aplicado el campo de movimientos')
}}

== Divergencia ==

La ∇·(u)=θ(2-1/ρ)/5 dada en el apartado anterior. Para calcular donde tiene los máximos,mínimos y nulos procederemos calculando su gradiente.
<big>∇(∇·u)=θ/(5*ρ²)gρ+(2-1/ρ)/5*(1/ρ²)gθ+0gz</big>

Como nos movemos entre ρ=(1,3) y θ=(0,pi/2) las derivada direccionales gρ y gθ son positivas, por lo que cuanto mayor sean ρ y θ mayor será su valor, alcanzándose el Máximo de ∇·u(ρ,θ)=0.5236 en (3,pi/2). Como ∇·(u) esta en función de θ, cuando θ=0 ∇·(u)=0 por lo que sera nulo y a su vez se alcanza su mínimo.

'''Dada la divergencia usada para definir el campo "∇·(u)=θ(2-1/ρ)/5", procederemos a evaluar sus máximos, mínimos y nulos con el gradiente<big>∇(∇·u)=θ/(5*ρ²)gρ+(2-1/ρ)/5*(1/ρ²)gθ+0gz</big>.
Al tener un dominio de ρ contenido en (1,3) y de θ en (0,pi/2), la derivada direccional gρ y gθ son positivas, también podemos sacar en claro que cuanto mayor sean ρ y θ mayor será su valor, llegando a su punto máximo en (3,pi/2) con: ∇·u(ρ,θ)=0.5236. A su vez el mínimo se encontrará donde los θ sean menores en nuestro dominio θ=0 dando lugar a: ∇·(u)=0 por estar ∇·(u) en función de θ'''

[[Archivo:DivergenciaGrupo10A.png|marco|derecha]]

{{matlab|codigo=
%Inicializamos las variables
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis([-1,4,-1,4]);

view(2);

%DIVERGENCIA
DIVERGENCIA=THETA.*((2-1./RHO)./5);

MaxDivergencia=max(max(DIVERGENCIA))

%Representación de la divergencia
surf(X,Y,DIVERGENCIA)

view(2)
colorbar
title('Divergencia')
}}

== Rotacional ==

Se ha calculado el rotacional del campo de desplazamientos u en todos los puntos del sólido y lo representamos. Hemos obtenido que el rotacional es: Rot [u(rho,thetha)] = (rho-1)/(5*rho)gz
Observando la gráfica se obtiene cuáles son los puntos que sufren un mayor rotacional: Los puntos con mayor rotacional son aquellos con rho=3 y su valor es 0.1333

El rotacional del campo de desplazamientos se ha calculado y aplicado a todos los puntos de la placa y han sido representados con colores. El ca

[[Archivo:RotacionalGrupo10A.png|marco|derecha]]

{{matlab|codigo=
%Inicializamos las variables
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis([-1,4,-1,4]);

view(2);

%ROTACIONAL
ROTACIONAL=(RHO-1)./(5.*RHO);

%Representación del rotacional
surf(X,Y,ROTACIONAL)

colorbar;
view(2)
title('Rotacional')

%Puntos con mayor rotacional
Maximo_Rotacional=max(max(ROTACIONAL))

%El máximo del rotacional corresponde a puntos con THETA = (PI/4)
title('Rotacional')

%Puntos con mayor rotacional
Maximo_Rotacional=max(max(ROTACIONAL))

}}

== Tensor de tensiones ==
Si consideramos que estamos en un medio elástico lineal, isótropo y homogéneo; en mecánica de medios continuos podemos definir el tensor de tensiones "σ"como:

[[Archivo:TensorDeTensionesGrupo10A.png|marco|centro]]

siendo λ y μ los coeficientes de Lamé, que dependen de las propiedades elásticas de cada material, considerando en nuestro caso que λ y μ son iguales a 1 y que ε es la parte simétrica del vector gradiente u.

Para hacer este apartado se requiere utilizar bases generales, en lo cual hemos llegado a la conclusión de que las componentes de ∇ū son:

[[Archivo:GradienteDeUGrupo10A.png|marco|centro]]

Siendo ε =(∇u + ∇ut)/2, con todo ello, calculamos la matriz σ:

[[Archivo:MatrizDeTensionesGrupo10A.png|marco|centro]]

Se van a representar las tensiones en las tres direcciones del espacio. En el caso de estar utilzando coordenadas cilíndricas, se representan las tensiones en las direcciones de la base cilíndrica ḡρ, ḡθ, ḡz. Para hallar las tensiones en estas direcciones debemos multiplicar escalarmente la matriz σ por cada uno de los vectores ḡi por ambos lados: Tensión normal a la dirección ḡi = ḡi*σ*ḡi

[[Archivo:DireccionGrupo10A.png|marco|centro]]

{{matlab|codigo=
%Inicializamos las variables
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%En 2D

%Dibujo de las tensiones respecto a g sub rho
subplot(2,3,1);

%Tomamos el valor (1,1) de la matriz sigma y lo dibujamos
r=2.*THETA./5+THETA.*(2-1./RHO)./5;

surf(X,Y,r)

axis equal
view(2)
title('Tensiones respecto a gRho')
colorbar

%Dibujo de las tensiones respecto a g sub tetha
subplot(2,3,2);

%Tomamos el valor (2,2) de la matriz sigma y lo dibujamos
t=2.*THETA.*(RHO-1)./(5.*RHO.^5)+THETA.*(2-1./RHO)./(5.*RHO.^2);

surf(X,Y,t)

axis equal
title('Tensiones respecto a gTheta')
view(2)

colorbar

%Dibujo de las tensiones respecto a g sub z
subplot(2,3,3)

%Tomamos el valor (3,3) de la matriz sigma y lo dibujamos
z=THETA.*(2-1./RHO)./5;

surf(X,Y,z)

axis equal
view(2)
title('Tensiones respecto a gz')
colorbar

%Cálculo en 3D:

%Dibujo de las tensiones respecto a g sub rho 3D
subplot(2,3,4);

%Tomamos el valor (1,1) de la matriz sigma y lo dibujamos

surf(X,Y,r)

axis equal
title('Tensiones respecto a gRho 3D')

%Dibujo de las tensiones respecto a g sub tetha 3D
subplot(2,3,5);

%Tomamos el valor (2,2) de la matriz sigma y lo dibujamos

surf(X,Y,t)

axis equal
title('Tensiones respecto a gTheta 3D')

%Dibujo de las tensiones respecto a g sub z 3D
subplot(2,3,6)

%Tomamos el valor (3,3) de la matriz sigma y lo dibujamos

surf(X,Y,z)

axis equal
title('Tensiones respecto a gz 3D')
}}

[[Archivo:tensionesNormalesgrupo10A.png|centro|thumb|1300px]]

== Tensiones tangenciales ==

Ambas tensiones tangenciales son iguales, pues la matriz σ es simétrica.
Las tensiones del apartado nos dan: σ(2,1)= σ(2,1)= 1/ρ*(1/5*sen(2θ)-1)

== Tensiones tangenciales respecto al plano ortogonal a ḡρ ==

[[Archivo:Tensionestangencialesgrupo10A.png||thumb|390px|derecha]]

{{matlab|codigo=
%Inicializamos las variables
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos las tensiones tangentes al plano g sub rho
Tang=(RHO-1)./(5.*RHO.^2);

%Visualización en 3D
subplot(2,1,1)

surf(X,Y,Tang)

colorbar
title('Tensiones tangenciales respecto al plano ortogonal a gRho en 3D')

%Visualización en 2D
subplot(2,1,2)

surf(X,Y,Tang)

colorbar
view(2)
title('Tensiones tangenciales respecto al plano ortogonal a gRho en 2D')
}}

== Tensión de Von Mises ==

La tensión de Von Mises se define mediante la expresión:

[[Archivo:FormulaVonMises.png|marco|centro]]

La tensión de Von Mises es una magnitud escalar que se suele utilizar como indicador para saber cuando un material inicia un comportamiento plástico y no elástico puro. En nuestro caso, el valor máximo de la Tensión de Von Mises es igual a 2,0783.

[[Archivo:tensionVonMisesgrupo10A.png|marco|derecha]]

{{matlab|codigo=
%Inicializamos las variables
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos la parametrización
[RHO,THETA]=meshgrid(rho,theta);

%Creamos la matriz sigma
Sigma=zeros(3,3);

%Creamos la matriz de Von Mises
MatrizVonMises=zeros(16,21);

%Obtenemos los valores de la tensión de Von Mises a partir de los autovalores
for i=1:16*21

r=RHO(i)';

t=THETA(i)';

Sigma(1,1)=2*t/5+t*(2-1/r)/5;

Sigma(1,2)=(r-1)/(5*r^2);

Sigma(2,1)= (r-1)/(5*r^2);

Sigma(2,2)=2*t*(r-1)/(5*r^3)+t*(2-1/r)/(5);
Sigma(3,3)=t*(2-1/r)/5;

%Obtención de los autovalores de la matriz sigma
[v,d]=eig(Sigma);

%Fórmula de Von Mises
MatrizVonMises(i)=sqrt(((d(1,1)-d(2,2))^2+(d(2,2)-d(3,3))^2+(d(3,3)-d(1,1))^2)/2);

end

%Cambiamos a coordenadas cartesianas
X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos la tensión de Von Mises
surf(X,Y,MatrizVonMises)

colorbar
view(2)
title('Tensión de Von Mises')

%Valor Máximo
[i,j]=find(MatrizVonMises==max(max((MatrizVonMises))))
z=max(max(MatrizVonMises))

%El valor máximo de la Tensión de Von Mises es igual a 0.6283
}}

Visualización de campos escalares y vectoriales en elasticidad. Grupo 10-A

2016-12-01T16:19:02Z

Alejandro de Benito Andrés:

{{ TrabajoED | Visualizaciones de campos escalares y vectoriales. Grupo 10-A | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC16/17|2016-17]] | Alejandro de Benito, Lucía Rico , Adrián Piñeiro, Daniel Apellániz, Guillermo Alfaya. }}

EDITAR Consideramos una placa plana que ocupa el anillo circular centrado en el origen que ocupa el cuarto
anillo circular centrado en el origen y comprendido entre los radios 1 y 3, en el primer cuadrante
x, y ≥ 0. Dibujar con los ejes en el cuadrado [−1, 4]×[−1, 4]. Trabajaremos en coordenadas cilíndricas, en la cual definimos la temperatura T(x,y)=(y+2)·x2, dependiente de las variables (x,y), y el campo de desplazamientos u(ρ,θ)=θf(ρ)gρ. producido por la acción de una fuerza. De ellos, estudiaremos el efecto de la temperatura dada por \(T(x,y)\) y el estudio del movimiento de la placa con el campo \(\vec u(\rho,\theta)\), con las consiguientes tensiones que producirá.

'''Dada una superficie plana en forma de anillo circular, concretamente la cuarta parte situada en el primer cuadrante, con centro en el origen de radio exterior 3 e interior 1. Se han utilizado como ejes el cuadrado [-1,4]x[-1,4] para una mejor visualización de la placa. Por comodidad han sido elegidas las coordenadas cilíndricas, a la que se aplicará una temperatura T(x,y)=(y+2)·x2, en la que se toman como variables independientes "x, y" y el campo de desplazamientos u(ρ,θ)=θf(ρ)gρ, conclusión de un agente exterior desconocido. En apartados posteriores veremos como afectan la temperatura o el campo de desplazamientos a la placa en cuestión con las tensiones ejercidas.
'''
== Mallado ==
Representaremos una malla de los puntos que componen la placa sólida con forma de corona semicircular.
El paso de muestreo dado es 0,1 tanto para la variable X como para la variable Y.
[[Archivo:Malladogrupo10A.png|thumb|400px|derecha]]

{{matlab|codigo=
%Definimos las variables con el paso de muestreo "h"
h=0.1
rho=[1:h:3];
theta=[0:h:pi/2];

%Creamos el cojunto de puntos que forman la malla y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis([-1,4,-1,4]);

view(2);
title('Mallado')
}}

== Temperatura ==

La temperatura de la corona semicircular cambia según la función T(x,y)=(y+2)·x2: La función de únicamente depende de la variable Y y se mantendrá constante para todo valor de X.
Se ha dibujado el campo de temperaturas tanto en 2D como en 3D.
El punto que ha resultado tener mayor temperatura T(x,y)=24.1902 es el punto x=1.1683,y=2.7632.

'''En la corona la temperatura varía según la función dada, "T(x,y)=(y+2)·x2". Para tener una idea más realista de la temperatura la hemos dibujado tanto como en 2D como curvas de nivel. El punto que ha resultado tener mayor temperatura T(x,y)=24.1902 es el punto x=1.1683,y=2.7632.

'''
[[File:temperaturaGrupo10A2.png|thumb|800px|derecha]]

{{matlab|codigo=
%Definimos las variables con el paso de muestreo "h"
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H); %Así se ajusta mejor a theta=pi/2

%Creamos el cojunto de puntos que forman la malla y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis equal;

view(2);

%CAMPO ESCALAR DE TEMPERATURAS
Temperatura=(Y+2).*X.^2
Maximo=max(max(Temperatura))

%Líneas de programa usadas para obtener posición del máximo
for i=1:16
for j=1:21
Temp=(Y(i,j)+2).*X(i,j).^2;
if Temp==Maximo
disp(X(i,j))
disp(Y(i,j))
end
end
end

%Representación de la malla en 2D en la subventana primera
subplot(1,2,1)

contour(X,Y,Temperatura,20);

view(2)
axis([-1,4,-1,4]);
title('Curvas de nivel de la Temperatura')

%Representación de la malla en 3D en la subventana segunda
subplot(1,2,2);

surf(X,Y,Temperatura);

view(2)
axis([-1,4,-1,4]);
title('Temperatura en 2D')
}}

== Gradiente de la temperatura ==

El gradiente del campo de temperaturas (∇T) muestra la variación de la temperatura al movernos por el plano, siendo el módulo del gradiente lo que nos indica la rapidez con la que aumenta la temperatura.
El gradiente apunta en la dirección en la cual la temperatura aumenta más rápido: donde la derivada direccional será máxima.
El vector gradiente es perpendicular a la superficie del sólido.
Cuanto más aumenta la Y, mayor es la temperatura: En el dibujo podemos observar el aumento de la temperatura a medida que los puntos del sólido avanzan en el eje Y.

'''Para conocer la variación de temperatura por la placa usaremos el gradiente del campo de temperaturas(∇T), el módulo de los vectores gradiente representa la rapidez con la que varía la temperatura, es decir apunta hacia la dirección en la que la temperatura tiene un aumento más notable, donde la derivada direccional es máxima. Se aprecia en la figura que el vector gradiente en perpendicular a la superficie de sólido.'''

[[Archivo:gradienteDeLaTemperaturaGrupo10A1.png|thumb|800px|Izquierda: Gráfica de la temperatura en 2D y vectores gradientes. Derecha: Vectores del gradiente|derecha]]
{{matlab|codigo=
%Definimos las variables con el paso de muestreo "h"
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el cojunto de puntos que forman la malla y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis([-1,4,-1,4]);
axis square

view(2);

%GRADIENTE
GradX=2*X.*Y+4*X;
GradY=X.^2;

%CAMPO ESCALAR DE TEMPERATURAS
Temperatura=(Y+2).*X.^2;

%Creación de la subventana primera: gradiente en el campo de temperaturas

subplot(1,2,1)
hold on
surf(X,Y,Temperatura);
axis([-1,4,-1,4]);
axis square

%Dibujo del gradiente
quiver3(X,Y,Temperatura,GradX,GradY,0*Temperatura);
axis([-1,4,-1,4]);
axis square

hold off

%Creación de la segunda subventana: campo vectorial del gradiente

subplot(1,2,2)
quiver(X,Y,GradX,GradY);
axis([-1,4,-1,4]);
axis equal
colorbar
}}

== Campo de desplazamientos ==

El campo de desplazamientos dado es: u(ρ,θ)=θf(ρ)gρ. Vamos a dibujar el campo de vectorial, sabiendo que: Los puntos situados en ρ=1 no sufren desplazamiento y div(u)=θ(2-1/ρ)/5 ,
lo que nos da la ecuación diferencial f'(ρ)+(1/ρ)*f(ρ)=(2-1/ρ)/5 con la condición de contorno f(1)=0. Resultando f(ρ)=(ρ-1)/5.
Lo metemos en el campo de desplazamiento dando u(ρ,θ)=θ*((ρ-1)/5)gρ .Cogiendo el campo en función de los dos componente i,j. Definiremos pues FX y FY. Sabiendo que gρ= cosθi+senθj.

'''En este apartado procederemos a dibujar el campo de desplazamientos "u(ρ,θ)=θf(ρ)gρ". Las condiciones que conocemos son: los puntos de ρ=1 no sufren desplazamiento y div(u)=θ(2-1/ρ)/5, lo que nos da la ecuación diferencial f'(ρ)+(1/ρ)*f(ρ)=(2-1/ρ)/5 con la condición de contorno f(1)=0. Resultando f(ρ)=(ρ-1)/5. Lo metemos en el campo de desplazamiento dando u(ρ,θ)=θ*((ρ-1)/5)gρ .Cogiendo el campo en función de los dos componente i,j. Definiremos pues FX y FY. Sabiendo que gρ= cosθi+senθj.'''

[[Archivo:CampoDesplazamientosGrupo10A1.png|marco|derecha]]

{{matlab|codigo=
%Definimos las variables con el paso de muestreo "h"
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el cojunto de puntos que forman la malla y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis([-1,4,-1,4]);

view(2);

%CAMPO VECTORIAL DE DESPLAZAMIENTOS
%Componente x(X,Y) del campo vectorial
FX=THETA.*((RHO-1)/5).*(cos(THETA));

%Componente y(X,Y) del campo vectorial
FY=THETA.*((RHO-1)/5).*(sin(THETA));

%Representamos el campo vectorial de desplazamientos
quiver(X,Y,FX,FY)

axis([-1,4,-1,4])
view(2)

xlabel('Eje x')
ylabel('Eje y')
title('Campo de desplazamientos sobre la placa')
}}

== Movimiento del sólido ==

Se ha representado el sólido antes y después del desplazamiento dado por el campo vectorial para así poder compararlos. Definiremos las funciones de desplazamiento en función de ambas componentes, sumando las componentes del vector de desplazamientos a las coordenadas respectivas de posición (Las de la posición incial).
Se puede ver que el desplazamiento es mínimo.

[[Archivo:movimientoDelSolidoGrupo10A1.png|thumb|800px|derecha]]

{{matlab|codigo=
%Inicializamos las variables
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis([-1,4,-1,4]);

view(2);

%CAMPO VECTORIAL DE DESPLAZAMIENTOS
%Componente x(X,Y) del campo vectorial
FX=THETA.*((RHO-1)/5).*(cos(THETA));

%Componente y(X,Y) del campo vectorial
FY=THETA.*((RHO-1)/5).*(sin(THETA));

%Dibujo de la placa sin movimiento
subplot(1,2,1)

surf(X,Y,0*X);

axis([-1,4,-1,4]);
view(2);
title('Reposo')

%Dibujo de la placa en movimiento
%Movimiento de la placa en el eje X
MovimientoX=X+FX;

%Movimiento de la placa en el eje Y
MovimientoY=Y+FY;

subplot(1,2,2)

surf(MovimientoX,MovimientoY,Y.*0)

axis([-1,4,-1,4]);
view(2);
title('Aplicado el campo de movimientos')
}}

== Divergencia ==

La ∇·(u)=θ(2-1/ρ)/5 dada en el apartado anterior. Para calcular donde tiene los máximos,mínimos y nulos procederemos calculando su gradiente.
<big>∇(∇·u)=θ/(5*ρ²)gρ+(2-1/ρ)/5*(1/ρ²)gθ+0gz</big>

Como nos movemos entre ρ=(1,3) y θ=(0,pi/2) las derivada direccionales gρ y gθ son positivas, por lo que cuanto mayor sean ρ y θ mayor será su valor, alcanzándose el Máximo de ∇·u(ρ,θ)=0.5236 en (3,pi/2). Como ∇·(u) esta en función de θ, cuando θ=0 ∇·(u)=0 por lo que sera nulo y a su vez se alcanza su mínimo.

[[Archivo:DivergenciaGrupo10A.png|marco|derecha]]

{{matlab|codigo=
%Inicializamos las variables
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis([-1,4,-1,4]);

view(2);

%DIVERGENCIA
DIVERGENCIA=THETA.*((2-1./RHO)./5);

MaxDivergencia=max(max(DIVERGENCIA))

%Representación de la divergencia
surf(X,Y,DIVERGENCIA)

view(2)
colorbar
title('Divergencia')
}}

== Rotacional ==

Se ha calculado el rotacional del campo de desplazamientos u en todos los puntos del sólido y lo representamos. Hemos obtenido que el rotacional es: Rot [u(rho,thetha)] = (rho-1)/(5*rho)gz
Observando la gráfica se obtiene cuáles son los puntos que sufren un mayor rotacional: Los puntos con mayor rotacional son aquellos con rho=3 y su valor es 0.1333
[[Archivo:RotacionalGrupo10A.png|marco|derecha]]

{{matlab|codigo=
%Inicializamos las variables
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis([-1,4,-1,4]);

view(2);

%ROTACIONAL
ROTACIONAL=(RHO-1)./(5.*RHO);

%Representación del rotacional
surf(X,Y,ROTACIONAL)

colorbar;
view(2)
title('Rotacional')

%Puntos con mayor rotacional
Maximo_Rotacional=max(max(ROTACIONAL))

%El máximo del rotacional corresponde a puntos con THETA = (PI/4)
title('Rotacional')

%Puntos con mayor rotacional
Maximo_Rotacional=max(max(ROTACIONAL))

}}

== Tensor de tensiones ==
Si consideramos que estamos en un medio elástico lineal, isótropo y homogéneo; en mecánica de medios continuos podemos definir el tensor de tensiones "σ"como:

[[Archivo:TensorDeTensionesGrupo10A.png|marco|centro]]

siendo λ y μ los coeficientes de Lamé, que dependen de las propiedades elásticas de cada material, considerando en nuestro caso que λ y μ son iguales a 1 y que ε es la parte simétrica del vector gradiente u.

Para hacer este apartado se requiere utilizar bases generales, en lo cual hemos llegado a la conclusión de que las componentes de ∇ū son:

[[Archivo:GradienteDeUGrupo10A.png|marco|centro]]

Siendo ε =(∇u + ∇ut)/2, con todo ello, calculamos la matriz σ:

[[Archivo:MatrizDeTensionesGrupo10A.png|marco|centro]]

Se van a representar las tensiones en las tres direcciones del espacio. En el caso de estar utilzando coordenadas cilíndricas, se representan las tensiones en las direcciones de la base cilíndrica ḡρ, ḡθ, ḡz. Para hallar las tensiones en estas direcciones debemos multiplicar escalarmente la matriz σ por cada uno de los vectores ḡi por ambos lados: Tensión normal a la dirección ḡi = ḡi*σ*ḡi

[[Archivo:DireccionGrupo10A.png|marco|centro]]

{{matlab|codigo=
%Inicializamos las variables
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%En 2D

%Dibujo de las tensiones respecto a g sub rho
subplot(2,3,1);

%Tomamos el valor (1,1) de la matriz sigma y lo dibujamos
r=2.*THETA./5+THETA.*(2-1./RHO)./5;

surf(X,Y,r)

axis equal
view(2)
title('Tensiones respecto a gRho')
colorbar

%Dibujo de las tensiones respecto a g sub tetha
subplot(2,3,2);

%Tomamos el valor (2,2) de la matriz sigma y lo dibujamos
t=2.*THETA.*(RHO-1)./(5.*RHO.^5)+THETA.*(2-1./RHO)./(5.*RHO.^2);

surf(X,Y,t)

axis equal
title('Tensiones respecto a gTheta')
view(2)

colorbar

%Dibujo de las tensiones respecto a g sub z
subplot(2,3,3)

%Tomamos el valor (3,3) de la matriz sigma y lo dibujamos
z=THETA.*(2-1./RHO)./5;

surf(X,Y,z)

axis equal
view(2)
title('Tensiones respecto a gz')
colorbar

%Cálculo en 3D:

%Dibujo de las tensiones respecto a g sub rho 3D
subplot(2,3,4);

%Tomamos el valor (1,1) de la matriz sigma y lo dibujamos

surf(X,Y,r)

axis equal
title('Tensiones respecto a gRho 3D')

%Dibujo de las tensiones respecto a g sub tetha 3D
subplot(2,3,5);

%Tomamos el valor (2,2) de la matriz sigma y lo dibujamos

surf(X,Y,t)

axis equal
title('Tensiones respecto a gTheta 3D')

%Dibujo de las tensiones respecto a g sub z 3D
subplot(2,3,6)

%Tomamos el valor (3,3) de la matriz sigma y lo dibujamos

surf(X,Y,z)

axis equal
title('Tensiones respecto a gz 3D')
}}

[[Archivo:tensionesNormalesgrupo10A.png|centro|thumb|1300px]]

== Tensiones tangenciales ==

Ambas tensiones tangenciales son iguales, pues la matriz σ es simétrica.
Las tensiones del apartado nos dan: σ(2,1)= σ(2,1)= 1/ρ*(1/5*sen(2θ)-1)

== Tensiones tangenciales respecto al plano ortogonal a ḡρ ==

[[Archivo:Tensionestangencialesgrupo10A.png||thumb|390px|derecha]]

{{matlab|codigo=
%Inicializamos las variables
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos las tensiones tangentes al plano g sub rho
Tang=(RHO-1)./(5.*RHO.^2);

%Visualización en 3D
subplot(2,1,1)

surf(X,Y,Tang)

colorbar
title('Tensiones tangenciales respecto al plano ortogonal a gRho en 3D')

%Visualización en 2D
subplot(2,1,2)

surf(X,Y,Tang)

colorbar
view(2)
title('Tensiones tangenciales respecto al plano ortogonal a gRho en 2D')
}}

== Tensión de Von Mises ==

La tensión de Von Mises se define mediante la expresión:

[[Archivo:FormulaVonMises.png|marco|centro]]

La tensión de Von Mises es una magnitud escalar que se suele utilizar como indicador para saber cuando un material inicia un comportamiento plástico y no elástico puro. En nuestro caso, el valor máximo de la Tensión de Von Mises es igual a 2,0783.

[[Archivo:tensionVonMisesgrupo10A.png|marco|derecha]]

{{matlab|codigo=
%Inicializamos las variables
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos la parametrización
[RHO,THETA]=meshgrid(rho,theta);

%Creamos la matriz sigma
Sigma=zeros(3,3);

%Creamos la matriz de Von Mises
MatrizVonMises=zeros(16,21);

%Obtenemos los valores de la tensión de Von Mises a partir de los autovalores
for i=1:16*21

r=RHO(i)';

t=THETA(i)';

Sigma(1,1)=2*t/5+t*(2-1/r)/5;

Sigma(1,2)=(r-1)/(5*r^2);

Sigma(2,1)= (r-1)/(5*r^2);

Sigma(2,2)=2*t*(r-1)/(5*r^3)+t*(2-1/r)/(5);
Sigma(3,3)=t*(2-1/r)/5;

%Obtención de los autovalores de la matriz sigma
[v,d]=eig(Sigma);

%Fórmula de Von Mises
MatrizVonMises(i)=sqrt(((d(1,1)-d(2,2))^2+(d(2,2)-d(3,3))^2+(d(3,3)-d(1,1))^2)/2);

end

%Cambiamos a coordenadas cartesianas
X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos la tensión de Von Mises
surf(X,Y,MatrizVonMises)

colorbar
view(2)
title('Tensión de Von Mises')

%Valor Máximo
[i,j]=find(MatrizVonMises==max(max((MatrizVonMises))))
z=max(max(MatrizVonMises))

%El valor máximo de la Tensión de Von Mises es igual a 0.6283
}}

Visualización de campos escalares y vectoriales en elasticidad. Grupo 10-A

2016-12-01T16:16:08Z

Alejandro de Benito Andrés:

{{ TrabajoED | Visualizaciones de campos escalares y vectoriales. Grupo 10-A | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC16/17|2016-17]] | Alejandro de Benito, Lucía Rico , Adrián Piñeiro, Daniel Apellániz, Guillermo Alfaya. }}

EDITAR Consideramos una placa plana que ocupa el anillo circular centrado en el origen que ocupa el cuarto
anillo circular centrado en el origen y comprendido entre los radios 1 y 3, en el primer cuadrante
x, y ≥ 0. Dibujar con los ejes en el cuadrado [−1, 4]×[−1, 4]. Trabajaremos en coordenadas cilíndricas, en la cual definimos la temperatura T(x,y)=(y+2)·x2, dependiente de las variables (x,y), y el campo de desplazamientos u(ρ,θ)=θf(ρ)gρ. producido por la acción de una fuerza. De ellos, estudiaremos el efecto de la temperatura dada por \(T(x,y)\) y el estudio del movimiento de la placa con el campo \(\vec u(\rho,\theta)\), con las consiguientes tensiones que producirá.

'''Dada una superficie plana en forma de anillo circular, concretamente la cuarta parte situada en el primer cuadrante, con centro en el origen de radio exterior 3 e interior 1. Se han utilizado como ejes el cuadrado [-1,4]x[-1,4] para una mejor visualización de la placa. Por comodidad han sido elegidas las coordenadas cilíndricas, a la que se aplicará una temperatura T(x,y)=(y+2)·x2, en la que se toman como variables independientes "x, y" y el campo de desplazamientos u(ρ,θ)=θf(ρ)gρ, conclusión de un agente exterior desconocido. En apartados posteriores veremos como afectan la temperatura o el campo de desplazamientos a la placa en cuestión con las tensiones ejercidas.
'''
== Mallado ==
Representaremos una malla de los puntos que componen la placa sólida con forma de corona semicircular.
El paso de muestreo dado es 0,1 tanto para la variable X como para la variable Y.
[[Archivo:Malladogrupo10A.png|thumb|400px|derecha]]

{{matlab|codigo=
%Definimos las variables con el paso de muestreo "h"
h=0.1
rho=[1:h:3];
theta=[0:h:pi/2];

%Creamos el cojunto de puntos que forman la malla y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis([-1,4,-1,4]);

view(2);
title('Mallado')
}}

== Temperatura ==

La temperatura de la corona semicircular cambia según la función T(x,y)=(y+2)·x2: La función de únicamente depende de la variable Y y se mantendrá constante para todo valor de X.
Se ha dibujado el campo de temperaturas tanto en 2D como en 3D.
El punto que ha resultado tener mayor temperatura T(x,y)=24.1902 es el punto x=1.1683,y=2.7632.

'''En la corona la temperatura varía según la función dada, "T(x,y)=(y+2)·x2". Para tener una idea más realista de la temperatura la hemos dibujado tanto como en 2D como curvas de nivel. El punto que ha resultado tener mayor temperatura T(x,y)=24.1902 es el punto x=1.1683,y=2.7632.

'''
[[File:temperaturaGrupo10A2.png|thumb|800px|derecha]]

{{matlab|codigo=
%Definimos las variables con el paso de muestreo "h"
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H); %Así se ajusta mejor a theta=pi/2

%Creamos el cojunto de puntos que forman la malla y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis equal;

view(2);

%CAMPO ESCALAR DE TEMPERATURAS
Temperatura=(Y+2).*X.^2
Maximo=max(max(Temperatura))

%Líneas de programa usadas para obtener posición del máximo
for i=1:16
for j=1:21
Temp=(Y(i,j)+2).*X(i,j).^2;
if Temp==Maximo
disp(X(i,j))
disp(Y(i,j))
end
end
end

%Representación de la malla en 2D en la subventana primera
subplot(1,2,1)

contour(X,Y,Temperatura,20);

view(2)
axis([-1,4,-1,4]);
title('Curvas de nivel de la Temperatura')

%Representación de la malla en 3D en la subventana segunda
subplot(1,2,2);

surf(X,Y,Temperatura);

view(2)
axis([-1,4,-1,4]);
title('Temperatura en 2D')
}}

== Gradiente de la temperatura ==

El gradiente del campo de temperaturas (∇T) muestra la variación de la temperatura al movernos por el plano, siendo el módulo del gradiente lo que nos indica la rapidez con la que aumenta la temperatura.
El gradiente apunta en la dirección en la cual la temperatura aumenta más rápido: donde la derivada direccional será máxima.
El vector gradiente es perpendicular a la superficie del sólido.
Cuanto más aumenta la Y, mayor es la temperatura: En el dibujo podemos observar el aumento de la temperatura a medida que los puntos del sólido avanzan en el eje Y.

'''Para conocer la variación de temperatura por la placa usaremos el gradiente del campo de temperaturas(∇T), el módulo de los vectores gradiente representa la rapidez con la que varía la temperatura, es decir apunta hacia la dirección en la que la temperatura tiene un aumento más notable, donde la derivada direccional es máxima. Se aprecia en la figura que el vector gradiente en perpendicular a la superficie de sólido.'''

[[Archivo:gradienteDeLaTemperaturaGrupo10A1.png|thumb|800px|Izquierda: Gráfica de la temperatura en 2D y vectores gradientes. Derecha: Vectores del gradiente|derecha]]
{{matlab|codigo=
%Definimos las variables con el paso de muestreo "h"
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el cojunto de puntos que forman la malla y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis([-1,4,-1,4]);
axis square

view(2);

%GRADIENTE
GradX=2*X.*Y+4*X;
GradY=X.^2;

%CAMPO ESCALAR DE TEMPERATURAS
Temperatura=(Y+2).*X.^2;

%Creación de la subventana primera: gradiente en el campo de temperaturas

subplot(1,2,1)
hold on
surf(X,Y,Temperatura);
axis([-1,4,-1,4]);
axis square

%Dibujo del gradiente
quiver3(X,Y,Temperatura,GradX,GradY,0*Temperatura);
axis([-1,4,-1,4]);
axis square

hold off

%Creación de la segunda subventana: campo vectorial del gradiente

subplot(1,2,2)
quiver(X,Y,GradX,GradY);
axis([-1,4,-1,4]);
axis equal
colorbar
}}

== Campo de desplazamientos ==

El campo de desplazamientos dado es: u(ρ,θ)=θf(ρ)gρ. Vamos a dibujar el campo de vectorial, sabiendo que: Los puntos situados en ρ=1 no sufren desplazamiento y div(u)=θ(2-1/ρ)/5 ,
lo que nos da la ecuación diferencial f'(ρ)+(1/ρ)*f(ρ)=(2-1/ρ)/5 con la condición de contorno f(1)=0. Resultando f(ρ)=(ρ-1)/5.
Lo metemos en el campo de desplazamiento dando u(ρ,θ)=θ*((ρ-1)/5)gρ .Cogiendo el campo en función de los dos componente i,j. Definiremos pues FX y FY. Sabiendo que gρ= cosθi+senθj.

'''En este apartado procederemos a dibujar el campo de desplazamientos "u(ρ,θ)=θf(ρ)gρ". Las condiciones que conocemos son: los puntos de ρ=1 no sufren desplazamiento y div(u)=θ(2-1/ρ)/5 ,
lo que nos da la ecuación diferencial f'(ρ)+(1/ρ)*f(ρ)=(2-1/ρ)/5 con la condición de contorno f(1)=0. Resultando f(ρ)=(ρ-1)/5.
Lo metemos en el campo de desplazamiento dando u(ρ,θ)=θ*((ρ-1)/5)gρ .Cogiendo el campo en función de los dos componente i,j. Definiremos pues FX y FY. Sabiendo que gρ= cosθi+senθj.'''

[[Archivo:CampoDesplazamientosGrupo10A1.png|marco|derecha]]

{{matlab|codigo=
%Definimos las variables con el paso de muestreo "h"
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el cojunto de puntos que forman la malla y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis([-1,4,-1,4]);

view(2);

%CAMPO VECTORIAL DE DESPLAZAMIENTOS
%Componente x(X,Y) del campo vectorial
FX=THETA.*((RHO-1)/5).*(cos(THETA));

%Componente y(X,Y) del campo vectorial
FY=THETA.*((RHO-1)/5).*(sin(THETA));

%Representamos el campo vectorial de desplazamientos
quiver(X,Y,FX,FY)

axis([-1,4,-1,4])
view(2)

xlabel('Eje x')
ylabel('Eje y')
title('Campo de desplazamientos sobre la placa')
}}

== Movimiento del sólido ==

Se ha representado el sólido antes y después del desplazamiento dado por el campo vectorial para así poder compararlos. Definiremos las funciones de desplazamiento en función de ambas componentes, sumando las componentes del vector de desplazamientos a las coordenadas respectivas de posición (Las de la posición incial).
Se puede ver que el desplazamiento es mínimo.

[[Archivo:movimientoDelSolidoGrupo10A1.png|thumb|800px|derecha]]

{{matlab|codigo=
%Inicializamos las variables
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis([-1,4,-1,4]);

view(2);

%CAMPO VECTORIAL DE DESPLAZAMIENTOS
%Componente x(X,Y) del campo vectorial
FX=THETA.*((RHO-1)/5).*(cos(THETA));

%Componente y(X,Y) del campo vectorial
FY=THETA.*((RHO-1)/5).*(sin(THETA));

%Dibujo de la placa sin movimiento
subplot(1,2,1)

surf(X,Y,0*X);

axis([-1,4,-1,4]);
view(2);
title('Reposo')

%Dibujo de la placa en movimiento
%Movimiento de la placa en el eje X
MovimientoX=X+FX;

%Movimiento de la placa en el eje Y
MovimientoY=Y+FY;

subplot(1,2,2)

surf(MovimientoX,MovimientoY,Y.*0)

axis([-1,4,-1,4]);
view(2);
title('Aplicado el campo de movimientos')
}}

== Divergencia ==

La ∇·(u)=θ(2-1/ρ)/5 dada en el apartado anterior. Para calcular donde tiene los máximos,mínimos y nulos procederemos calculando su gradiente.
<big>∇(∇·u)=θ/(5*ρ²)gρ+(2-1/ρ)/5*(1/ρ²)gθ+0gz</big>

Como nos movemos entre ρ=(1,3) y θ=(0,pi/2) las derivada direccionales gρ y gθ son positivas, por lo que cuanto mayor sean ρ y θ mayor será su valor, alcanzándose el Máximo de ∇·u(ρ,θ)=0.5236 en (3,pi/2). Como ∇·(u) esta en función de θ, cuando θ=0 ∇·(u)=0 por lo que sera nulo y a su vez se alcanza su mínimo.

[[Archivo:DivergenciaGrupo10A.png|marco|derecha]]

{{matlab|codigo=
%Inicializamos las variables
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis([-1,4,-1,4]);

view(2);

%DIVERGENCIA
DIVERGENCIA=THETA.*((2-1./RHO)./5);

MaxDivergencia=max(max(DIVERGENCIA))

%Representación de la divergencia
surf(X,Y,DIVERGENCIA)

view(2)
colorbar
title('Divergencia')
}}

== Rotacional ==

Se ha calculado el rotacional del campo de desplazamientos u en todos los puntos del sólido y lo representamos. Hemos obtenido que el rotacional es: Rot [u(rho,thetha)] = (rho-1)/(5*rho)gz
Observando la gráfica se obtiene cuáles son los puntos que sufren un mayor rotacional: Los puntos con mayor rotacional son aquellos con rho=3 y su valor es 0.1333
[[Archivo:RotacionalGrupo10A.png|marco|derecha]]

{{matlab|codigo=
%Inicializamos las variables
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis([-1,4,-1,4]);

view(2);

%ROTACIONAL
ROTACIONAL=(RHO-1)./(5.*RHO);

%Representación del rotacional
surf(X,Y,ROTACIONAL)

colorbar;
view(2)
title('Rotacional')

%Puntos con mayor rotacional
Maximo_Rotacional=max(max(ROTACIONAL))

%El máximo del rotacional corresponde a puntos con THETA = (PI/4)
title('Rotacional')

%Puntos con mayor rotacional
Maximo_Rotacional=max(max(ROTACIONAL))

}}

== Tensor de tensiones ==
Si consideramos que estamos en un medio elástico lineal, isótropo y homogéneo; en mecánica de medios continuos podemos definir el tensor de tensiones "σ"como:

[[Archivo:TensorDeTensionesGrupo10A.png|marco|centro]]

siendo λ y μ los coeficientes de Lamé, que dependen de las propiedades elásticas de cada material, considerando en nuestro caso que λ y μ son iguales a 1 y que ε es la parte simétrica del vector gradiente u.

Para hacer este apartado se requiere utilizar bases generales, en lo cual hemos llegado a la conclusión de que las componentes de ∇ū son:

[[Archivo:GradienteDeUGrupo10A.png|marco|centro]]

Siendo ε =(∇u + ∇ut)/2, con todo ello, calculamos la matriz σ:

[[Archivo:MatrizDeTensionesGrupo10A.png|marco|centro]]

Se van a representar las tensiones en las tres direcciones del espacio. En el caso de estar utilzando coordenadas cilíndricas, se representan las tensiones en las direcciones de la base cilíndrica ḡρ, ḡθ, ḡz. Para hallar las tensiones en estas direcciones debemos multiplicar escalarmente la matriz σ por cada uno de los vectores ḡi por ambos lados: Tensión normal a la dirección ḡi = ḡi*σ*ḡi

[[Archivo:DireccionGrupo10A.png|marco|centro]]

{{matlab|codigo=
%Inicializamos las variables
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%En 2D

%Dibujo de las tensiones respecto a g sub rho
subplot(2,3,1);

%Tomamos el valor (1,1) de la matriz sigma y lo dibujamos
r=2.*THETA./5+THETA.*(2-1./RHO)./5;

surf(X,Y,r)

axis equal
view(2)
title('Tensiones respecto a gRho')
colorbar

%Dibujo de las tensiones respecto a g sub tetha
subplot(2,3,2);

%Tomamos el valor (2,2) de la matriz sigma y lo dibujamos
t=2.*THETA.*(RHO-1)./(5.*RHO.^5)+THETA.*(2-1./RHO)./(5.*RHO.^2);

surf(X,Y,t)

axis equal
title('Tensiones respecto a gTheta')
view(2)

colorbar

%Dibujo de las tensiones respecto a g sub z
subplot(2,3,3)

%Tomamos el valor (3,3) de la matriz sigma y lo dibujamos
z=THETA.*(2-1./RHO)./5;

surf(X,Y,z)

axis equal
view(2)
title('Tensiones respecto a gz')
colorbar

%Cálculo en 3D:

%Dibujo de las tensiones respecto a g sub rho 3D
subplot(2,3,4);

%Tomamos el valor (1,1) de la matriz sigma y lo dibujamos

surf(X,Y,r)

axis equal
title('Tensiones respecto a gRho 3D')

%Dibujo de las tensiones respecto a g sub tetha 3D
subplot(2,3,5);

%Tomamos el valor (2,2) de la matriz sigma y lo dibujamos

surf(X,Y,t)

axis equal
title('Tensiones respecto a gTheta 3D')

%Dibujo de las tensiones respecto a g sub z 3D
subplot(2,3,6)

%Tomamos el valor (3,3) de la matriz sigma y lo dibujamos

surf(X,Y,z)

axis equal
title('Tensiones respecto a gz 3D')
}}

[[Archivo:tensionesNormalesgrupo10A.png|centro|thumb|1300px]]

== Tensiones tangenciales ==

Ambas tensiones tangenciales son iguales, pues la matriz σ es simétrica.
Las tensiones del apartado nos dan: σ(2,1)= σ(2,1)= 1/ρ*(1/5*sen(2θ)-1)

== Tensiones tangenciales respecto al plano ortogonal a ḡρ ==

[[Archivo:Tensionestangencialesgrupo10A.png||thumb|390px|derecha]]

{{matlab|codigo=
%Inicializamos las variables
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos las tensiones tangentes al plano g sub rho
Tang=(RHO-1)./(5.*RHO.^2);

%Visualización en 3D
subplot(2,1,1)

surf(X,Y,Tang)

colorbar
title('Tensiones tangenciales respecto al plano ortogonal a gRho en 3D')

%Visualización en 2D
subplot(2,1,2)

surf(X,Y,Tang)

colorbar
view(2)
title('Tensiones tangenciales respecto al plano ortogonal a gRho en 2D')
}}

== Tensión de Von Mises ==

La tensión de Von Mises se define mediante la expresión:

[[Archivo:FormulaVonMises.png|marco|centro]]

La tensión de Von Mises es una magnitud escalar que se suele utilizar como indicador para saber cuando un material inicia un comportamiento plástico y no elástico puro. En nuestro caso, el valor máximo de la Tensión de Von Mises es igual a 2,0783.

[[Archivo:tensionVonMisesgrupo10A.png|marco|derecha]]

{{matlab|codigo=
%Inicializamos las variables
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos la parametrización
[RHO,THETA]=meshgrid(rho,theta);

%Creamos la matriz sigma
Sigma=zeros(3,3);

%Creamos la matriz de Von Mises
MatrizVonMises=zeros(16,21);

%Obtenemos los valores de la tensión de Von Mises a partir de los autovalores
for i=1:16*21

r=RHO(i)';

t=THETA(i)';

Sigma(1,1)=2*t/5+t*(2-1/r)/5;

Sigma(1,2)=(r-1)/(5*r^2);

Sigma(2,1)= (r-1)/(5*r^2);

Sigma(2,2)=2*t*(r-1)/(5*r^3)+t*(2-1/r)/(5);
Sigma(3,3)=t*(2-1/r)/5;

%Obtención de los autovalores de la matriz sigma
[v,d]=eig(Sigma);

%Fórmula de Von Mises
MatrizVonMises(i)=sqrt(((d(1,1)-d(2,2))^2+(d(2,2)-d(3,3))^2+(d(3,3)-d(1,1))^2)/2);

end

%Cambiamos a coordenadas cartesianas
X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos la tensión de Von Mises
surf(X,Y,MatrizVonMises)

colorbar
view(2)
title('Tensión de Von Mises')

%Valor Máximo
[i,j]=find(MatrizVonMises==max(max((MatrizVonMises))))
z=max(max(MatrizVonMises))

%El valor máximo de la Tensión de Von Mises es igual a 0.6283
}}

Visualización de campos escalares y vectoriales en elasticidad. Grupo 10-A

2016-12-01T11:48:41Z

Alejandro de Benito Andrés:

{{ TrabajoED | Visualizaciones de campos escalares y vectoriales. Grupo 10-A | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC16/17|2016-17]] | Alejandro de Benito, Lucía Rico , Adrián Piñeiro, Daniel Apellániz, Guillermo Alfaya. }}

EDITAR Consideramos una placa plana que ocupa el anillo circular centrado en el origen que ocupa el cuarto
anillo circular centrado en el origen y comprendido entre los radios 1 y 3, en el primer cuadrante
x, y ≥ 0. Dibujar con los ejes en el cuadrado [−1, 4]×[−1, 4]. Trabajaremos en coordenadas cilíndricas, en la cual definimos la temperatura T(x,y)=(y+2)·x2, dependiente de las variables (x,y), y el campo de desplazamientos u(ρ,θ)=θf(ρ)gρ. producido por la acción de una fuerza. De ellos, estudiaremos el efecto de la temperatura dada por \(T(x,y)\) y el estudio del movimiento de la placa con el campo \(\vec u(\rho,\theta)\), con las consiguientes tensiones que producirá.

'''Dada una superficie plana en forma de anillo circular, concretamente la cuarta parte situada en el primer cuadrante, con centro en el origen de radio exterior 3 e interior 1. Se han utilizado como ejes el cuadrado [-1,4]x[-1,4] para una mejor visualización de la placa. Por comodidad han sido elegidas las coordenadas cilíndricas, a la que se aplicará una temperatura T(x,y)=(y+2)·x2, en la que se toman como variables independientes "x, y" y el campo de desplazamientos u(ρ,θ)=θf(ρ)gρ, conclusión de un agente exterior desconocido. En apartados posteriores veremos como afectan la temperatura o el campo de desplazamientos a la placa en cuestión con las tensiones ejercidas.
'''
== Mallado ==
Representaremos una malla de los puntos que componen la placa sólida con forma de corona semicircular.
El paso de muestreo dado es 0,1 tanto para la variable X como para la variable Y.
[[Archivo:Malladogrupo10A.png|thumb|400px|derecha]]

{{matlab|codigo=
%Definimos las variables con el paso de muestreo "h"
h=0.1
rho=[1:h:3];
theta=[0:h:pi/2];

%Creamos el cojunto de puntos que forman la malla y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis([-1,4,-1,4]);

view(2);
title('Mallado')
}}

== Temperatura ==

La temperatura de la corona semicircular cambia según la función T(x,y)=(y+2)·x2: La función de únicamente depende de la variable Y y se mantendrá constante para todo valor de X.
Se ha dibujado el campo de temperaturas tanto en 2D como en 3D.
El punto que ha resultado tener mayor temperatura T(x,y)=24.1902 es el punto x=1.1683,y=2.7632.

'''En la corona la temperatura varía según la función dada, "T(x,y)=(y+2)·x2". Para tener una idea más realista de la temperatura la hemos dibujado tanto como en 2D como en 3D.El punto que ha resultado tener mayor temperatura T(x,y)=24.1902 es el punto x=1.1683,y=2.7632.

'''
[[File:temperaturaGrupo10A2.png|thumb|800px|derecha]]

{{matlab|codigo=
%Definimos las variables con el paso de muestreo "h"
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H); %Así se ajusta mejor a theta=pi/2

%Creamos el cojunto de puntos que forman la malla y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis equal;

view(2);

%CAMPO ESCALAR DE TEMPERATURAS
Temperatura=(Y+2).*X.^2
Maximo=max(max(Temperatura))

%Líneas de programa usadas para obtener posición del máximo
for i=1:16
for j=1:21
Temp=(Y(i,j)+2).*X(i,j).^2;
if Temp==Maximo
disp(X(i,j))
disp(Y(i,j))
end
end
end

%Representación de la malla en 2D en la subventana primera
subplot(1,2,1)

contour(X,Y,Temperatura,20);

view(2)
axis([-1,4,-1,4]);
title('Curvas de nivel de la Temperatura')

%Representación de la malla en 3D en la subventana segunda
subplot(1,2,2);

surf(X,Y,Temperatura);

view(2)
axis([-1,4,-1,4]);
title('Temperatura en 2D')
}}

== Gradiente de la temperatura ==

El gradiente del campo de temperaturas (∇T) muestra la variación de la temperatura al movernos por el plano, siendo el módulo del gradiente lo que nos indica la rapidez con la que aumenta la temperatura.
El gradiente apunta en la dirección en la cual la temperatura aumenta más rápido: donde la derivada direccional será máxima.
El vector gradiente es perpendicular a la superficie del sólido.
Cuanto más aumenta la Y, mayor es la temperatura: En el dibujo podemos observar el aumento de la temperatura a medida que los puntos del sólido avanzan en el eje Y.

'''Para conocer la variación de temperaturaamovernos por la placa usaremos el gradiente del campo de temperaturas(∇T), el módulo de los vectores gradiente representa la rapidez con la que varía la temperatura, es decir apunta hacia la dirección en la que la temperatura tiene un aumento más notable, donde la derivada direccional es máxima'''

[[Archivo:gradienteDeLaTemperaturaGrupo10A1.png|thumb|800px|Izquierda: Gráfica de la temperatura en 2D y vectores gradientes. Derecha: Vectores del gradiente|derecha]]
{{matlab|codigo=
%Inicializamos las variables
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis([-1,4,-1,4]);
axis square

view(2);

%GRADIENTE
GradX=2*X.*Y+4*X;
GradY=X.^2;

%CAMPO ESCALAR DE TEMPERATURAS
Temperatura=(Y+2).*X.^2;

%Creación de la subventana primera: gradiente en el campo de temperaturas

subplot(1,2,1)
hold on
surf(X,Y,Temperatura);
axis([-1,4,-1,4]);
axis square

%Dibujo del gradiente
quiver3(X,Y,Temperatura,GradX,GradY,0*Temperatura);
axis([-1,4,-1,4]);
axis square

hold off

%Creación de la segunda subventana: campo vectorial del gradiente

subplot(1,2,2)
quiver(X,Y,GradX,GradY);
axis([-1,4,-1,4]);
axis equal
colorbar
}}

== Campo de desplazamientos ==

El campo de desplazamientos dado es: u(ρ,θ)=θf(ρ)gρ. Vamos a dibujar el campo de vectorial, sabiendo que: Los puntos situados en ρ=1 no sufren desplazamiento y div(u)=θ(2-1/ρ)/5 ,
lo que nos da la ecuación diferencial f'(ρ)+(1/ρ)*f(ρ)=(2-1/ρ)/5 con la condición de contorno f(1)=0. Resultando f(ρ)=(ρ-1)/5.
Lo metemos en el campo de desplazamiento dando u(ρ,θ)=θ*((ρ-1)/5)gρ .Cogiendo el campo en función de los dos componente i,j. Definiremos pues FX y FY. Sabiendo que gρ= cosθi+senθj.

[[Archivo:CampoDesplazamientosGrupo10A1.png|marco|derecha]]

{{matlab|codigo=
%Inicializamos las variables
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis([-1,4,-1,4]);

view(2);

%CAMPO VECTORIAL DE DESPLAZAMIENTOS
%Componente x(X,Y) del campo vectorial
FX=THETA.*((RHO-1)/5).*(cos(THETA));

%Componente y(X,Y) del campo vectorial
FY=THETA.*((RHO-1)/5).*(sin(THETA));

%Representamos el campo vectorial de desplazamientos
quiver(X,Y,FX,FY)

axis([-1,4,-1,4])
view(2)

xlabel('Eje x')
ylabel('Eje y')
title('Campo de desplazamientos sobre la placa')
}}

== Movimiento del sólido ==

Se ha representado el sólido antes y después del desplazamiento dado por el campo vectorial para así poder compararlos. Definiremos las funciones de desplazamiento en función de ambas componentes, sumando las componentes del vector de desplazamientos a las coordenadas respectivas de posición (Las de la posición incial).
Se puede ver que el desplazamiento es mínimo.

[[Archivo:movimientoDelSolidoGrupo10A1.png|thumb|800px|derecha]]

{{matlab|codigo=
%Inicializamos las variables
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis([-1,4,-1,4]);

view(2);

%CAMPO VECTORIAL DE DESPLAZAMIENTOS
%Componente x(X,Y) del campo vectorial
FX=THETA.*((RHO-1)/5).*(cos(THETA));

%Componente y(X,Y) del campo vectorial
FY=THETA.*((RHO-1)/5).*(sin(THETA));

%Dibujo de la placa sin movimiento
subplot(1,2,1)

surf(X,Y,0*X);

axis([-1,4,-1,4]);
view(2);
title('Reposo')

%Dibujo de la placa en movimiento
%Movimiento de la placa en el eje X
MovimientoX=X+FX;

%Movimiento de la placa en el eje Y
MovimientoY=Y+FY;

subplot(1,2,2)

surf(MovimientoX,MovimientoY,Y.*0)

axis([-1,4,-1,4]);
view(2);
title('Aplicado el campo de movimientos')
}}

== Divergencia ==

La ∇·(u)=θ(2-1/ρ)/5 dada en el apartado anterior. Para calcular donde tiene los máximos,mínimos y nulos procederemos calculando su gradiente.
<big>∇(∇·u)=θ/(5*ρ²)gρ+(2-1/ρ)/5*(1/ρ²)gθ+0gz</big>

Como nos movemos entre ρ=(1,3) y θ=(0,pi/2) las derivada direccionales gρ y gθ son positivas, por lo que cuanto mayor sean ρ y θ mayor será su valor, alcanzándose el Máximo de ∇·u(ρ,θ)=0.5236 en (3,pi/2). Como ∇·(u) esta en función de θ, cuando θ=0 ∇·(u)=0 por lo que sera nulo y a su vez se alcanza su mínimo.

[[Archivo:DivergenciaGrupo10A.png|marco|derecha]]

{{matlab|codigo=
%Inicializamos las variables
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis([-1,4,-1,4]);

view(2);

%DIVERGENCIA
DIVERGENCIA=THETA.*((2-1./RHO)./5);

MaxDivergencia=max(max(DIVERGENCIA))

%Representación de la divergencia
surf(X,Y,DIVERGENCIA)

view(2)
colorbar
title('Divergencia')
}}

== Rotacional ==

Se ha calculado el rotacional del campo de desplazamientos u en todos los puntos del sólido y lo representamos. Hemos obtenido que el rotacional es: Rot [u(rho,thetha)] = (rho-1)/(5*rho)gz
Observando la gráfica se obtiene cuáles son los puntos que sufren un mayor rotacional: Los puntos con mayor rotacional son aquellos con rho=3 y su valor es 0.1333
[[Archivo:RotacionalGrupo10A.png|marco|derecha]]

{{matlab|codigo=
%Inicializamos las variables
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos el mallado
mesh(X,Y,0*X)
axis([-1,4,-1,4]);

view(2);

%ROTACIONAL
ROTACIONAL=(RHO-1)./(5.*RHO);

%Representación del rotacional
surf(X,Y,ROTACIONAL)

colorbar;
view(2)
title('Rotacional')

%Puntos con mayor rotacional
Maximo_Rotacional=max(max(ROTACIONAL))

%El máximo del rotacional corresponde a puntos con THETA = (PI/4)
title('Rotacional')

%Puntos con mayor rotacional
Maximo_Rotacional=max(max(ROTACIONAL))

}}

== Tensor de tensiones ==
Si consideramos que estamos en un medio elástico lineal, isótropo y homogéneo; en mecánica de medios continuos podemos definir el tensor de tensiones "σ"como:

[[Archivo:TensorDeTensionesGrupo10A.png|marco|centro]]

siendo λ y μ los coeficientes de Lamé, que dependen de las propiedades elásticas de cada material, considerando en nuestro caso que λ y μ son iguales a 1 y que ε es la parte simétrica del vector gradiente u.

Para hacer este apartado se requiere utilizar bases generales, en lo cual hemos llegado a la conclusión de que las componentes de ∇ū son:

[[Archivo:GradienteDeUGrupo10A.png|marco|centro]]

Siendo ε =(∇u + ∇ut)/2, con todo ello, calculamos la matriz σ:

[[Archivo:MatrizDeTensionesGrupo10A.png|marco|centro]]

Se van a representar las tensiones en las tres direcciones del espacio. En el caso de estar utilzando coordenadas cilíndricas, se representan las tensiones en las direcciones de la base cilíndrica ḡρ, ḡθ, ḡz. Para hallar las tensiones en estas direcciones debemos multiplicar escalarmente la matriz σ por cada uno de los vectores ḡi por ambos lados: Tensión normal a la dirección ḡi = ḡi*σ*ḡi

[[Archivo:DireccionGrupo10A.png|marco|centro]]

{{matlab|codigo=
%Inicializamos las variables
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%En 2D

%Dibujo de las tensiones respecto a g sub rho
subplot(2,3,1);

%Tomamos el valor (1,1) de la matriz sigma y lo dibujamos
r=2.*THETA./5+THETA.*(2-1./RHO)./5;

surf(X,Y,r)

axis equal
view(2)
title('Tensiones respecto a gRho')
colorbar

%Dibujo de las tensiones respecto a g sub tetha
subplot(2,3,2);

%Tomamos el valor (2,2) de la matriz sigma y lo dibujamos
t=2.*THETA.*(RHO-1)./(5.*RHO.^5)+THETA.*(2-1./RHO)./(5.*RHO.^2);

surf(X,Y,t)

axis equal
title('Tensiones respecto a gTheta')
view(2)

colorbar

%Dibujo de las tensiones respecto a g sub z
subplot(2,3,3)

%Tomamos el valor (3,3) de la matriz sigma y lo dibujamos
z=THETA.*(2-1./RHO)./5;

surf(X,Y,z)

axis equal
view(2)
title('Tensiones respecto a gz')
colorbar

%Cálculo en 3D:

%Dibujo de las tensiones respecto a g sub rho 3D
subplot(2,3,4);

%Tomamos el valor (1,1) de la matriz sigma y lo dibujamos

surf(X,Y,r)

axis equal
title('Tensiones respecto a gRho 3D')

%Dibujo de las tensiones respecto a g sub tetha 3D
subplot(2,3,5);

%Tomamos el valor (2,2) de la matriz sigma y lo dibujamos

surf(X,Y,t)

axis equal
title('Tensiones respecto a gTheta 3D')

%Dibujo de las tensiones respecto a g sub z 3D
subplot(2,3,6)

%Tomamos el valor (3,3) de la matriz sigma y lo dibujamos

surf(X,Y,z)

axis equal
title('Tensiones respecto a gz 3D')
}}

[[Archivo:tensionesNormalesgrupo10A.png|centro|thumb|1300px]]

== Tensiones tangenciales ==

Ambas tensiones tangenciales son iguales, pues la matriz σ es simétrica.
Las tensiones del apartado nos dan: σ(2,1)= σ(2,1)= 1/ρ*(1/5*sen(2θ)-1)

== Tensiones tangenciales respecto al plano ortogonal a ḡρ ==

[[Archivo:Tensionestangencialesgrupo10A.png||thumb|390px|derecha]]

{{matlab|codigo=
%Inicializamos las variables
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos el mallado y la parametrización
[RHO,THETA]=meshgrid(rho,theta);

X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos las tensiones tangentes al plano g sub rho
Tang=(RHO-1)./(5.*RHO.^2);

%Visualización en 3D
subplot(2,1,1)

surf(X,Y,Tang)

colorbar
title('Tensiones tangenciales respecto al plano ortogonal a gRho en 3D')

%Visualización en 2D
subplot(2,1,2)

surf(X,Y,Tang)

colorbar
view(2)
title('Tensiones tangenciales respecto al plano ortogonal a gRho en 2D')
}}

== Tensión de Von Mises ==

La tensión de Von Mises se define mediante la expresión:

[[Archivo:FormulaVonMises.png|marco|centro]]

La tensión de Von Mises es una magnitud escalar que se suele utilizar como indicador para saber cuando un material inicia un comportamiento plástico y no elástico puro. En nuestro caso, el valor máximo de la Tensión de Von Mises es igual a 2,0783.

[[Archivo:tensionVonMisesgrupo10A.png|marco|derecha]]

{{matlab|codigo=
%Inicializamos las variables
h=0.1;
H=round(pi/2/h);

rho=[1:h:3];
theta=linspace(0,pi/2,H);

%Creamos la parametrización
[RHO,THETA]=meshgrid(rho,theta);

%Creamos la matriz sigma
Sigma=zeros(3,3);

%Creamos la matriz de Von Mises
MatrizVonMises=zeros(16,21);

%Obtenemos los valores de la tensión de Von Mises a partir de los autovalores
for i=1:16*21

r=RHO(i)';

t=THETA(i)';

Sigma(1,1)=2*t/5+t*(2-1/r)/5;

Sigma(1,2)=(r-1)/(5*r^2);

Sigma(2,1)= (r-1)/(5*r^2);

Sigma(2,2)=2*t*(r-1)/(5*r^3)+t*(2-1/r)/(5);
Sigma(3,3)=t*(2-1/r)/5;

%Obtención de los autovalores de la matriz sigma
[v,d]=eig(Sigma);

%Fórmula de Von Mises
MatrizVonMises(i)=sqrt(((d(1,1)-d(2,2))^2+(d(2,2)-d(3,3))^2+(d(3,3)-d(1,1))^2)/2);

end

%Cambiamos a coordenadas cartesianas
X=RHO.*cos(THETA);
Y=RHO.*sin(THETA);

%Dibujamos la tensión de Von Mises
surf(X,Y,MatrizVonMises)

colorbar
view(2)
title('Tensión de Von Mises')

%Valor Máximo
[i,j]=find(MatrizVonMises==max(max((MatrizVonMises))))
z=max(max(MatrizVonMises))

%El valor máximo de la Tensión de Von Mises es igual a 0.6283
}}