MateWiki - Contribuciones del usuario [es]

Modelos Epidemológicos Grupo 4-C

2015-03-11T13:51:35Z

FernandoSancha:

{{ TrabajoED | Modelos Epidemológicos. Grupo 20-C | [[:Categoría:Ecuaciones Diferenciales|Ecuaciones Diferenciales]]|[[:Categoría:ED14/15|Curso 2014-15]] | Angela Béjar Gómez, Eduardo Bonet García, Gonzalo Ubeda, Elisa Pamplona Frey, Alberto Rojas, Fernando Sancha Domínguez }}
En el desarrollo de una epidemia se distingue dos tipos de individuos: Los que ya han contraido de la enfermedad I, y los que son susceptibles de contraerla por encontrarse en zona de riesgo S.
Consideramos las variables: t tiempo, S(t) población de individuos susceptibles a contraer la enfermedad, I(t) población de individuos infectados; y el sistema::
<math> \left \{ \begin{matrix} \frac{dS}{dt}=-aSI \\
\frac{dI}{dt}=aSI-bI-cI\end{matrix}\right. </math>
Suponemos que :
# La población de personas infectadas se altera por el fallecimiento o la cura de las mismas. En ambos casos, la tasa de cambio depende del numero de personas infectadas;
# La tasa de individuos que pasan de ser susceptibles a contraer la enfermedad a estar infectados es proporcional a la interacción entre el numero de individuos en ambas clases

Los parámetros a, b y c como:
#"a":Proporción de la interacción entre personas susceptibles a ser infectados y personas ya infectadas.
#"b":Proporción de personas que superan la enfermedad.
#"c":Proporción de personas que mueren a causa de la enfermedad.
== Método de Euler y Trapecio==
Consideramos la ecuación:
<math>\frac{dI}{dt}=aSI-bI-cI</math>
Con '''S=0''', es decir, <math>\frac{dI}{dt}=-bI-cI</math>. Conocidos los valores b=0.3y c=0.01, nos queda el problema de valor inicial siguiente:
<math>\left \{ \begin{matrix}\frac{dI}{dt}=-bI-cI; \\ I_{0}=2000 \end{matrix}\right.</math>: <math>b=0.3 \\ c=0.01</math>
Aplicamos el método de Euler para resolver la ecuación diferencial utilizamos un intervalo de tiempo de 0 a 50 porque a los 41 días el mínimo de infectados se reduce a cero. La gráfica muestra que el numero de infectados se reduce exponencialmente con el tiempo, en concreto,(ayudándonos de un bucle while de matlab) se observa que el número de infectados se reduce a la cuarta parte en 4,5 días.
El código matlab utilizado es:
[[Archivo:fotofoto2.png|400px|thumb|right|Comparación del método de Euler y del Trapecio]]
{{matlab|codigo=
t0=0;
tN=50;
% consideramos el intervalo de tiempo [0,50] ya que el numero
%de infectados será cero a los 41 días
I0=2000;
S=0
h=0.1;
N=round((tN-t0)/h);
t=t0:h:tN;
I=zeros(1,N+1); %euler
I(1)=I0;
Z=zeros(1,N+1); % Trapecio
Z(1)=I0;
G=zeros(1,N+1);
% Vector para cáculo del tiempo que tarda en reducirse I(i)a la cuarta parte
G(1)=I0;
a=0.003;
b=0.3;
c=0.01;
for i=1:N
I(i+1)=I(i)+h*(a*S.*I(i)-b.*I(i)-c.*I(i)); %euler
Z(i+1)=(Z(i).*(1+h/2*(-b-c)))/(1+h/2*(b+c)); %trapecio
end
i=1;
tg0=0;
tg(1)=tg0;
while G(i)>(I0/4)
G(i+1)=G(i)+h*(a*S.*G(i)-b.*G(i)-c.*G(i));
%bucle que calcula el tiempo que tarda en reducirse I(i) a la cuarta parte
tg(i+1)=tg(i)+h;
i=i+1;
end
%sacamos tabla de resultados
[t',I',Z']
hold on
plot(t,I)
plot(t,Z,'r')
hold off
legend('Euler','Trapecio','Location','best');
disp('Tiempo final:')
disp(tg(end))
}}

Consideramos ahora '''S=100''' constante a lo largo del tiempo y a=0,003, hacemos uso de un código de Matlab, como el del apartado anterior para resolver esta ecuación y se interpreta en la gráfica que cuanto mayor sea la población susceptible a contraer la enfermedad, mayor será el tiempo que tarde en erradicarse. Calculamos a su vez cuanto tarda en reducirse el numero de infectados a la cuarta parte y resulta 138.6000 días.

[[Archivo:Graficazorra.png|400px|thumb|right|Comparación del método de Euler y del Trapecio]]
{{matlab|codigo=

t0=0;
tN=600; % consideramos el intervalo de tiempo [0,600]
I0=2000;
S=100
h=0.1;
N=round((tN-t0)/h);
t=t0:h:tN;
I=zeros(1,N+1); %euler
I(1)=I0;
Z=zeros(1,N+1); % Trapecio
Z(1)=I0;
G=zeros(1,N+1);
% Vector para cáculo del tiempo que tarda en reducirse I(i)a la cuarta parte
G(1)=I0;
a=0.003;
b=0.3;
c=0.01;
for i=1:N
I(i+1)=I(i)+h*(a*S.*I(i)-b.*I(i)-c.*I(i)); %euler
Z(i+1)=(Z(i).*(1+h/2*(a*S-b-c)))/(1-h/2*(a*S-b-c)); %trapecio
end
i=1;
tg0=0;
tg(1)=tg0;
while G(i)>500
G(i+1)=G(i)+h*(a*S.*G(i)-b.*G(i)-c.*G(i));
%bucle que calcula el tiempo que tarda en reducirse I(i)a la cuarta parte
tg(i+1)=tg(i)+h;
i=i+1;
end
%sacamos tabla de resultados
[t',I',Z']
hold on
plot(t,I)
plot(t,Z,'r')
hold off
legend('Euler','Trapecio','Location','best');
disp('Tiempo final:')
disp(tg(end))
}}

Repetimos el proceso para '''S=200''' constante a lo largo del tiempo y vemos en la gráfica que la ecuación ya no es valida para este valor de S debido a que sobrepasa el límite y la enfermedad no podrá ser erradicada, ya que la ecuación crece exponencialmente con el tiempo, es decir que el numero de infectados no para de crecer. Por tanto, no tiene sentido calcular el tiempo que tarda en reducirse el numero de infectados a la cuarta parte.
[[Archivo:Graficazorra2.png|400px|thumb|right|Comparación del método de Euler y del Trapecio]]
{{matlab|codigo=
t0=0;
tN=600; % consideramos el intervalo de tiempo [0,600]
I0=2000;
S=200
h=0.1;
N=round((tN-t0)/h);
t=t0:h:tN;
I=zeros(1,N+1); %euler
I(1)=I0;
Z=zeros(1,N+1); % Trapecio
Z(1)=I0;
G=zeros(1,N+1);
% Vector para cáculo del tiempo que tarda en reducirse I(i)
%a la cuarta parte
G(1)=I0;
a=0.003;
b=0.3;
c=0.01;
for i=1:N
I(i+1)=I(i)+h*(a*S.*I(i)-b.*I(i)-c.*I(i)); %euler
Z(i+1)=(Z(i).*(1+h/2*(a*S-b-c)))/(1-h/2*(a*S-b-c)); %trapecio
end
i=1;
tg0=0;
tg(1)=tg0;
while G(i)>500
G(i+1)=G(i)+h*(a*S.*G(i)-b.*G(i)-c.*G(i));
%bucle que calcula el tiempo que tarda en reducirse I(i)
%a la cuarta parte
tg(i+1)=tg(i)+h;
i=i+1;
end
%sacamos tabla de resultados
[t',I',Z']
hold on
plot(t,I)
plot(t,Z,'r')
hold off
legend('Euler','Trapecio','Location','best');
disp('Tiempo final:')
disp(tg(end))
}}
== Modelo Completo==
En este caso obtendremos valores de S e I en función del tiempo por el método de Euler resolviendo el sistema
debido a que el sistema no es lineal aplicamos Euler para cada una de las ecuaciones, lo que nos dará el
resultado gráficamente de S e I y sus evoluciones en el tiempo, además utilizaremos cuatro longitudes de paso
lo cual nos reportará cuatro gráficas en cada caso. Lo resolvemos para unos valores iniciales I(0)=20 S(0)=800
e I(0)=40 S(0)=10000.
En el primer caso utilizamos el siguiente código MatLab:
[[Archivo:graficafer1.png|400px|thumb|right|Modelo Completo]]

{{matlab|codigo=

% Euler
clear all
%DATOS DEL PROBLEMA
t0=0;
tN=40;
I0=20;
S0=800;
h1=0.1;
h2=0.01;
h3=0.001;
h4=0.0001;
a=0.003;
b=0.3;
c=0.01;
%Calculamos número de subintervalos
N1=round((tN-t0)/h1);
N2=round((tN-t0)/h2);
N3=round((tN-t0)/h3);
N4=round((tN-t0)/h4);
% Definimos la variable independiente
t1=t0:h1:tN;
t2=t0:h2:tN;
t3=t0:h3:tN;
t4=t0:h4:tN;
%Ahora vamos a guardar los valores de la solución aproximada en el vector y
I=zeros(1,N1+1);
S=zeros(1,N1+1);
F=zeros(1,N2+1);
R=zeros(1,N2+1);
A=zeros(1,N3+1);
L=zeros(1,N3+1);
G=zeros(1,N4+1);
U=zeros(1,N4+1);
I(1)=I0;
S(1)=S0;
F(1)=I0;
R(1)=S0;
A(1)=I0;
L(1)=S0;
G(1)=I0;
U(1)=S0;
for i=1:N1
S(i+1)=S(i)+h1.*(-a.*S(i).*I(i)); %eulerS
I(i+1)=I(i)+h1.*(a.*S(i).*I(i)-b.*I(i)-c.*I(i)) ; %eulerI
end
for i=1:N2
R(i+1)=R(i)+h2.*(-a.*R(i).*F(i)); %eulerS
F(i+1)=F(i)+h2.*(a.*R(i).*F(i)-b.*F(i)-c.*F(i)) ; %eulerI
end
for i=1:N3
L(i+1)=L(i)+h3.*(-a.*L(i).*A(i)); %eulerS
A(i+1)=A(i)+h3.*(a.*L(i).*A(i)-b.*A(i)-c.*A(i)) ; %eulerI
end
for i=1:N4
U(i+1)=U(i)+h4.*(-a.*U(i).*G(i)); %eulerS
G(i+1)=G(i)+h4.*(a.*U(i).*G(i)-b.*G(i)-c.*G(i)) ; %eulerI
end
J1=max(I);
J2=max(F);
J3=max(A);
J4=max(G);
disp('El valor maximo en la grafica 1 es')
J1
disp('El valor maximo en la grafica 2 es')
J2
disp('El valor maximo en la grafica 3 es')
J3
disp('El valor maximo en la grafica 4 es')
J4
subplot(2,2,1)
hold on
plot(t1,I);
plot(t1,S,'r');
legend('INFECTADOS','SUSCEPTIBLES')
hold off
subplot(2,2,2)
hold on
plot(t2,F);
plot(t2,R,'g');
legend('INFECTADOS','SUSCEPTIBLES')
hold off
subplot(2,2,3)
hold on
plot(t3,A);
plot(t3,L,'y');
legend('INFECTADOS','SUSCEPTIBLES')
hold off
subplot(2,2,4)
hold on
plot(t4,G);
plot(t4,U,'c');
legend('INFECTADOS','SUSCEPTIBLES')
hold off }}

Éste programa de MatLab nos dará además el valor máximo de población infectada en cada gráfica que son:
1.I=517.45
2.I=506.37
3.I=505.3
4.I=505.19
En el segundo caso utilizamos el código MatLab:
[[Archivo:graficafer2.png|400px|thumb|right|Modelo Completo]]
{{matlab|codigo=
% Euler
clear all
%DATOS DEL PROBLEMA
t0=0;
tN=40;
I0=40;
S0=10000;
h1=0.1;
h2=0.01;
h3=0.001;
h4=0.0001;
a=0.003;
b=0.3;
c=0.01;
%Calculamos número de subintervalos
N1=round((tN-t0)/h1); % también N=length(t)-1;
N2=round((tN-t0)/h2); % también N=length(t)-1;
N3=round((tN-t0)/h3); % también N=length(t)-1;
N4=round((tN-t0)/h4); % también N=length(t)-1;
% Definimos la variable independiente
t1=t0:h1:tN;
t2=t0:h2:tN;
t3=t0:h3:tN;
t4=t0:h4:tN;
%Ahora vamos a guardar los valores de la solución aproximada en el vector y
I=zeros(1,N1+1);
S=zeros(1,N1+1);
F=zeros(1,N2+1);
R=zeros(1,N2+1);
A=zeros(1,N3+1);
L=zeros(1,N3+1);
G=zeros(1,N4+1);
U=zeros(1,N4+1);
I(1)=I0;
S(1)=S0;
F(1)=I0;
R(1)=S0;
A(1)=I0;
L(1)=S0;
G(1)=I0;
U(1)=S0;
for i=1:N1
S(i+1)=S(i)+h1.*(-a.*S(i).*I(i)); %eulerS
I(i+1)=I(i)+h1.*(a.*S(i).*I(i)-b.*I(i)-c.*I(i)) ; %eulerI
end
for i=1:N2
R(i+1)=R(i)+h2.*(-a.*R(i).*F(i)); %eulerS
F(i+1)=F(i)+h2.*(a.*R(i).*F(i)-b.*F(i)-c.*F(i)) ; %eulerI
end
for i=1:N3
L(i+1)=L(i)+h3.*(-a.*L(i).*A(i)); %eulerS
A(i+1)=A(i)+h3.*(a.*L(i).*A(i)-b.*A(i)-c.*A(i)) ; %eulerI
end
for i=1:N4
U(i+1)=U(i)+h4.*(-a.*U(i).*G(i)); %eulerS
G(i+1)=G(i)+h4.*(a.*U(i).*G(i)-b.*G(i)-c.*G(i)) ; %eulerI
end
J1=max(I);
J2=max(F);
J3=max(A);
J4=max(G);
disp('El valor maximo en la grafica 1 es')
J1
disp('El valor maximo en la grafica 2 es')
J2
disp('El valor maximo en la grafica 3 es')
J3
disp('El valor maximo en la grafica 4 es')
J4
subplot(2,2,1)
hold on
plot(t1,I);
plot(t1,S,'r');
legend('INFECTADOS','SUSCEPTIBLES')
hold off
subplot(2,2,2)
hold on
plot(t2,F);
plot(t2,R,'g');
legend('INFECTADOS','SUSCEPTIBLES')
hold off
subplot(2,2,3)
hold on
plot(t3,A);
plot(t3,L,'y');
legend('INFECTADOS','SUSCEPTIBLES')
hold off
subplot(2,2,4)
hold on
plot(t4,G);
plot(t4,U,'c');
legend('INFECTADOS','SUSCEPTIBLES')
hold off }}

Éste programa de MatLab nos dará el valor máximo de población infectada en cada gráfica que son:
1.I=1.26+e004
2.I=9521.1
3.I=9469.6
4.I=9464.7
Observamos en este segundo programa que las gráficas y resultados no parecen concordar con lo que sería la realidad
ésto puede ser debido a que la alta diferencia entre los valores iniciales de I y S no permite elcálculo efectivo de
su variación respecto del tiempo y por tanto la infección no podrá ser erradicada.

== Método Runge-Kutta==
Haciendo uso de este método de cuarto orden para resolver la ecuación <math>\frac{dI}{dt}=-bI-cI</math>
Al comparar el resultado con el método de Euler observamos que los resultados para ambos métodos son muy próximos. No es aconsejable usar un metodo implicito como el trapezoidal debido a que es demasiado complejo despejar la ecuación analíticamente Por tanto se usa normalmente un método como el de Runge-Kutta.
Según vemos en la gráfica la ecuación crece exponencialmente con el tiempo, esto es, que el numero de infectados no para de crecer por lo que la enfermedad no podrá ser erradicada
[[Archivo:Graficazorra3.png|400px|thumb|right|Comparación del método Runge-Kutta y de Euler]]
{{matlab|codigo=
t0=0;
tN=25;
% consideramos el intervalo de tiempo [0,50]
%ya que el numero de infectados será cero a los 41 días
I0=2000;
S=0
h=0.1;
N=round((tN-t0)/h);
t=t0:h:tN;
I=zeros(1,N+1); %euler
I(1)=I0;
a=0.003;
b=0.3;
c=0.01;
r=zeros(1,N+1); %RK4
r(1)=I0;
for i=1:N
I(i+1)=I(i)+h*(a*S.*I(i)-b.*I(i)-c.*I(i)); %euler
% RK4
K1=a*S.*r(i)-b.*r(i)-c.*r(i);
K2=-b.*(r(i)+(1/2).*K1.*h)-c.*(r(i)+(1/2).*K1.*h);
K3=-b.*(r(i)+(1/2).*K2.*h)-c.*(r(i)+(1/2).*K2.*h);
K4=-b.*(r(i)+K3*h)-c.*(r(i)+K3*h);
r(i+1)= r(i)+h/6*(K1+2*K2+2*K3+K4)
end

%sacamos tabla de resultados
[t',I',r']
hold on
plot(t,I)
plot(t,r,'r')
hold off
legend('Euler','Runge Kutta','Location','best');
}}

== Método de Heun==
Suponemos ahora que <math>a\left ( t \right )=\frac{0.003}{1+t}</math>,lo que significa que la interaccion entre infectados y susceptibles disminuye a medida que aumenta el tiempo.
Tambien tomamos <math> I_{O}=40</math> y <math> S_{O}=1640</math>. Con lo que el sistema nos quedaría asi:

<math> \left \{ \begin{matrix}\frac{dS}{dt}= \frac{-0.003}{1+t}\cdot S\cdot I \\ \frac{dI}{dt}= \frac{0.003}{1+t}\cdot S\cdot I-0.3\cdot I-0.01\cdot I\end{matrix}\right.</math>
En la gráfica vemos que el número de infectados alcanza un máximo entorno a t=2, y luego decrece exponencialmente con el tiempo, sin embargo el número de susceptibles empieza en 1640 y decrece exponencialmente con el tiempo más rápido que el número de infectados

Para resolver este sistema hacemos uso del metodo de Heun, ayudandonos del siguiente código :
[[Archivo:ejercicio69.png|400px|thumb|right|Método de Heun]]
{{matlab|codigo=
%sistemas no lineales, no puedo usar matrices
clear ; clc;
t0=0;
tN=10;
h=0.0001;
N=round((tN-t0)/h); %numero de subintervalos

t=t0:h:tN;
t(1)=t0;
I=zeros(1,N+1); %INFECTADOS
S=zeros(1,N+1); %SUBCEPTIBLES
I(1)=40;
S(1)=1640;
a=zeros(1,N+1);
a(1)=0.003/(1+t0);

for i=1:N
a(i+1)=0.003/(1+t(i));
K1S=-a(i).*S(i).*I(i);
K2S=-a(i)*S(i).*I(i)+h*K1S;
S(i+1)=S(i)+(h/2).*(K1S+K2S);
K1I=a(i).*S(i)*I(i)-0.3.*I(i)-0.01*I(i);
K2I=a(i).*S(i)*I(i)-0.31*I(i)+K1I*h;
I(i+1)=I(i)+(h/2)*(K1I+K2I);
end
[S,I,t]';
hold on
plot(t,S);
plot(t,I,'r');
legend('SUSCEPTIBLES','INFECTADOS');
hold off
}}

Basándonos en los resultados obtenidos, resolvemos el sistema para diferentes valores de la constante a. Con el programa de matlab desarrollado a continuación obtenemos una tabla con todos los valores de tiempo en los que se alcanza el número máximo de infectados para cada valor de a. Concluimos que el valor de a para el cual el máximo está más cerca de 5 es 0.0005, es decir, el menor del intervalo, y a medida que crece la I máxima se alcanza antes que t=5. A continuación se muestra el programa utilizado así como la tabla obtenida con este.
[[Archivo:bonethijoputa.png|400px|thumb|right|Tabla de datos obtenidos]]

{{matlab|codigo=

clear all;
clc;
h=0.0001;
t0=0;
tN=5;
N=(tN-t0)/h; %numero de subintervalos
% Definimos la variable independiente
t=t0:h:tN;
t(1)=t0;
I=zeros(1,N+1); %INFECTADOS
S=zeros(1,N+1); %SUBCEPTIBLES
I(1)=40;
S(1)=1640;
c=[]
for a=0.0005:h:0.002
for i=1:N
K1S=a*S(i)*I(i);
K2S=a*(S(i)+K1S*h)*(I(i)+h*K1S);
K1I=a*S(i)*I(i)-0.3*I(i)-0.01*I(i);
K2I=a*(S(i)+h*K1S)*(I(i)+K1I*h)-0.3*(I(i)+K1I*h)-(0.01*(I(i)+K1I*h));
I(i+1)=I(i)+(h/2)*(K1I+K2I);
S(i+1)=S(i)+(h/2)*(K1S+K2S);
[v,w]=max(I);
end
c=[c;[a,t(w)]];
end

}}

[[Categoría:Ecuaciones Diferenciales]]

[[Categoría:ED14/15]]

[[Categoría:Trabajos 2014-15]]

Visualización de campos escalares y vectoriales en elasticidad. Grupo C4

2014-12-05T13:36:28Z

FernandoSancha: /* Masa total de la placa */

{{ TrabajoED |Visualización de campos escalares y vectoriales en elasticidad. Grupo C4 | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC14/15|2013-14]] | Alberto Rojas Rivero, Eduardo Bonet García, Fernando Sancha Domínguez, Carlos Fernández Bermejo, Ricardo Mazón Cabrera, Emilio Valero Muñoz-Rojas }}

Consideramos una placa circular comprendida entre las circunferencias de radio 1 y radio 2, centradas en el origen.
En ella vamos a suponer que tenemos definidas dos cantidades físicas: la temperatura <math>T(\rho,\theta,t)</math> ,que depende de las dos variables espaciales <math>(x,y)</math> , y el tiempo <math>t</math>, y los desplazamientos <math>\vec u(\rho,\theta)</math>.De esta forma, si definimos <math>r_0(\rho,\theta)</math> el vector de posición de los puntos de la placa antes de la deformación , la posición de cada punto <math>(\rho,\theta)</math> de la placa después de la deformación viene dada por:
<math>
\vec r (\rho,\theta)= \vec r_{0}(\rho,\theta)+\vec u(\rho,\theta).
</math>
Vamos a suponer que la fuerza aplicada sobre la placa ha provocado una desplazamiento de manera que los desplazamientos en un instante <math>t_0</math> vienen dados por:
<math>\vec u(\rho,\theta)=(1-\rho)^2 \vec g_{\theta}.
</math>

== Mallado del sólido ==
Nuestro mallado consiste en la representación de los puntos interiores de una placa solida, en forma de corona circular. Tomando como paso de muestreo h=1/10, los comandos en Matlab de la visualizacion de nuestra placa es la siguiente:
[[Archivo:Sin título.jpg|350px|miniaturadeimagen|derecha|Mallado de la corona circular.]]
{{matlab|codigo=
h=0.1 % Paso de muestreo
u=1:h:2; % Intervalo [1,2]
v=0:h:2*pi+h; % Intervalo [0,2*pi]
[uu,vv]=meshgrid(u,v); % Matrices de las coordenadas u y v
figure(1)
xx=uu.*cos(vv); % parametrizacion
yy=uu.*sin(vv);
mesh(xx,yy,0*xx) % dibuja la malla
axis([-3,3,-3,3]) % selecciona la región
view(2) }}

== Temperatura del sólido ==
La temperatura es un campo escalar representado por <math>T(x,y)=e^{-y}</math>,vamos a tratar de ver como varía la temperatura en la placa que estamos estudiando . Nuestra función depende únicamente de la variable Y, esto quiere decir que la temperatura se mantendrá constante para todo valor de X. Observamos que cuanto mayor es el valor de la Y ,menor es la temperatura , luego los puntos de máxima temperatura son los más alejados del foco( el origen) y con menor valor de la variable Y, como podemos observar en la figura de la derecha.
[[Archivo:2.jpg|325px|derecha|Distribución de la temperatura en la placa]]
{{matlab|codigo=
h=0.1 % paso de muestreo
u=1:h:2; % Intervalo[1,2]
v=0:h:2*pi+h; % Intervalo [0,2*pi]
[uu,vv]=meshgrid(u,v); % matrices
figure(1)
xx=uu.*cos(vv); % parametrización
yy=uu.*sin(vv);
f=exp(-yy); % Campo escalar
surf(xx,yy,f) % visualización
axis([-3,3,-3,3]) % región de visualización
view(2)}}

==Cálculo y representación del <math>∇T </math>==

Para estudiar como varía la temperatura al movernos por el plano calculamos el gradiente de T.
<math>∇T= -e^{-y} j </math>
En la figura de la derecha podemos observar la representación del campo vectorial <math>∇T </math> y que este es perpendicular a las curvas de nivel del campo escalar T.
[[Archivo:Anillo doblado.jpg|400px|derecha| Curvas de nivel de T y gradiente de T]]
{{matlab|codigo=
g=-exp(-yy);
subplot(1,2,1);
hold on
contour(xx,yy,f,20)
surf(xx,yy,f);
quiver3(xx,yy,f,0*xx,g,0*f);
hold off
subplot(1,2,2)
hold on
quiver(xx,yy,0*xx,Gradiente); colorbar;
contour(xx,yy,f,20)
hold off }}

==Campo de desplazamientos==
Supongamos la aplicación de una fuerza que provoca una vibración sobre la placa, de manera que, en un t aleatorio los desplazamientos vienen dados por:
<math>\vec u(\rho,\theta)=(1-\rho)^2 \vec g_{\theta}.
</math>
Para distinguir el campo de vectores en los puntos del mallado del sólido utilizaremos el código matlab siguiente:
[[Archivo:campodesplazam.jpg|500px|derecha|Distribución de la temperatura en la placCampo de desplazamientos de la placa]]
{{matlab|codigo=
u=1:h:2; % definir intervalo rho [1,2]

v=0:h:2*pi+h; % definir intervalo theta [0,2*pi]

[uu,vv]=meshgrid(u,v); % definir mallado

xx=uu.*cos(vv); % parametrizacion

yy=uu.*sin(vv);

ux=(xx.*(1-uu).^2);

uy=(yy.*(1-uu).^2);

subplot(1,2,1);

mesh(xx,yy,0.*xx)

view(2)

axis([-3,3,-3,3])

hold on

quiver(xx,yy,ux,uy)

title('campo con vectores')

subplot(1,2,2)

plot3(ux,uy,0.*xx)

title('campo de deformaciones')

}}

==Sólido antes y despues de desplazamiento==
El desplazamiento realizado por la placa viene determinado por una serie de movimientos transversales en direcciones distintas. En las imagenes adjuntas calculadas mediante Matlab, se aprecia que el desplazamiento es practicamente nulo.
[[Archivo:antesydespues.jpg|500px|derecha|Placa antes y después del desplazamiento]]
{{matlab|codigo=
h=0.1;
ro=1:h:2; % Valor de ro según el muestreo del intervalo
teta=0:h:2*pi+h; % Valor de teta según el muestreo del intervalo
[roro,tetateta]=meshgrid(ro,teta); %Matrices de la placa en coordenadas polares
figure(1)
xx=roro.*cos(tetateta); % parametrización
yy=roro.*sin(tetateta);
subplot(1,3,1) % figura solido momento inicial
mesh(xx,yy,0*xx) % mallado
axis([-4,4,-4,4]); % ejes
view(2); %visualización
desx=xx+((1-roro).^2).*(-yy); % movimiento parametrizado
desy=yy+((1-roro).^2).*(xx);
subplot(1,3,2) % figura solido en movimiento
mesh(desx,desy,yy.*0) % mallado
axis([-4,4,-4,4]);
view(2);
subplot(1,3,3); % figura solido en movimiento y campo vectorial
hold on
mesh(desx,desy,yy.*0)
axis([-4,4,-4,4]); % ejes
view(2); %
quiver(xx,yy,fx,fy) % campo vectorial del gradiente de temperaturas
axis([-4,4,-4,4]) % ejes
view(2)
hold off

}}

==Divergencia de u==
La Divergencia mide la diferencia de flujo de un campo vectorial sobre una superficie. Como vemos en la gráfica la divergencia está calculada en todos los puntos de la placa y es siempre mayor que cero, lo que nos indica que no tiene pérdidas de flujo, sino que la diferencia es cero por tanto la placa solamente tendrá fuentes.

[[Archivo:divergence.jpg|350px|derecha|Divergencia en la placa]]
{{matlab|codigo=
h=(0.1);
ro=1:h:2;
teta=0:h:2*pi+h;
[roro,tetateta]=meshgrid(ro,teta);
figure(1)
xx=roro.*cos(tetateta);
yy=roro.*sin(tetateta)
f=(1./roro)-4+3.*roro; % El campo divergencia
surf(xx,yy,f) % Dibujmos la malla
axis([-2,2,-2,2])
colorbar
title('Divergencia del campo del desplazamiento')
view(2)
}}

==Rotacional de u==
A traves del programa de Matlab, que podemos observar a continuacion, hemos calculado el rotacional del campo desplazamiento u en los puntos del solido y los hemos representado, observandose en la grafica cuales son los puntos que sufren mayor rotacional.
[[Archivo:rotacionalc4.jpg|400px|derecha|Rotacional en la placa]]
{{matlab|codigo=
h=(0.1);
ro=1:h:2;
teta=0:h:2*pi+h;
[roro,tetateta]=meshgrid(ro,teta);
figure(1)
xx=roro.*cos(tetateta);
yy=roro.*sin(tetateta);
f=(2+4.*roro.^2-6.*roro);
surf(xx,yy,f)
axis([-2,2,-2,2])
colorbar
title('Rotacional del campo del desplazamiento')
view(2)

}}

== Estudio del tensor de tensiones ==
Definamos <math>\epsilon(\vec u)=(\nabla \vec u + \nabla \vec u^t)/2</math>, la parte simétrica del tensor gradiente de <math>\vec u</math>, que se denomina tensor de deformaciones. En un medio elástico lineal, isótropo y homogéneo los desplazamientos permiten escribir el tensor de tensiones <math>\sigma_{ij}</math> a través de la fórmula::
<math>
\sigma_{ij}=\lambda \nabla \cdot \vec u \delta_{ij} + 2\mu \epsilon_{ij},
</math>
Siendo <math>\nabla \cdot \vec u = \frac{1}{\sqrt{g}}\cdot \ \frac{\partial(\sqrt{g}u^{i})}{\partial x^{i}}=\frac{1}{\rho}\cdot (\frac{\partial \rho - 2\rho ^{2}+ \rho ^{3}}{\partial \rho }+\frac{\partial 0}{\partial \theta })=\frac{1}{\rho }-4+3\rho
</math>

y calculando el gradiente de <math>\vec u</math> como:: <math>\nabla\vec u=u^i_{.j}\vec g_i \otimes \vec g^j</math> donde

<math>u^i_{.j}= \vec g^i \nabla\vec u \vec g_j=\frac{\partial u^i}{\partial x^j}+\Gamma^i_{kj}u^k.= \left(\begin{array}{ccc}-2+2\rho& 0 \\0& \ 0\end{array}\right)+ (1-2\rho +\rho ^{2}).\left(\begin{array}{ccc}0& 0 \\0& \ 1\rho\end{array}\right)
</math>

Podemos calcular el tensor gradiente como 1-contravariante 1-covariante:
:<math>\nabla\vec u=u^i_{.j}=\left(\begin{array}{ccc}-2+2\rho& 0 \\0& \ 1/\rho -2 +\rho\end{array}\right)
</math>

Las coordenadas 1-covariantes 1-contravariante coinciden por lo que la transposicion seá directa y el tensor quedara como <math>(\epsilon^i_{.j})=(\nabla \vec u + \nabla \vec u^t)/2</math> , que corrresponde con la matriz mostrada anteriormente

Sustituyendo ahora los valores en la fórmula del tensor de tensiones,y tomando \(\lambda=\mu=1\), siendo \(\lambda\) y \(\mu\) los conocidos como coeficientes de Lamé, obtenemos:

:<math>\left(\begin{array}{ccc}1/\rho-8+7\rho& 0 \\0& \ 3/\rho -8 +5\rho\end{array}\right)</math>

Las tensiones normales de ambos vectores de la base nos quedaran como:

Las tensiones normales en la dirección \(\vec g_{\rho}\) serían:\[\vec g_{\rho} \cdot \sigma \cdot \vec g_\rho =1/\rho-8+7\rho\]

Las tensiones normales en la dirección \(\vec g_\theta/\rho\) serían:\[\vec g_\theta/\rho \cdot \sigma \cdot \vec g_\theta/\rho = 3/\rho -8 +5\rho\]

A continuación se muestra el código Matlab empleado para dibujar estas tensiones normales en ambas direcciones:
[[Archivo:tensionesnormales.jpg|500px|derecha|Tensiones normales en la placa]]
{{matlab|codigo=
h=0.1; % Intervalo de separación

u=[1:0.1:2]; % Intervalo de rho [1,2]

v=[0:h:2*pi+h]; % Intervalo de theta [0,2*pi]

[uu,vv]=meshgrid(u,v); % Matrices de rho and theta

figure(1)

subplot(1,2,1)

xx=uu.*cos(vv); % Parametrización

yy=uu.*sin(vv);

f=(1./uu)-8+7.*uu; % Campo escalar

surf(xx,yy,f) % Dibujar el mallado

axis([-2,2,-2,2]) % Región seleccionada

axis equal

colorbar

title('Tension normal en dirección g sub rho')

view(2)

subplot(1,2,2)

f=(3./uu)-8+5.*uu; % Campo escalar

surf(xx,yy,f) % Dibujar el mallado

axis([-2,2,-2,2]) % Región del dibujo

axis equal

colorbar

title('Tension normal en dirección g sub theta')

view(2)

}}

Vemos como la tensiones serán mayores a medida que nos alejamos del centro de masas de nuestra placa y simétricas respecto de este ya que dependerá de la distancia de cada uno de los puntos al origen. Además se aprecia que en los bordes interiores de las placas dicha tensión será nula.

===Tensiones tangenciales===
Al realizar el calculo de las tensiones tangenciales de nuestra placa para las direcciones <math>\vec g_\rho</math> y a <math>\vec g_\theta/\rho</math>::
<math>|\sigma \cdot \vec g_\rho-(\vec g_\rho \cdot \sigma \cdot \vec g_\rho) \vec g_\rho| = 0
</math><math>|\sigma \cdot \vec g_\theta/\rho-(\vec g_\theta/\rho \cdot \sigma \cdot \vec g_\theta/\rho) \vec g_\theta/\rho|= 0 </math>

Vemos que ambas tensiones son nulas, algebraicamente hablando esto es debido a que nuestro campo de desplazamientos depende únicamente del valor de la distancia al origen. Esto implica que al estar definido dicho tensor por la divergencia y el tensor unitario, será nula la tensión tangencial siempre que el campo no dependa del angulo, ya que sera igual en todas las direcciones.

==Tensión de Von Mises==
Se define por la fórmula: <math> \sigma _{VM}=\sqrt{\frac{(\sigma 1-\sigma 2)^{2}+(\sigma 2-\sigma 3)^{2}+(\sigma 3-\sigma 1)^{2}}{2}}</math>
Esta tensión viene dada proporcionalmente a la Energía de Distorsión. En ingeniería se utiliza como un indicador de tensiones a las cuales el material falla (teoría de fallo), este cálculo se realiza sobre todo a los materiales dúctiles.
En nuestro caso el campo de tensiones de Von Misses viene dado por la fórmula <math> \sigma _{vm}=\frac{1}{\sqrt{2}}\cdot \left ( \frac{-2}{\rho }+2\cdot \rho \right ) </math>.

La representaremos mediante el siguiente programa en MATLAB:
[[Archivo:tensionesvon.jpg|500px|derecha|Tensiones Von Misses sobre la placa]]
{{matlab|codigo=
h=(0.1);
ro=1:h:2;
teta=0:h:2*pi+h;
[roro,tetateta]=meshgrid(ro,teta);
figure(1)
xx=roro.*cos(tetateta);
yy=roro.*sin(tetateta)
f=((-2./roro)-4+2.*roro)/sqrt(2); % El campo de tensiones
surf(xx,yy,f) % Dibujamos la malla
axis([-2,2,-2,2])
colorbar
title('Tensiones de Von Misses')
view(2) }}

==Masa total de la placa==
Si la densidad de la placa es <math> d(x,y,z)=xye^{-1/x^{2}} </math> y el área del anillo es <math> \pi.R^{2}-\pi.r^{2} </math> ,siendo R y r los radios exteriores e interiores respectivamente, la masa total de la placa será el producto del área y la densidad. Al tener la funcíon de la densidad, necesitamos integrarla para calcularla en todo el dominio.
Para una mejor resolución, efectuaremos el cambio a coordenadas polares:
:<math> x=\rho \cdot \cos \theta ; y=\rho \cdot \sin \theta </math>
Tendremos entonces que resolver la integral:
<math>\int_{0 }^{2\pi} \int_{1 }^{2}\rho ^{3}\cos \theta . \sin \theta .\log_{10}(\rho .\cos \theta +2)d\rho d\theta </math>
En matlab sería:
{{matlab|codigo=
h=1/500; % Número de particiones del trapecio
u=1:h:2;
v=0:h:2*pi+h;
[uu,vv]=meshgrid(ro,teta); % Malla de la placa en polares
xx=uu.*cos(vv); % Referencia en cartesianas
yy=uu.*sin(vv);
d=xx.*yy.*log(xx+2); % Función de densidad de la placa
a=h^2*d; % Masa en cada punto
mas=sum(sum(a)) % Suma de todas las masas puntuales)
}}
La masa de la placa es <math> 2,097\cdot 10^{-6} </math>

[[Categoría:Teoría de Campos]]
[[Categoría:TC14/15]]

Visualización de campos escalares y vectoriales en elasticidad. Grupo C4

2014-12-05T13:35:50Z

FernandoSancha: /* Masa total de la placa */

{{ TrabajoED |Visualización de campos escalares y vectoriales en elasticidad. Grupo C4 | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC14/15|2013-14]] | Alberto Rojas Rivero, Eduardo Bonet García, Fernando Sancha Domínguez, Carlos Fernández Bermejo, Ricardo Mazón Cabrera, Emilio Valero Muñoz-Rojas }}

Consideramos una placa circular comprendida entre las circunferencias de radio 1 y radio 2, centradas en el origen.
En ella vamos a suponer que tenemos definidas dos cantidades físicas: la temperatura <math>T(\rho,\theta,t)</math> ,que depende de las dos variables espaciales <math>(x,y)</math> , y el tiempo <math>t</math>, y los desplazamientos <math>\vec u(\rho,\theta)</math>.De esta forma, si definimos <math>r_0(\rho,\theta)</math> el vector de posición de los puntos de la placa antes de la deformación , la posición de cada punto <math>(\rho,\theta)</math> de la placa después de la deformación viene dada por:
<math>
\vec r (\rho,\theta)= \vec r_{0}(\rho,\theta)+\vec u(\rho,\theta).
</math>
Vamos a suponer que la fuerza aplicada sobre la placa ha provocado una desplazamiento de manera que los desplazamientos en un instante <math>t_0</math> vienen dados por:
<math>\vec u(\rho,\theta)=(1-\rho)^2 \vec g_{\theta}.
</math>

== Mallado del sólido ==
Nuestro mallado consiste en la representación de los puntos interiores de una placa solida, en forma de corona circular. Tomando como paso de muestreo h=1/10, los comandos en Matlab de la visualizacion de nuestra placa es la siguiente:
[[Archivo:Sin título.jpg|350px|miniaturadeimagen|derecha|Mallado de la corona circular.]]
{{matlab|codigo=
h=0.1 % Paso de muestreo
u=1:h:2; % Intervalo [1,2]
v=0:h:2*pi+h; % Intervalo [0,2*pi]
[uu,vv]=meshgrid(u,v); % Matrices de las coordenadas u y v
figure(1)
xx=uu.*cos(vv); % parametrizacion
yy=uu.*sin(vv);
mesh(xx,yy,0*xx) % dibuja la malla
axis([-3,3,-3,3]) % selecciona la región
view(2) }}

== Temperatura del sólido ==
La temperatura es un campo escalar representado por <math>T(x,y)=e^{-y}</math>,vamos a tratar de ver como varía la temperatura en la placa que estamos estudiando . Nuestra función depende únicamente de la variable Y, esto quiere decir que la temperatura se mantendrá constante para todo valor de X. Observamos que cuanto mayor es el valor de la Y ,menor es la temperatura , luego los puntos de máxima temperatura son los más alejados del foco( el origen) y con menor valor de la variable Y, como podemos observar en la figura de la derecha.
[[Archivo:2.jpg|325px|derecha|Distribución de la temperatura en la placa]]
{{matlab|codigo=
h=0.1 % paso de muestreo
u=1:h:2; % Intervalo[1,2]
v=0:h:2*pi+h; % Intervalo [0,2*pi]
[uu,vv]=meshgrid(u,v); % matrices
figure(1)
xx=uu.*cos(vv); % parametrización
yy=uu.*sin(vv);
f=exp(-yy); % Campo escalar
surf(xx,yy,f) % visualización
axis([-3,3,-3,3]) % región de visualización
view(2)}}

==Cálculo y representación del <math>∇T </math>==

Para estudiar como varía la temperatura al movernos por el plano calculamos el gradiente de T.
<math>∇T= -e^{-y} j </math>
En la figura de la derecha podemos observar la representación del campo vectorial <math>∇T </math> y que este es perpendicular a las curvas de nivel del campo escalar T.
[[Archivo:Anillo doblado.jpg|400px|derecha| Curvas de nivel de T y gradiente de T]]
{{matlab|codigo=
g=-exp(-yy);
subplot(1,2,1);
hold on
contour(xx,yy,f,20)
surf(xx,yy,f);
quiver3(xx,yy,f,0*xx,g,0*f);
hold off
subplot(1,2,2)
hold on
quiver(xx,yy,0*xx,Gradiente); colorbar;
contour(xx,yy,f,20)
hold off }}

==Campo de desplazamientos==
Supongamos la aplicación de una fuerza que provoca una vibración sobre la placa, de manera que, en un t aleatorio los desplazamientos vienen dados por:
<math>\vec u(\rho,\theta)=(1-\rho)^2 \vec g_{\theta}.
</math>
Para distinguir el campo de vectores en los puntos del mallado del sólido utilizaremos el código matlab siguiente:
[[Archivo:campodesplazam.jpg|500px|derecha|Distribución de la temperatura en la placCampo de desplazamientos de la placa]]
{{matlab|codigo=
u=1:h:2; % definir intervalo rho [1,2]

v=0:h:2*pi+h; % definir intervalo theta [0,2*pi]

[uu,vv]=meshgrid(u,v); % definir mallado

xx=uu.*cos(vv); % parametrizacion

yy=uu.*sin(vv);

ux=(xx.*(1-uu).^2);

uy=(yy.*(1-uu).^2);

subplot(1,2,1);

mesh(xx,yy,0.*xx)

view(2)

axis([-3,3,-3,3])

hold on

quiver(xx,yy,ux,uy)

title('campo con vectores')

subplot(1,2,2)

plot3(ux,uy,0.*xx)

title('campo de deformaciones')

}}

==Sólido antes y despues de desplazamiento==
El desplazamiento realizado por la placa viene determinado por una serie de movimientos transversales en direcciones distintas. En las imagenes adjuntas calculadas mediante Matlab, se aprecia que el desplazamiento es practicamente nulo.
[[Archivo:antesydespues.jpg|500px|derecha|Placa antes y después del desplazamiento]]
{{matlab|codigo=
h=0.1;
ro=1:h:2; % Valor de ro según el muestreo del intervalo
teta=0:h:2*pi+h; % Valor de teta según el muestreo del intervalo
[roro,tetateta]=meshgrid(ro,teta); %Matrices de la placa en coordenadas polares
figure(1)
xx=roro.*cos(tetateta); % parametrización
yy=roro.*sin(tetateta);
subplot(1,3,1) % figura solido momento inicial
mesh(xx,yy,0*xx) % mallado
axis([-4,4,-4,4]); % ejes
view(2); %visualización
desx=xx+((1-roro).^2).*(-yy); % movimiento parametrizado
desy=yy+((1-roro).^2).*(xx);
subplot(1,3,2) % figura solido en movimiento
mesh(desx,desy,yy.*0) % mallado
axis([-4,4,-4,4]);
view(2);
subplot(1,3,3); % figura solido en movimiento y campo vectorial
hold on
mesh(desx,desy,yy.*0)
axis([-4,4,-4,4]); % ejes
view(2); %
quiver(xx,yy,fx,fy) % campo vectorial del gradiente de temperaturas
axis([-4,4,-4,4]) % ejes
view(2)
hold off

}}

==Divergencia de u==
La Divergencia mide la diferencia de flujo de un campo vectorial sobre una superficie. Como vemos en la gráfica la divergencia está calculada en todos los puntos de la placa y es siempre mayor que cero, lo que nos indica que no tiene pérdidas de flujo, sino que la diferencia es cero por tanto la placa solamente tendrá fuentes.

[[Archivo:divergence.jpg|350px|derecha|Divergencia en la placa]]
{{matlab|codigo=
h=(0.1);
ro=1:h:2;
teta=0:h:2*pi+h;
[roro,tetateta]=meshgrid(ro,teta);
figure(1)
xx=roro.*cos(tetateta);
yy=roro.*sin(tetateta)
f=(1./roro)-4+3.*roro; % El campo divergencia
surf(xx,yy,f) % Dibujmos la malla
axis([-2,2,-2,2])
colorbar
title('Divergencia del campo del desplazamiento')
view(2)
}}

==Rotacional de u==
A traves del programa de Matlab, que podemos observar a continuacion, hemos calculado el rotacional del campo desplazamiento u en los puntos del solido y los hemos representado, observandose en la grafica cuales son los puntos que sufren mayor rotacional.
[[Archivo:rotacionalc4.jpg|400px|derecha|Rotacional en la placa]]
{{matlab|codigo=
h=(0.1);
ro=1:h:2;
teta=0:h:2*pi+h;
[roro,tetateta]=meshgrid(ro,teta);
figure(1)
xx=roro.*cos(tetateta);
yy=roro.*sin(tetateta);
f=(2+4.*roro.^2-6.*roro);
surf(xx,yy,f)
axis([-2,2,-2,2])
colorbar
title('Rotacional del campo del desplazamiento')
view(2)

}}

== Estudio del tensor de tensiones ==
Definamos <math>\epsilon(\vec u)=(\nabla \vec u + \nabla \vec u^t)/2</math>, la parte simétrica del tensor gradiente de <math>\vec u</math>, que se denomina tensor de deformaciones. En un medio elástico lineal, isótropo y homogéneo los desplazamientos permiten escribir el tensor de tensiones <math>\sigma_{ij}</math> a través de la fórmula::
<math>
\sigma_{ij}=\lambda \nabla \cdot \vec u \delta_{ij} + 2\mu \epsilon_{ij},
</math>
Siendo <math>\nabla \cdot \vec u = \frac{1}{\sqrt{g}}\cdot \ \frac{\partial(\sqrt{g}u^{i})}{\partial x^{i}}=\frac{1}{\rho}\cdot (\frac{\partial \rho - 2\rho ^{2}+ \rho ^{3}}{\partial \rho }+\frac{\partial 0}{\partial \theta })=\frac{1}{\rho }-4+3\rho
</math>

y calculando el gradiente de <math>\vec u</math> como:: <math>\nabla\vec u=u^i_{.j}\vec g_i \otimes \vec g^j</math> donde

<math>u^i_{.j}= \vec g^i \nabla\vec u \vec g_j=\frac{\partial u^i}{\partial x^j}+\Gamma^i_{kj}u^k.= \left(\begin{array}{ccc}-2+2\rho& 0 \\0& \ 0\end{array}\right)+ (1-2\rho +\rho ^{2}).\left(\begin{array}{ccc}0& 0 \\0& \ 1\rho\end{array}\right)
</math>

Podemos calcular el tensor gradiente como 1-contravariante 1-covariante:
:<math>\nabla\vec u=u^i_{.j}=\left(\begin{array}{ccc}-2+2\rho& 0 \\0& \ 1/\rho -2 +\rho\end{array}\right)
</math>

Las coordenadas 1-covariantes 1-contravariante coinciden por lo que la transposicion seá directa y el tensor quedara como <math>(\epsilon^i_{.j})=(\nabla \vec u + \nabla \vec u^t)/2</math> , que corrresponde con la matriz mostrada anteriormente

Sustituyendo ahora los valores en la fórmula del tensor de tensiones,y tomando \(\lambda=\mu=1\), siendo \(\lambda\) y \(\mu\) los conocidos como coeficientes de Lamé, obtenemos:

:<math>\left(\begin{array}{ccc}1/\rho-8+7\rho& 0 \\0& \ 3/\rho -8 +5\rho\end{array}\right)</math>

Las tensiones normales de ambos vectores de la base nos quedaran como:

Las tensiones normales en la dirección \(\vec g_{\rho}\) serían:\[\vec g_{\rho} \cdot \sigma \cdot \vec g_\rho =1/\rho-8+7\rho\]

Las tensiones normales en la dirección \(\vec g_\theta/\rho\) serían:\[\vec g_\theta/\rho \cdot \sigma \cdot \vec g_\theta/\rho = 3/\rho -8 +5\rho\]

A continuación se muestra el código Matlab empleado para dibujar estas tensiones normales en ambas direcciones:
[[Archivo:tensionesnormales.jpg|500px|derecha|Tensiones normales en la placa]]
{{matlab|codigo=
h=0.1; % Intervalo de separación

u=[1:0.1:2]; % Intervalo de rho [1,2]

v=[0:h:2*pi+h]; % Intervalo de theta [0,2*pi]

[uu,vv]=meshgrid(u,v); % Matrices de rho and theta

figure(1)

subplot(1,2,1)

xx=uu.*cos(vv); % Parametrización

yy=uu.*sin(vv);

f=(1./uu)-8+7.*uu; % Campo escalar

surf(xx,yy,f) % Dibujar el mallado

axis([-2,2,-2,2]) % Región seleccionada

axis equal

colorbar

title('Tension normal en dirección g sub rho')

view(2)

subplot(1,2,2)

f=(3./uu)-8+5.*uu; % Campo escalar

surf(xx,yy,f) % Dibujar el mallado

axis([-2,2,-2,2]) % Región del dibujo

axis equal

colorbar

title('Tension normal en dirección g sub theta')

view(2)

}}

Vemos como la tensiones serán mayores a medida que nos alejamos del centro de masas de nuestra placa y simétricas respecto de este ya que dependerá de la distancia de cada uno de los puntos al origen. Además se aprecia que en los bordes interiores de las placas dicha tensión será nula.

===Tensiones tangenciales===
Al realizar el calculo de las tensiones tangenciales de nuestra placa para las direcciones <math>\vec g_\rho</math> y a <math>\vec g_\theta/\rho</math>::
<math>|\sigma \cdot \vec g_\rho-(\vec g_\rho \cdot \sigma \cdot \vec g_\rho) \vec g_\rho| = 0
</math><math>|\sigma \cdot \vec g_\theta/\rho-(\vec g_\theta/\rho \cdot \sigma \cdot \vec g_\theta/\rho) \vec g_\theta/\rho|= 0 </math>

Vemos que ambas tensiones son nulas, algebraicamente hablando esto es debido a que nuestro campo de desplazamientos depende únicamente del valor de la distancia al origen. Esto implica que al estar definido dicho tensor por la divergencia y el tensor unitario, será nula la tensión tangencial siempre que el campo no dependa del angulo, ya que sera igual en todas las direcciones.

==Tensión de Von Mises==
Se define por la fórmula: <math> \sigma _{VM}=\sqrt{\frac{(\sigma 1-\sigma 2)^{2}+(\sigma 2-\sigma 3)^{2}+(\sigma 3-\sigma 1)^{2}}{2}}</math>
Esta tensión viene dada proporcionalmente a la Energía de Distorsión. En ingeniería se utiliza como un indicador de tensiones a las cuales el material falla (teoría de fallo), este cálculo se realiza sobre todo a los materiales dúctiles.
En nuestro caso el campo de tensiones de Von Misses viene dado por la fórmula <math> \sigma _{vm}=\frac{1}{\sqrt{2}}\cdot \left ( \frac{-2}{\rho }+2\cdot \rho \right ) </math>.

La representaremos mediante el siguiente programa en MATLAB:
[[Archivo:tensionesvon.jpg|500px|derecha|Tensiones Von Misses sobre la placa]]
{{matlab|codigo=
h=(0.1);
ro=1:h:2;
teta=0:h:2*pi+h;
[roro,tetateta]=meshgrid(ro,teta);
figure(1)
xx=roro.*cos(tetateta);
yy=roro.*sin(tetateta)
f=((-2./roro)-4+2.*roro)/sqrt(2); % El campo de tensiones
surf(xx,yy,f) % Dibujamos la malla
axis([-2,2,-2,2])
colorbar
title('Tensiones de Von Misses')
view(2) }}

==Masa total de la placa==
Si la densidad de la placa es <math> d(x,y,z)=xye^{-1/x^{2}} </math> y el área del anillo es <math> \pi.R^{2}-\pi.r^{2} </math> ,siendo R y r los radios exteriores e interiores respectivamente, la masa total de la placa será el producto del área y la densidad. Al tener la funcíon de la densidad, necesitamos integrarla para calcularla en todo el dominio.
Para una mejor resolución, efectuaremos el cambio a coordenadas polares:
:<math> x=\rho \cdot \cos \theta ; y=\rho \cdot \sin \theta </math>
Tendremos entonces que resolver la integral:
<math>\int_{0 }^{2\pi} \int_{1 }^{2}\rho ^{3}\cos \theta . \sin \theta .\log_{10}(\rho .\cos \theta +2)d\rho d\theta </math>
En matlab sería:
{{matlab|codigo=
h=1/500; % Número de particiones del trapecio
u=1:h:2;
v=0:h:2*pi+h;
[uu,vv]=meshgrid(ro,teta); % Malla de la placa en polares
xx=uu.*cos(vv); % Referencia en cartesianas
yy=uu.*sin(vv);
d=xx.*yy.*log(xx+2); % Función de densidad de la placa
a=h^2*d; % Masa en cada punto
mas=sum(sum(a)) % Suma de todas las masas puntuales)
}}
La masa de la placa es <math> m=2,2097.10^{-6}</math>

[[Categoría:Teoría de Campos]]
[[Categoría:TC14/15]]

Visualización de campos escalares y vectoriales en elasticidad. Grupo C4

2014-12-05T13:34:45Z

FernandoSancha: /* Masa total de la placa */

{{ TrabajoED |Visualización de campos escalares y vectoriales en elasticidad. Grupo C4 | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC14/15|2013-14]] | Alberto Rojas Rivero, Eduardo Bonet García, Fernando Sancha Domínguez, Carlos Fernández Bermejo, Ricardo Mazón Cabrera, Emilio Valero Muñoz-Rojas }}

Consideramos una placa circular comprendida entre las circunferencias de radio 1 y radio 2, centradas en el origen.
En ella vamos a suponer que tenemos definidas dos cantidades físicas: la temperatura <math>T(\rho,\theta,t)</math> ,que depende de las dos variables espaciales <math>(x,y)</math> , y el tiempo <math>t</math>, y los desplazamientos <math>\vec u(\rho,\theta)</math>.De esta forma, si definimos <math>r_0(\rho,\theta)</math> el vector de posición de los puntos de la placa antes de la deformación , la posición de cada punto <math>(\rho,\theta)</math> de la placa después de la deformación viene dada por:
<math>
\vec r (\rho,\theta)= \vec r_{0}(\rho,\theta)+\vec u(\rho,\theta).
</math>
Vamos a suponer que la fuerza aplicada sobre la placa ha provocado una desplazamiento de manera que los desplazamientos en un instante <math>t_0</math> vienen dados por:
<math>\vec u(\rho,\theta)=(1-\rho)^2 \vec g_{\theta}.
</math>

== Mallado del sólido ==
Nuestro mallado consiste en la representación de los puntos interiores de una placa solida, en forma de corona circular. Tomando como paso de muestreo h=1/10, los comandos en Matlab de la visualizacion de nuestra placa es la siguiente:
[[Archivo:Sin título.jpg|350px|miniaturadeimagen|derecha|Mallado de la corona circular.]]
{{matlab|codigo=
h=0.1 % Paso de muestreo
u=1:h:2; % Intervalo [1,2]
v=0:h:2*pi+h; % Intervalo [0,2*pi]
[uu,vv]=meshgrid(u,v); % Matrices de las coordenadas u y v
figure(1)
xx=uu.*cos(vv); % parametrizacion
yy=uu.*sin(vv);
mesh(xx,yy,0*xx) % dibuja la malla
axis([-3,3,-3,3]) % selecciona la región
view(2) }}

== Temperatura del sólido ==
La temperatura es un campo escalar representado por <math>T(x,y)=e^{-y}</math>,vamos a tratar de ver como varía la temperatura en la placa que estamos estudiando . Nuestra función depende únicamente de la variable Y, esto quiere decir que la temperatura se mantendrá constante para todo valor de X. Observamos que cuanto mayor es el valor de la Y ,menor es la temperatura , luego los puntos de máxima temperatura son los más alejados del foco( el origen) y con menor valor de la variable Y, como podemos observar en la figura de la derecha.
[[Archivo:2.jpg|325px|derecha|Distribución de la temperatura en la placa]]
{{matlab|codigo=
h=0.1 % paso de muestreo
u=1:h:2; % Intervalo[1,2]
v=0:h:2*pi+h; % Intervalo [0,2*pi]
[uu,vv]=meshgrid(u,v); % matrices
figure(1)
xx=uu.*cos(vv); % parametrización
yy=uu.*sin(vv);
f=exp(-yy); % Campo escalar
surf(xx,yy,f) % visualización
axis([-3,3,-3,3]) % región de visualización
view(2)}}

==Cálculo y representación del <math>∇T </math>==

Para estudiar como varía la temperatura al movernos por el plano calculamos el gradiente de T.
<math>∇T= -e^{-y} j </math>
En la figura de la derecha podemos observar la representación del campo vectorial <math>∇T </math> y que este es perpendicular a las curvas de nivel del campo escalar T.
[[Archivo:Anillo doblado.jpg|400px|derecha| Curvas de nivel de T y gradiente de T]]
{{matlab|codigo=
g=-exp(-yy);
subplot(1,2,1);
hold on
contour(xx,yy,f,20)
surf(xx,yy,f);
quiver3(xx,yy,f,0*xx,g,0*f);
hold off
subplot(1,2,2)
hold on
quiver(xx,yy,0*xx,Gradiente); colorbar;
contour(xx,yy,f,20)
hold off }}

==Campo de desplazamientos==
Supongamos la aplicación de una fuerza que provoca una vibración sobre la placa, de manera que, en un t aleatorio los desplazamientos vienen dados por:
<math>\vec u(\rho,\theta)=(1-\rho)^2 \vec g_{\theta}.
</math>
Para distinguir el campo de vectores en los puntos del mallado del sólido utilizaremos el código matlab siguiente:
[[Archivo:campodesplazam.jpg|500px|derecha|Distribución de la temperatura en la placCampo de desplazamientos de la placa]]
{{matlab|codigo=
u=1:h:2; % definir intervalo rho [1,2]

v=0:h:2*pi+h; % definir intervalo theta [0,2*pi]

[uu,vv]=meshgrid(u,v); % definir mallado

xx=uu.*cos(vv); % parametrizacion

yy=uu.*sin(vv);

ux=(xx.*(1-uu).^2);

uy=(yy.*(1-uu).^2);

subplot(1,2,1);

mesh(xx,yy,0.*xx)

view(2)

axis([-3,3,-3,3])

hold on

quiver(xx,yy,ux,uy)

title('campo con vectores')

subplot(1,2,2)

plot3(ux,uy,0.*xx)

title('campo de deformaciones')

}}

==Sólido antes y despues de desplazamiento==
El desplazamiento realizado por la placa viene determinado por una serie de movimientos transversales en direcciones distintas. En las imagenes adjuntas calculadas mediante Matlab, se aprecia que el desplazamiento es practicamente nulo.
[[Archivo:antesydespues.jpg|500px|derecha|Placa antes y después del desplazamiento]]
{{matlab|codigo=
h=0.1;
ro=1:h:2; % Valor de ro según el muestreo del intervalo
teta=0:h:2*pi+h; % Valor de teta según el muestreo del intervalo
[roro,tetateta]=meshgrid(ro,teta); %Matrices de la placa en coordenadas polares
figure(1)
xx=roro.*cos(tetateta); % parametrización
yy=roro.*sin(tetateta);
subplot(1,3,1) % figura solido momento inicial
mesh(xx,yy,0*xx) % mallado
axis([-4,4,-4,4]); % ejes
view(2); %visualización
desx=xx+((1-roro).^2).*(-yy); % movimiento parametrizado
desy=yy+((1-roro).^2).*(xx);
subplot(1,3,2) % figura solido en movimiento
mesh(desx,desy,yy.*0) % mallado
axis([-4,4,-4,4]);
view(2);
subplot(1,3,3); % figura solido en movimiento y campo vectorial
hold on
mesh(desx,desy,yy.*0)
axis([-4,4,-4,4]); % ejes
view(2); %
quiver(xx,yy,fx,fy) % campo vectorial del gradiente de temperaturas
axis([-4,4,-4,4]) % ejes
view(2)
hold off

}}

==Divergencia de u==
La Divergencia mide la diferencia de flujo de un campo vectorial sobre una superficie. Como vemos en la gráfica la divergencia está calculada en todos los puntos de la placa y es siempre mayor que cero, lo que nos indica que no tiene pérdidas de flujo, sino que la diferencia es cero por tanto la placa solamente tendrá fuentes.

[[Archivo:divergence.jpg|350px|derecha|Divergencia en la placa]]
{{matlab|codigo=
h=(0.1);
ro=1:h:2;
teta=0:h:2*pi+h;
[roro,tetateta]=meshgrid(ro,teta);
figure(1)
xx=roro.*cos(tetateta);
yy=roro.*sin(tetateta)
f=(1./roro)-4+3.*roro; % El campo divergencia
surf(xx,yy,f) % Dibujmos la malla
axis([-2,2,-2,2])
colorbar
title('Divergencia del campo del desplazamiento')
view(2)
}}

==Rotacional de u==
A traves del programa de Matlab, que podemos observar a continuacion, hemos calculado el rotacional del campo desplazamiento u en los puntos del solido y los hemos representado, observandose en la grafica cuales son los puntos que sufren mayor rotacional.
[[Archivo:rotacionalc4.jpg|400px|derecha|Rotacional en la placa]]
{{matlab|codigo=
h=(0.1);
ro=1:h:2;
teta=0:h:2*pi+h;
[roro,tetateta]=meshgrid(ro,teta);
figure(1)
xx=roro.*cos(tetateta);
yy=roro.*sin(tetateta);
f=(2+4.*roro.^2-6.*roro);
surf(xx,yy,f)
axis([-2,2,-2,2])
colorbar
title('Rotacional del campo del desplazamiento')
view(2)

}}

== Estudio del tensor de tensiones ==
Definamos <math>\epsilon(\vec u)=(\nabla \vec u + \nabla \vec u^t)/2</math>, la parte simétrica del tensor gradiente de <math>\vec u</math>, que se denomina tensor de deformaciones. En un medio elástico lineal, isótropo y homogéneo los desplazamientos permiten escribir el tensor de tensiones <math>\sigma_{ij}</math> a través de la fórmula::
<math>
\sigma_{ij}=\lambda \nabla \cdot \vec u \delta_{ij} + 2\mu \epsilon_{ij},
</math>
Siendo <math>\nabla \cdot \vec u = \frac{1}{\sqrt{g}}\cdot \ \frac{\partial(\sqrt{g}u^{i})}{\partial x^{i}}=\frac{1}{\rho}\cdot (\frac{\partial \rho - 2\rho ^{2}+ \rho ^{3}}{\partial \rho }+\frac{\partial 0}{\partial \theta })=\frac{1}{\rho }-4+3\rho
</math>

y calculando el gradiente de <math>\vec u</math> como:: <math>\nabla\vec u=u^i_{.j}\vec g_i \otimes \vec g^j</math> donde

<math>u^i_{.j}= \vec g^i \nabla\vec u \vec g_j=\frac{\partial u^i}{\partial x^j}+\Gamma^i_{kj}u^k.= \left(\begin{array}{ccc}-2+2\rho& 0 \\0& \ 0\end{array}\right)+ (1-2\rho +\rho ^{2}).\left(\begin{array}{ccc}0& 0 \\0& \ 1\rho\end{array}\right)
</math>

Podemos calcular el tensor gradiente como 1-contravariante 1-covariante:
:<math>\nabla\vec u=u^i_{.j}=\left(\begin{array}{ccc}-2+2\rho& 0 \\0& \ 1/\rho -2 +\rho\end{array}\right)
</math>

Las coordenadas 1-covariantes 1-contravariante coinciden por lo que la transposicion seá directa y el tensor quedara como <math>(\epsilon^i_{.j})=(\nabla \vec u + \nabla \vec u^t)/2</math> , que corrresponde con la matriz mostrada anteriormente

Sustituyendo ahora los valores en la fórmula del tensor de tensiones,y tomando \(\lambda=\mu=1\), siendo \(\lambda\) y \(\mu\) los conocidos como coeficientes de Lamé, obtenemos:

:<math>\left(\begin{array}{ccc}1/\rho-8+7\rho& 0 \\0& \ 3/\rho -8 +5\rho\end{array}\right)</math>

Las tensiones normales de ambos vectores de la base nos quedaran como:

Las tensiones normales en la dirección \(\vec g_{\rho}\) serían:\[\vec g_{\rho} \cdot \sigma \cdot \vec g_\rho =1/\rho-8+7\rho\]

Las tensiones normales en la dirección \(\vec g_\theta/\rho\) serían:\[\vec g_\theta/\rho \cdot \sigma \cdot \vec g_\theta/\rho = 3/\rho -8 +5\rho\]

A continuación se muestra el código Matlab empleado para dibujar estas tensiones normales en ambas direcciones:
[[Archivo:tensionesnormales.jpg|500px|derecha|Tensiones normales en la placa]]
{{matlab|codigo=
h=0.1; % Intervalo de separación

u=[1:0.1:2]; % Intervalo de rho [1,2]

v=[0:h:2*pi+h]; % Intervalo de theta [0,2*pi]

[uu,vv]=meshgrid(u,v); % Matrices de rho and theta

figure(1)

subplot(1,2,1)

xx=uu.*cos(vv); % Parametrización

yy=uu.*sin(vv);

f=(1./uu)-8+7.*uu; % Campo escalar

surf(xx,yy,f) % Dibujar el mallado

axis([-2,2,-2,2]) % Región seleccionada

axis equal

colorbar

title('Tension normal en dirección g sub rho')

view(2)

subplot(1,2,2)

f=(3./uu)-8+5.*uu; % Campo escalar

surf(xx,yy,f) % Dibujar el mallado

axis([-2,2,-2,2]) % Región del dibujo

axis equal

colorbar

title('Tension normal en dirección g sub theta')

view(2)

}}

Vemos como la tensiones serán mayores a medida que nos alejamos del centro de masas de nuestra placa y simétricas respecto de este ya que dependerá de la distancia de cada uno de los puntos al origen. Además se aprecia que en los bordes interiores de las placas dicha tensión será nula.

===Tensiones tangenciales===
Al realizar el calculo de las tensiones tangenciales de nuestra placa para las direcciones <math>\vec g_\rho</math> y a <math>\vec g_\theta/\rho</math>::
<math>|\sigma \cdot \vec g_\rho-(\vec g_\rho \cdot \sigma \cdot \vec g_\rho) \vec g_\rho| = 0
</math><math>|\sigma \cdot \vec g_\theta/\rho-(\vec g_\theta/\rho \cdot \sigma \cdot \vec g_\theta/\rho) \vec g_\theta/\rho|= 0 </math>

Vemos que ambas tensiones son nulas, algebraicamente hablando esto es debido a que nuestro campo de desplazamientos depende únicamente del valor de la distancia al origen. Esto implica que al estar definido dicho tensor por la divergencia y el tensor unitario, será nula la tensión tangencial siempre que el campo no dependa del angulo, ya que sera igual en todas las direcciones.

==Tensión de Von Mises==
Se define por la fórmula: <math> \sigma _{VM}=\sqrt{\frac{(\sigma 1-\sigma 2)^{2}+(\sigma 2-\sigma 3)^{2}+(\sigma 3-\sigma 1)^{2}}{2}}</math>
Esta tensión viene dada proporcionalmente a la Energía de Distorsión. En ingeniería se utiliza como un indicador de tensiones a las cuales el material falla (teoría de fallo), este cálculo se realiza sobre todo a los materiales dúctiles.
En nuestro caso el campo de tensiones de Von Misses viene dado por la fórmula <math> \sigma _{vm}=\frac{1}{\sqrt{2}}\cdot \left ( \frac{-2}{\rho }+2\cdot \rho \right ) </math>.

La representaremos mediante el siguiente programa en MATLAB:
[[Archivo:tensionesvon.jpg|500px|derecha|Tensiones Von Misses sobre la placa]]
{{matlab|codigo=
h=(0.1);
ro=1:h:2;
teta=0:h:2*pi+h;
[roro,tetateta]=meshgrid(ro,teta);
figure(1)
xx=roro.*cos(tetateta);
yy=roro.*sin(tetateta)
f=((-2./roro)-4+2.*roro)/sqrt(2); % El campo de tensiones
surf(xx,yy,f) % Dibujamos la malla
axis([-2,2,-2,2])
colorbar
title('Tensiones de Von Misses')
view(2) }}

==Masa total de la placa==
Si la densidad de la placa es <math> d(x,y,z)=xye^{-1/x^{2}} </math> y el área del anillo es <math> \pi.R^{2}-\pi.r^{2} </math> ,siendo R y r los radios exteriores e interiores respectivamente, la masa total de la placa será el producto del área y la densidad. Al tener la funcíon de la densidad, necesitamos integrarla para calcularla en todo el dominio.
Para una mejor resolución, efectuaremos el cambio a coordenadas polares:
:<math> x=\rho \cdot \cos \theta ; y=\rho \cdot \sin \theta </math>
Tendremos entonces que resolver la integral:
<math>\int_{0 }^{2\pi} \int_{1 }^{2}\rho ^{3}\cos \theta . \sin \theta .\log_{10}(\rho .\cos \theta +2)d\rho d\theta </math>
En matlab sería:
{{matlab|codigo=
h=1/500; % Número de particiones del trapecio
u=1:h:2;
v=0:h:2*pi+h;
[uu,vv]=meshgrid(ro,teta); % Malla de la placa en polares
xx=uu.*cos(vv); % Referencia en cartesianas
yy=uu.*sin(vv);
d=xx.*yy.*log(xx+2); % Función de densidad de la placa
a=h^2*d; % Masa en cada punto
mas=sum(sum(a)) % Suma de todas las masas puntuales)
}}
La masa de la placa es <math>m=2,2097.10^{-6}</math>

[[Categoría:Teoría de Campos]]
[[Categoría:TC14/15]]

Visualización de campos escalares y vectoriales en elasticidad. Grupo C4

2014-12-05T13:07:52Z

FernandoSancha: /* Divergencia de u */

{{ TrabajoED |Visualización de campos escalares y vectoriales en elasticidad. Grupo C4 | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC14/15|2013-14]] | Alberto Rojas Rivero, Eduardo Bonet García, Fernando Sancha Domínguez, Carlos Fernández Bermejo, Ricardo Mazón Cabrera, Emilio Valero Muñoz-Rojas }}

Consideramos una placa circular comprendida entre las circunferencias de radio 1 y radio 2, centradas en el origen.
En ella vamos a suponer que tenemos definidas dos cantidades físicas: la temperatura <math>T(\rho,\theta,t)</math> ,que depende de las dos variables espaciales <math>(x,y)</math> , y el tiempo <math>t</math>, y los desplazamientos <math>\vec u(\rho,\theta)</math>.De esta forma, si definimos <math>r_0(\rho,\theta)</math> el vector de posición de los puntos de la placa antes de la deformación , la posición de cada punto <math>(\rho,\theta)</math> de la placa después de la deformación viene dada por:
<math>
\vec r (\rho,\theta)= \vec r_{0}(\rho,\theta)+\vec u(\rho,\theta).
</math>
Vamos a suponer que la fuerza aplicada sobre la placa ha provocado una desplazamiento de manera que los desplazamientos en un instante <math>t_0</math> vienen dados por:
<math>\vec u(\rho,\theta)=(1-\rho)^2 \vec g_{\theta}.
</math>

== Mallado del sólido ==
Nuestro mallado consiste en la representación de los puntos interiores de una placa solida, en forma de corona circular. Tomando como paso de muestreo h=1/10, los comandos en Matlab de la visualizacion de nuestra placa es la siguiente:
[[Archivo:Sin título.jpg|350px|miniaturadeimagen|derecha|Mallado de la corona circular.]]
{{matlab|codigo=
h=0.1 % Paso de muestreo
u=1:h:2; % Intervalo [1,2]
v=0:h:2*pi+h; % Intervalo [0,2*pi]
[uu,vv]=meshgrid(u,v); % Matrices de las coordenadas u y v
figure(1)
xx=uu.*cos(vv); % parametrizacion
yy=uu.*sin(vv);
mesh(xx,yy,0*xx) % dibuja la malla
axis([-3,3,-3,3]) % selecciona la región
view(2) }}

== Temperatura del sólido ==
La temperatura es un campo escalar representado por <math>T(x,y)=e^{-y}</math>,vamos a tratar de ver como varía la temperatura en la placa que estamos estudiando . Nuestra función depende únicamente de la variable Y, esto quiere decir que la temperatura se mantendrá constante para todo valor de X. Observamos que cuanto mayor es el valor de la Y ,menor es la temperatura , luego los puntos de máxima temperatura son los más alejados del foco( el origen) y con menor valor de la variable Y, como podemos observar en la figura de la derecha.
[[Archivo:2.jpg|325px|derecha|Distribución de la temperatura en la placa]]
{{matlab|codigo=
h=0.1 % paso de muestreo
u=1:h:2; % Intervalo[1,2]
v=0:h:2*pi+h; % Intervalo [0,2*pi]
[uu,vv]=meshgrid(u,v); % matrices
figure(1)
xx=uu.*cos(vv); % parametrización
yy=uu.*sin(vv);
f=exp(-yy); % Campo escalar
surf(xx,yy,f) % visualización
axis([-3,3,-3,3]) % región de visualización
view(2)}}

==Cálculo y representación del <math>∇T </math>==

Para estudiar como varía la temperatura al movernos por el plano calculamos el gradiente de T.
<math>∇T= -e^{-y} j </math>
En la figura de la derecha podemos observar la representación del campo vectorial <math>∇T </math> y que este es perpendicular a las curvas de nivel del campo escalar T.
[[Archivo:Anillo doblado.jpg|400px|derecha| Curvas de nivel de T y gradiente de T]]
{{matlab|codigo=
g=-exp(-yy);
subplot(1,2,1);
hold on
contour(xx,yy,f,20)
surf(xx,yy,f);
quiver3(xx,yy,f,0*xx,g,0*f);
hold off
subplot(1,2,2)
hold on
quiver(xx,yy,0*xx,Gradiente); colorbar;
contour(xx,yy,f,20)
hold off }}

==Campo de desplazamientos==
Supongamos la aplicación de una fuerza que provoca una vibración sobre la placa, de manera que, en un t aleatorio los desplazamientos vienen dados por:
<math>\vec u(\rho,\theta)=(1-\rho)^2 \vec g_{\theta}.
</math>
Para distinguir el campo de vectores en los puntos del mallado del sólido utilizaremos el código matlab siguiente:
[[Archivo:campodesplazam.jpg|500px|derecha|Distribución de la temperatura en la placCampo de desplazamientos de la placa]]
{{matlab|codigo=
u=1:h:2; % definir intervalo rho [1,2]

v=0:h:2*pi+h; % definir intervalo theta [0,2*pi]

[uu,vv]=meshgrid(u,v); % definir mallado

xx=uu.*cos(vv); % parametrizacion

yy=uu.*sin(vv);

ux=(xx.*(1-uu).^2);

uy=(yy.*(1-uu).^2);

subplot(1,2,1);

mesh(xx,yy,0.*xx)

view(2)

axis([-3,3,-3,3])

hold on

quiver(xx,yy,ux,uy)

title('campo con vectores')

subplot(1,2,2)

plot3(ux,uy,0.*xx)

title('campo de deformaciones')

}}

==Sólido antes y despues de desplazamiento==
El desplazamiento realizado por la placa viene determinado por una serie de movimientos transversales en direcciones distintas. En las imagenes adjuntas calculadas mediante Matlab, se aprecia que el desplazamiento es practicamente nulo.
[[Archivo:antesydespues.jpg|500px|derecha|Placa antes y después del desplazamiento]]
{{matlab|codigo=
h=0.1;
ro=1:h:2; % Valor de ro según el muestreo del intervalo
teta=0:h:2*pi+h; % Valor de teta según el muestreo del intervalo
[roro,tetateta]=meshgrid(ro,teta); %Matrices de la placa en coordenadas polares
figure(1)
xx=roro.*cos(tetateta); % parametrización
yy=roro.*sin(tetateta);
subplot(1,3,1) % figura solido momento inicial
mesh(xx,yy,0*xx) % mallado
axis([-4,4,-4,4]); % ejes
view(2); %visualización
desx=xx+((1-roro).^2).*(-yy); % movimiento parametrizado
desy=yy+((1-roro).^2).*(xx);
subplot(1,3,2) % figura solido en movimiento
mesh(desx,desy,yy.*0) % mallado
axis([-4,4,-4,4]);
view(2);
subplot(1,3,3); % figura solido en movimiento y campo vectorial
hold on
mesh(desx,desy,yy.*0)
axis([-4,4,-4,4]); % ejes
view(2); %
quiver(xx,yy,fx,fy) % campo vectorial del gradiente de temperaturas
axis([-4,4,-4,4]) % ejes
view(2)
hold off

}}

==Divergencia de u==
La Divergencia mide la diferencia de flujo de un campo vectorial sobre una superficie. Como vemos en la gráfica la divergencia está calculada en todos los puntos de la placa y es siempre mayor que cero, lo que nos indica que no tiene pérdidas de flujo, sino que la diferencia es cero por tanto la placa solamente tendrá fuentes.

[[Archivo:divergence.jpg|350px|derecha|Divergencia en la placa]]
{{matlab|codigo=
h=(0.1);
ro=1:h:2;
teta=0:h:2*pi+h;
[roro,tetateta]=meshgrid(ro,teta);
figure(1)
xx=roro.*cos(tetateta);
yy=roro.*sin(tetateta)
f=(1./roro)-4+3.*roro; % El campo divergencia
surf(xx,yy,f) % Dibujmos la malla
axis([-2,2,-2,2])
colorbar
title('Divergencia del campo del desplazamiento')
view(2)
}}

==Rotacional de u==
A traves del programa de Matlab, que podemos observar a continuacion, hemos calculado el rotacional del campo desplazamiento u en los puntos del solido y los hemos representado, observandose en la grafica cuales son los puntos que sufren mayor rotacional.
[[Archivo:rotacionalc4.jpg|400px|derecha|Rotacional en la placa]]
{{matlab|codigo=
h=(0.1);
ro=1:h:2;
teta=0:h:2*pi+h;
[roro,tetateta]=meshgrid(ro,teta);
figure(1)
xx=roro.*cos(tetateta);
yy=roro.*sin(tetateta);
f=(2+4.*roro.^2-6.*roro);
surf(xx,yy,f)
axis([-2,2,-2,2])
colorbar
title('Rotacional del campo del desplazamiento')
view(2)

}}

== Estudio del tensor de tensiones ==
Definamos <math>\epsilon(\vec u)=(\nabla \vec u + \nabla \vec u^t)/2</math>, la parte simétrica del tensor gradiente de <math>\vec u</math>, que se denomina tensor de deformaciones. En un medio elástico lineal, isótropo y homogéneo los desplazamientos permiten escribir el tensor de tensiones <math>\sigma_{ij}</math> a través de la fórmula::
<math>
\sigma_{ij}=\lambda \nabla \cdot \vec u \delta_{ij} + 2\mu \epsilon_{ij},
</math>
Siendo <math>\nabla \cdot \vec u = \frac{1}{\sqrt{g}}\cdot \ \frac{\partial(\sqrt{g}u^{i})}{\partial x^{i}}=\frac{1}{\rho}\cdot (\frac{\partial \rho - 2\rho ^{2}+ \rho ^{3}}{\partial \rho }+\frac{\partial 0}{\partial \theta })=\frac{1}{\rho }-4+3\rho
</math>

y calculando el gradiente de <math>\vec u</math> como:: <math>\nabla\vec u=u^i_{.j}\vec g_i \otimes \vec g^j</math> donde

<math>u^i_{.j}= \vec g^i \nabla\vec u \vec g_j=\frac{\partial u^i}{\partial x^j}+\Gamma^i_{kj}u^k.= \left(\begin{array}{ccc}-2+2\rho& 0 \\0& \ 0\end{array}\right)+ (1-2\rho +\rho ^{2}).\left(\begin{array}{ccc}0& 0 \\0& \ 1\rho\end{array}\right)
</math>

Podemos calcular el tensor gradiente como 1-contravariante 1-covariante:
:<math>\nabla\vec u=u^i_{.j}=\left(\begin{array}{ccc}-2+2\rho& 0 \\0& \ 1/\rho -2 +\rho\end{array}\right)
</math>

Las coordenadas 1-covariantes 1-contravariante coinciden por lo que la transposicion seá directa y el tensor quedara como <math>(\epsilon^i_{.j})=(\nabla \vec u + \nabla \vec u^t)/2</math> , que corrresponde con la matriz mostrada anteriormente

Sustituyendo ahora los valores en la fórmula del tensor de tensiones,y tomando \(\lambda=\mu=1\), siendo \(\lambda\) y \(\mu\) los conocidos como coeficientes de Lamé, obtenemos:

:<math>\left(\begin{array}{ccc}1/\rho-8+7\rho& 0 \\0& \ 3/\rho -8 +5\rho\end{array}\right)</math>

Las tensiones normales de ambos vectores de la base nos quedaran como:

Las tensiones normales en la dirección \(\vec g_{\rho}\) serían:\[\vec g_{\rho} \cdot \sigma \cdot \vec g_\rho =1/\rho-8+7\rho\]

Las tensiones normales en la dirección \(\vec g_\theta/\rho\) serían:\[\vec g_\theta/\rho \cdot \sigma \cdot \vec g_\theta/\rho = 3/\rho -8 +5\rho\]

A continuación se muestra el código Matlab empleado para dibujar estas tensiones normales en ambas direcciones:
[[Archivo:tensionesnormales.jpg|500px|derecha|Tensiones normales en la placa]]
{{matlab|codigo=
h=0.1; % Intervalo de separación

u=[1:0.1:2]; % Intervalo de rho [1,2]

v=[0:h:2*pi+h]; % Intervalo de theta [0,2*pi]

[uu,vv]=meshgrid(u,v); % Matrices de rho and theta

figure(1)

subplot(1,2,1)

xx=uu.*cos(vv); % Parametrización

yy=uu.*sin(vv);

f=(1./uu)-8+7.*uu; % Campo escalar

surf(xx,yy,f) % Dibujar el mallado

axis([-2,2,-2,2]) % Región seleccionada

axis equal

colorbar

title('Tension normal en dirección g sub rho')

view(2)

subplot(1,2,2)

f=(3./uu)-8+5.*uu; % Campo escalar

surf(xx,yy,f) % Dibujar el mallado

axis([-2,2,-2,2]) % Región del dibujo

axis equal

colorbar

title('Tension normal en dirección g sub theta')

view(2)

}}

Vemos como la tensiones serán mayores a medida que nos alejamos del centro de masas de nuestra placa y simétricas respecto de este ya que dependerá de la distancia de cada uno de los puntos al origen. Además se aprecia que en los bordes interiores de las placas dicha tensión será nula.

===Tensiones tangenciales===
Al realizar el calculo de las tensiones tangenciales de nuestra placa para las direcciones <math>\vec g_\rho</math> y a <math>\vec g_\theta/\rho</math>::
<math>|\sigma \cdot \vec g_\rho-(\vec g_\rho \cdot \sigma \cdot \vec g_\rho) \vec g_\rho| = 0
</math><math>|\sigma \cdot \vec g_\theta/\rho-(\vec g_\theta/\rho \cdot \sigma \cdot \vec g_\theta/\rho) \vec g_\theta/\rho|= 0 </math>

Vemos que ambas tensiones son nulas, algebraicamente hablando esto es debido a que nuestro campo de desplazamientos depende únicamente del valor de la distancia al origen. Esto implica que al estar definido dicho tensor por la divergencia y el tensor unitario, será nula la tensión tangencial siempre que el campo no dependa del angulo, ya que sera igual en todas las direcciones.

==Tensión de Von Mises==
Se define por la fórmula: <math> \sigma _{VM}=\sqrt{\frac{(\sigma 1-\sigma 2)^{2}+(\sigma 2-\sigma 3)^{2}+(\sigma 3-\sigma 1)^{2}}{2}}</math>
Esta tensión viene dada proporcionalmente a la Energía de Distorsión. En ingeniería se utiliza como un indicador de tensiones a las cuales el material falla (teoría de fallo), este cálculo se realiza sobre todo a los materiales dúctiles.
En nuestro caso el campo de tensiones de Von Misses viene dado por la fórmula <math> \sigma _{vm}=\frac{1}{\sqrt{2}}\cdot \left ( \frac{-2}{\rho }+2\cdot \rho \right ) </math>.

La representaremos mediante el siguiente programa en MATLAB:
[[Archivo:tensionesvon.jpg|500px|derecha|Tensiones Von Misses sobre la placa]]
{{matlab|codigo=
h=(0.1);
ro=1:h:2;
teta=0:h:2*pi+h;
[roro,tetateta]=meshgrid(ro,teta);
figure(1)
xx=roro.*cos(tetateta);
yy=roro.*sin(tetateta)
f=((-2./roro)-4+2.*roro)/sqrt(2); % El campo de tensiones
surf(xx,yy,f) % Dibujamos la malla
axis([-2,2,-2,2])
colorbar
title('Tensiones de Von Misses')
view(2) }}

==Masa total de la placa==
Si la densidad de la placa es <math> d(x,y,z)=xye^{-1/x^{2}} </math> y el área del anillo es <math> \pi.R^{2}-\pi.r^{2} </math> ,siendo R y r los radios exteriores e interiores respectivamente, la masa total de la placa será el producto del área y la densidad. Al tener la funcíon de la densidad, necesitamos integrarla para calcularla en todo el dominio.
Para una mejor resolución, efectuaremos el cambio a coordenadas polares:
:<math> x=\rho \cdot \cos \theta ; y=\rho \cdot \sin \theta </math>
Tendremos entonces que resolver la integral:
<math>\int_{0 }^{2\pi} \int_{1 }^{2}\rho ^{3}\cos \theta . \sin \theta .\log_{10}(\rho .\cos \theta +2)d\rho d\theta </math>
En matlab sería:
{{matlab|codigo=
h=1/100
u=1:h:2
v=0:h:2*pi
[uu,vv]=meshgrid(u,v)
f=(uu.^3).*sin(vv).*cos(vv).*log(uu.*cos(vv)+2)
a=(h^2)*f
masa=sum(sum(a))
}}
La masa de la placa es <math>m=5,0709.10^{-5}</math>

[[Categoría:Teoría de Campos]]
[[Categoría:TC14/15]]

Visualización de campos escalares y vectoriales en elasticidad. Grupo C4

2014-12-05T12:59:30Z

FernandoSancha: /* Tensión de Von Mises */

{{ TrabajoED |Visualización de campos escalares y vectoriales en elasticidad. Grupo C4 | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC14/15|2013-14]] | Alberto Rojas Rivero, Eduardo Bonet García, Fernando Sancha Domínguez, Carlos Fernández Bermejo, Ricardo Mazón Cabrera, Emilio Valero Muñoz-Rojas }}

Consideramos una placa circular comprendida entre las circunferencias de radio 1 y radio 2, centradas en el origen.
En ella vamos a suponer que tenemos definidas dos cantidades físicas: la temperatura <math>T(\rho,\theta,t)</math> ,que depende de las dos variables espaciales <math>(x,y)</math> , y el tiempo <math>t</math>, y los desplazamientos <math>\vec u(\rho,\theta)</math>.De esta forma, si definimos <math>r_0(\rho,\theta)</math> el vector de posición de los puntos de la placa antes de la deformación , la posición de cada punto <math>(\rho,\theta)</math> de la placa después de la deformación viene dada por:
<math>
\vec r (\rho,\theta)= \vec r_{0}(\rho,\theta)+\vec u(\rho,\theta).
</math>
Vamos a suponer que la fuerza aplicada sobre la placa ha provocado una desplazamiento de manera que los desplazamientos en un instante <math>t_0</math> vienen dados por:
<math>\vec u(\rho,\theta)=(1-\rho)^2 \vec g_{\theta}.
</math>

== Mallado del sólido ==
Nuestro mallado consiste en la representación de los puntos interiores de una placa solida, en forma de corona circular. Tomando como paso de muestreo h=1/10, los comandos en Matlab de la visualizacion de nuestra placa es la siguiente:
[[Archivo:Sin título.jpg|350px|miniaturadeimagen|derecha|Mallado de la corona circular.]]
{{matlab|codigo=
h=0.1 % Paso de muestreo
u=1:h:2; % Intervalo [1,2]
v=0:h:2*pi+h; % Intervalo [0,2*pi]
[uu,vv]=meshgrid(u,v); % Matrices de las coordenadas u y v
figure(1)
xx=uu.*cos(vv); % parametrizacion
yy=uu.*sin(vv);
mesh(xx,yy,0*xx) % dibuja la malla
axis([-3,3,-3,3]) % selecciona la región
view(2) }}

== Temperatura del sólido ==
La temperatura es un campo escalar representado por <math>T(x,y)=e^{-y}</math>,vamos a tratar de ver como varía la temperatura en la placa que estamos estudiando . Nuestra función depende únicamente de la variable Y, esto quiere decir que la temperatura se mantendrá constante para todo valor de X. Observamos que cuanto mayor es el valor de la Y ,menor es la temperatura , luego los puntos de máxima temperatura son los más alejados del foco( el origen) y con menor valor de la variable Y, como podemos observar en la figura de la derecha.
[[Archivo:2.jpg|325px|derecha|Distribución de la temperatura en la placa]]
{{matlab|codigo=
h=0.1 % paso de muestreo
u=1:h:2; % Intervalo[1,2]
v=0:h:2*pi+h; % Intervalo [0,2*pi]
[uu,vv]=meshgrid(u,v); % matrices
figure(1)
xx=uu.*cos(vv); % parametrización
yy=uu.*sin(vv);
f=exp(-yy); % Campo escalar
surf(xx,yy,f) % visualización
axis([-3,3,-3,3]) % región de visualización
view(2)}}

==Cálculo y representación del <math>∇T </math>==

Para estudiar como varía la temperatura al movernos por el plano calculamos el gradiente de T.
<math>∇T= -e^{-y} j </math>
En la figura de la derecha podemos observar la representación del campo vectorial <math>∇T </math> y que este es perpendicular a las curvas de nivel del campo escalar T.
[[Archivo:Anillo doblado.jpg|400px|derecha| Curvas de nivel de T y gradiente de T]]
{{matlab|codigo=
g=-exp(-yy);
subplot(1,2,1);
hold on
contour(xx,yy,f,20)
surf(xx,yy,f);
quiver3(xx,yy,f,0*xx,g,0*f);
hold off
subplot(1,2,2)
hold on
quiver(xx,yy,0*xx,Gradiente); colorbar;
contour(xx,yy,f,20)
hold off }}

==Campo de desplazamientos==
Supongamos la aplicación de una fuerza que provoca una vibración sobre la placa, de manera que, en un t aleatorio los desplazamientos vienen dados por:
<math>\vec u(\rho,\theta)=(1-\rho)^2 \vec g_{\theta}.
</math>
Para distinguir el campo de vectores en los puntos del mallado del sólido utilizaremos el código matlab siguiente:
[[Archivo:campodesplazam.jpg|500px|derecha|Distribución de la temperatura en la placCampo de desplazamientos de la placa]]
{{matlab|codigo=
u=1:h:2; % definir intervalo rho [1,2]

v=0:h:2*pi+h; % definir intervalo theta [0,2*pi]

[uu,vv]=meshgrid(u,v); % definir mallado

xx=uu.*cos(vv); % parametrizacion

yy=uu.*sin(vv);

ux=(xx.*(1-uu).^2);

uy=(yy.*(1-uu).^2);

subplot(1,2,1);

mesh(xx,yy,0.*xx)

view(2)

axis([-3,3,-3,3])

hold on

quiver(xx,yy,ux,uy)

title('campo con vectores')

subplot(1,2,2)

plot3(ux,uy,0.*xx)

title('campo de deformaciones')

}}

==Sólido antes y despues de desplazamiento==
El desplazamiento realizado por la placa viene determinado por una serie de movimientos transversales en direcciones distintas. En las imagenes adjuntas calculadas mediante Matlab, se aprecia que el desplazamiento es practicamente nulo.
[[Archivo:antesydespues.jpg|500px|derecha|Placa antes y después del desplazamiento]]
{{matlab|codigo=
h=0.1;
ro=1:h:2; % Valor de ro según el muestreo del intervalo
teta=0:h:2*pi+h; % Valor de teta según el muestreo del intervalo
[roro,tetateta]=meshgrid(ro,teta); %Matrices de la placa en coordenadas polares
figure(1)
xx=roro.*cos(tetateta); % parametrización
yy=roro.*sin(tetateta);
subplot(1,3,1) % figura solido momento inicial
mesh(xx,yy,0*xx) % mallado
axis([-4,4,-4,4]); % ejes
view(2); %visualización
desx=xx+((1-roro).^2).*(-yy); % movimiento parametrizado
desy=yy+((1-roro).^2).*(xx);
subplot(1,3,2) % figura solido en movimiento
mesh(desx,desy,yy.*0) % mallado
axis([-4,4,-4,4]);
view(2);
subplot(1,3,3); % figura solido en movimiento y campo vectorial
hold on
mesh(desx,desy,yy.*0)
axis([-4,4,-4,4]); % ejes
view(2); %
quiver(xx,yy,fx,fy) % campo vectorial del gradiente de temperaturas
axis([-4,4,-4,4]) % ejes
view(2)
hold off

}}

==Divergencia de u==
La divergencia mide el cambio de volumen local, en este caso en nuestra placa, debido al desplazamiento.
Con el grafico adjunto observamos que todos los puntos tienen el mismo valor, cero.
Como demostramos en el apartado anterior, nuestra placa esta sometida a desplazamientos transversales sin variacion de volumen, por lo que no hay cambio del mismo.
[[Archivo:divergence.jpg|350px|derecha|Divergencia en la placa]]
{{matlab|codigo=
h=(0.1);
ro=1:h:2;
teta=0:h:2*pi+h;
[roro,tetateta]=meshgrid(ro,teta);
figure(1)
xx=roro.*cos(tetateta);
yy=roro.*sin(tetateta)
f=(1./roro)-4+3.*roro; % El campo divergencia
surf(xx,yy,f) % Dibujmos la malla
axis([-2,2,-2,2])
colorbar
title('Divergencia del campo del desplazamiento')
view(2)
}}

==Rotacional de u==
A traves del programa de Matlab, que podemos observar a continuacion, hemos calculado el rotacional del campo desplazamiento u en los puntos del solido y los hemos representado, observandose en la grafica cuales son los puntos que sufren mayor rotacional.
[[Archivo:rotacionalc4.jpg|400px|derecha|Rotacional en la placa]]
{{matlab|codigo=
h=(0.1);
ro=1:h:2;
teta=0:h:2*pi+h;
[roro,tetateta]=meshgrid(ro,teta);
figure(1)
xx=roro.*cos(tetateta);
yy=roro.*sin(tetateta);
f=(2+4.*roro.^2-6.*roro);
surf(xx,yy,f)
axis([-2,2,-2,2])
colorbar
title('Rotacional del campo del desplazamiento')
view(2)

}}

== Estudio del tensor de tensiones ==
Definamos <math>\epsilon(\vec u)=(\nabla \vec u + \nabla \vec u^t)/2</math>, la parte simétrica del tensor gradiente de <math>\vec u</math>, que se denomina tensor de deformaciones. En un medio elástico lineal, isótropo y homogéneo los desplazamientos permiten escribir el tensor de tensiones <math>\sigma_{ij}</math> a través de la fórmula::
<math>
\sigma_{ij}=\lambda \nabla \cdot \vec u \delta_{ij} + 2\mu \epsilon_{ij},
</math>
Siendo <math>\nabla \cdot \vec u = \frac{1}{\sqrt{g}}\cdot \ \frac{\partial(\sqrt{g}u^{i})}{\partial x^{i}}=\frac{1}{\rho}\cdot (\frac{\partial \rho - 2\rho ^{2}+ \rho ^{3}}{\partial \rho }+\frac{\partial 0}{\partial \theta })=\frac{1}{\rho }-4+3\rho
</math>

y calculando el gradiente de <math>\vec u</math> como:: <math>\nabla\vec u=u^i_{.j}\vec g_i \otimes \vec g^j</math> donde

<math>u^i_{.j}= \vec g^i \nabla\vec u \vec g_j=\frac{\partial u^i}{\partial x^j}+\Gamma^i_{kj}u^k.= \left(\begin{array}{ccc}-2+2\rho& 0 \\0& \ 0\end{array}\right)+ (1-2\rho +\rho ^{2}).\left(\begin{array}{ccc}0& 0 \\0& \ 1\rho\end{array}\right)
</math>

Podemos calcular el tensor gradiente como 1-contravariante 1-covariante:
:<math>\nabla\vec u=u^i_{.j}=\left(\begin{array}{ccc}-2+2\rho& 0 \\0& \ 1/\rho -2 +\rho\end{array}\right)
</math>

Las coordenadas 1-covariantes 1-contravariante coinciden por lo que la transposicion seá directa y el tensor quedara como <math>(\epsilon^i_{.j})=(\nabla \vec u + \nabla \vec u^t)/2</math> , que corrresponde con la matriz mostrada anteriormente

Sustituyendo ahora los valores en la fórmula del tensor de tensiones,y tomando \(\lambda=\mu=1\), siendo \(\lambda\) y \(\mu\) los conocidos como coeficientes de Lamé, obtenemos:

:<math>\left(\begin{array}{ccc}1/\rho-8+7\rho& 0 \\0& \ 3/\rho -8 +5\rho\end{array}\right)</math>

Las tensiones normales de ambos vectores de la base nos quedaran como:

Las tensiones normales en la dirección \(\vec g_{\rho}\) serían:\[\vec g_{\rho} \cdot \sigma \cdot \vec g_\rho =1/\rho-8+7\rho\]

Las tensiones normales en la dirección \(\vec g_\theta/\rho\) serían:\[\vec g_\theta/\rho \cdot \sigma \cdot \vec g_\theta/\rho = 3/\rho -8 +5\rho\]

A continuación se muestra el código Matlab empleado para dibujar estas tensiones normales en ambas direcciones:
[[Archivo:tensionesnormales.jpg|500px|derecha|Tensiones normales en la placa]]
{{matlab|codigo=
h=0.1; % Intervalo de separación

u=[1:0.1:2]; % Intervalo de rho [1,2]

v=[0:h:2*pi+h]; % Intervalo de theta [0,2*pi]

[uu,vv]=meshgrid(u,v); % Matrices de rho and theta

figure(1)

subplot(1,2,1)

xx=uu.*cos(vv); % Parametrización

yy=uu.*sin(vv);

f=(1./uu)-8+7.*uu; % Campo escalar

surf(xx,yy,f) % Dibujar el mallado

axis([-2,2,-2,2]) % Región seleccionada

axis equal

colorbar

title('Tension normal en dirección g sub rho')

view(2)

subplot(1,2,2)

f=(3./uu)-8+5.*uu; % Campo escalar

surf(xx,yy,f) % Dibujar el mallado

axis([-2,2,-2,2]) % Región del dibujo

axis equal

colorbar

title('Tension normal en dirección g sub theta')

view(2)

}}

Vemos como la tensiones serán mayores a medida que nos alejamos del centro de masas de nuestra placa y simétricas respecto de este ya que dependerá de la distancia de cada uno de los puntos al origen. Además se aprecia que en los bordes interiores de las placas dicha tensión será nula.

===Tensiones tangenciales===
Al realizar el calculo de las tensiones tangenciales de nuestra placa para las direcciones <math>\vec g_\rho</math> y a <math>\vec g_\theta/\rho</math>::
<math>|\sigma \cdot \vec g_\rho-(\vec g_\rho \cdot \sigma \cdot \vec g_\rho) \vec g_\rho| = 0
</math><math>|\sigma \cdot \vec g_\theta/\rho-(\vec g_\theta/\rho \cdot \sigma \cdot \vec g_\theta/\rho) \vec g_\theta/\rho|= 0 </math>

Vemos que ambas tensiones son nulas, algebraicamente hablando esto es debido a que nuestro campo de desplazamientos depende únicamente del valor de la distancia al origen. Esto implica que al estar definido dicho tensor por la divergencia y el tensor unitario, será nula la tensión tangencial siempre que el campo no dependa del angulo, ya que sera igual en todas las direcciones.

==Tensión de Von Mises==
Se define por la fórmula: <math> \sigma _{VM}=\sqrt{\frac{(\sigma 1-\sigma 2)^{2}+(\sigma 2-\sigma 3)^{2}+(\sigma 3-\sigma 1)^{2}}{2}}</math>
Esta tensión viene dada proporcionalmente a la Energía de Distorsión. En ingeniería se utiliza como un indicador de tensiones a las cuales el material falla (teoría de fallo), este cálculo se realiza sobre todo a los materiales dúctiles.
En nuestro caso el campo de tensiones de Von Misses viene dado por la fórmula <math> \sigma _{vm}=\frac{1}{\sqrt{2}}\cdot \left ( \frac{-2}{\rho }+2\cdot \rho \right ) </math>.

La representaremos mediante el siguiente programa en MATLAB:
[[Archivo:tensionesvon.jpg|500px|derecha|Tensiones Von Misses sobre la placa]]
{{matlab|codigo=
h=(0.1);
ro=1:h:2;
teta=0:h:2*pi+h;
[roro,tetateta]=meshgrid(ro,teta);
figure(1)
xx=roro.*cos(tetateta);
yy=roro.*sin(tetateta)
f=((-2./roro)-4+2.*roro)/sqrt(2); % El campo de tensiones
surf(xx,yy,f) % Dibujamos la malla
axis([-2,2,-2,2])
colorbar
title('Tensiones de Von Misses')
view(2) }}

==Masa total de la placa==
Si la densidad de la placa es <math> d(x,y,z)=xye^{-1/x^{2}} </math> y el área del anillo es <math> \pi.R^{2}-\pi.r^{2} </math> ,siendo R y r los radios exteriores e interiores respectivamente, la masa total de la placa será el producto del área y la densidad. Al tener la funcíon de la densidad, necesitamos integrarla para calcularla en todo el dominio.
Para una mejor resolución, efectuaremos el cambio a coordenadas polares:
:<math> x=\rho \cdot \cos \theta ; y=\rho \cdot \sin \theta </math>
Tendremos entonces que resolver la integral:
<math>\int_{0 }^{2\pi} \int_{1 }^{2}\rho ^{3}\cos \theta . \sin \theta .\log_{10}(\rho .\cos \theta +2)d\rho d\theta </math>
En matlab sería:
{{matlab|codigo=
h=1/100
u=1:h:2
v=0:h:2*pi
[uu,vv]=meshgrid(u,v)
f=(uu.^3).*sin(vv).*cos(vv).*log(uu.*cos(vv)+2)
a=(h^2)*f
masa=sum(sum(a))
}}
La masa de la placa es <math>m=5,0709.10^{-5}</math>

[[Categoría:Teoría de Campos]]
[[Categoría:TC14/15]]

Visualización de campos escalares y vectoriales en elasticidad. Grupo C4

2014-12-05T12:27:41Z

FernandoSancha:

{{ TrabajoED |Visualización de campos escalares y vectoriales en elasticidad. Grupo C4 | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC14/15|2013-14]] | Alberto Rojas Rivero, Eduardo Bonet García, Fernando Sancha Domínguez, Carlos Fernández Bermejo, Ricardo Mazón Cabrera, Emilio Valero Muñoz-Rojas }}

Consideramos una placa circular comprendida entre las circunferencias de radio 1 y radio 2, centradas en el origen.
En ella vamos a suponer que tenemos definidas dos cantidades físicas: la temperatura <math>T(\rho,\theta,t)</math> ,que depende de las dos variables espaciales <math>(x,y)</math> , y el tiempo <math>t</math>, y los desplazamientos <math>\vec u(\rho,\theta)</math>.De esta forma, si definimos <math>r_0(\rho,\theta)</math> el vector de posición de los puntos de la placa antes de la deformación , la posición de cada punto <math>(\rho,\theta)</math> de la placa después de la deformación viene dada por:
<math>
\vec r (\rho,\theta)= \vec r_{0}(\rho,\theta)+\vec u(\rho,\theta).
</math>
Vamos a suponer que la fuerza aplicada sobre la placa ha provocado una desplazamiento de manera que los desplazamientos en un instante <math>t_0</math> vienen dados por:
<math>\vec u(\rho,\theta)=(1-\rho)^2 \vec g_{\theta}.
</math>

== Mallado del sólido ==
Nuestro mallado consiste en la representación de los puntos interiores de una placa solida, en forma de corona circular. Tomando como paso de muestreo h=1/10, los comandos en Matlab de la visualizacion de nuestra placa es la siguiente:
[[Archivo:Sin título.jpg|350px|miniaturadeimagen|derecha|Mallado de la corona circular.]]
{{matlab|codigo=
h=0.1 % Paso de muestreo
u=1:h:2; % Intervalo [1,2]
v=0:h:2*pi+h; % Intervalo [0,2*pi]
[uu,vv]=meshgrid(u,v); % Matrices de las coordenadas u y v
figure(1)
xx=uu.*cos(vv); % parametrizacion
yy=uu.*sin(vv);
mesh(xx,yy,0*xx) % dibuja la malla
axis([-3,3,-3,3]) % selecciona la región
view(2) }}

== Temperatura del sólido ==
La temperatura es un campo escalar representado por <math>T(x,y)=e^{-y}</math>,vamos a tratar de ver como varía la temperatura en la placa que estamos estudiando . Nuestra función depende únicamente de la variable Y, esto quiere decir que la temperatura se mantendrá constante para todo valor de X. Observamos que cuanto mayor es el valor de la Y ,menor es la temperatura , luego los puntos de máxima temperatura son los más alejados del foco( el origen) y con menor valor de la variable Y, como podemos observar en la figura de la derecha.
[[Archivo:2.jpg|325px|derecha|Distribución de la temperatura en la placa]]
{{matlab|codigo=
h=0.1 % paso de muestreo
u=1:h:2; % Intervalo[1,2]
v=0:h:2*pi+h; % Intervalo [0,2*pi]
[uu,vv]=meshgrid(u,v); % matrices
figure(1)
xx=uu.*cos(vv); % parametrización
yy=uu.*sin(vv);
f=exp(-yy); % Campo escalar
surf(xx,yy,f) % visualización
axis([-3,3,-3,3]) % región de visualización
view(2)}}

==Cálculo y representación del <math>∇T </math>==

Para estudiar como varía la temperatura al movernos por el plano calculamos el gradiente de T.
<math>∇T= -e^{-y} j </math>
En la figura de la derecha podemos observar la representación del campo vectorial <math>∇T </math> y que este es perpendicular a las curvas de nivel del campo escalar T.
[[Archivo:Anillo doblado.jpg|400px|derecha| Curvas de nivel de T y gradiente de T]]
{{matlab|codigo=
g=-exp(-yy);
subplot(1,2,1);
hold on
contour(xx,yy,f,20)
surf(xx,yy,f);
quiver3(xx,yy,f,0*xx,g,0*f);
hold off
subplot(1,2,2)
hold on
quiver(xx,yy,0*xx,Gradiente); colorbar;
contour(xx,yy,f,20)
hold off }}

==Campo de desplazamientos==
Supongamos la aplicación de una fuerza que provoca una vibración sobre la placa, de manera que, en un t aleatorio los desplazamientos vienen dados por:
<math>\vec u(\rho,\theta)=(1-\rho)^2 \vec g_{\theta}.
</math>
Para distinguir el campo de vectores en los puntos del mallado del sólido utilizaremos el código matlab siguiente:
[[Archivo:campodesplazam.jpg|500px|derecha|Distribución de la temperatura en la placCampo de desplazamientos de la placa]]
{{matlab|codigo=
u=1:h:2; % definir intervalo rho [1,2]

v=0:h:2*pi+h; % definir intervalo theta [0,2*pi]

[uu,vv]=meshgrid(u,v); % definir mallado

xx=uu.*cos(vv); % parametrizacion

yy=uu.*sin(vv);

ux=(xx.*(1-uu).^2);

uy=(yy.*(1-uu).^2);

subplot(1,2,1);

mesh(xx,yy,0.*xx)

view(2)

axis([-3,3,-3,3])

hold on

quiver(xx,yy,ux,uy)

title('campo con vectores')

subplot(1,2,2)

plot3(ux,uy,0.*xx)

title('campo de deformaciones')

}}

==Sólido antes y despues de desplazamiento==
El desplazamiento realizado por la placa viene determinado por una serie de movimientos transversales en direcciones distintas. En las imagenes adjuntas calculadas mediante Matlab, se aprecia que el desplazamiento es practicamente nulo.
[[Archivo:antesydespues.jpg|500px|derecha|Placa antes y después del desplazamiento]]
{{matlab|codigo=
h=0.1;
ro=1:h:2; % Valor de ro según el muestreo del intervalo
teta=0:h:2*pi+h; % Valor de teta según el muestreo del intervalo
[roro,tetateta]=meshgrid(ro,teta); %Matrices de la placa en coordenadas polares
figure(1)
xx=roro.*cos(tetateta); % parametrización
yy=roro.*sin(tetateta);
subplot(1,3,1) % figura solido momento inicial
mesh(xx,yy,0*xx) % mallado
axis([-4,4,-4,4]); % ejes
view(2); %visualización
desx=xx+((1-roro).^2).*(-yy); % movimiento parametrizado
desy=yy+((1-roro).^2).*(xx);
subplot(1,3,2) % figura solido en movimiento
mesh(desx,desy,yy.*0) % mallado
axis([-4,4,-4,4]);
view(2);
subplot(1,3,3); % figura solido en movimiento y campo vectorial
hold on
mesh(desx,desy,yy.*0)
axis([-4,4,-4,4]); % ejes
view(2); %
quiver(xx,yy,fx,fy) % campo vectorial del gradiente de temperaturas
axis([-4,4,-4,4]) % ejes
view(2)
hold off

}}

==Divergencia de u==
La divergencia mide el cambio de volumen local, en este caso en nuestra placa, debido al desplazamiento.
Con el grafico adjunto observamos que todos los puntos tienen el mismo valor, cero.
Como demostramos en el apartado anterior, nuestra placa esta sometida a desplazamientos transversales sin variacion de volumen, por lo que no hay cambio del mismo.
[[Archivo:divergence.jpg|350px|derecha|Divergencia en la placa]]
{{matlab|codigo=
h=(0.1);
ro=1:h:2;
teta=0:h:2*pi+h;
[roro,tetateta]=meshgrid(ro,teta);
figure(1)
xx=roro.*cos(tetateta);
yy=roro.*sin(tetateta)
f=(1./roro)-4+3.*roro; % El campo divergencia
surf(xx,yy,f) % Dibujmos la malla
axis([-2,2,-2,2])
colorbar
title('Divergencia del campo del desplazamiento')
view(2)
}}

==Rotacional de u==
A traves del programa de Matlab, que podemos observar a continuacion, hemos calculado el rotacional del campo desplazamiento u en los puntos del solido y los hemos representado, observandose en la grafica cuales son los puntos que sufren mayor rotacional.
[[Archivo:rotacionalc4.jpg|400px|derecha|Rotacional en la placa]]
{{matlab|codigo=
h=(0.1);
ro=1:h:2;
teta=0:h:2*pi+h;
[roro,tetateta]=meshgrid(ro,teta);
figure(1)
xx=roro.*cos(tetateta);
yy=roro.*sin(tetateta);
f=(2+4.*roro.^2-6.*roro);
surf(xx,yy,f)
axis([-2,2,-2,2])
colorbar
title('Rotacional del campo del desplazamiento')
view(2)

}}

== Estudio del tensor de tensiones ==
Definamos <math>\epsilon(\vec u)=(\nabla \vec u + \nabla \vec u^t)/2</math>, la parte simétrica del tensor gradiente de <math>\vec u</math>, que se denomina tensor de deformaciones. En un medio elástico lineal, isótropo y homogéneo los desplazamientos permiten escribir el tensor de tensiones <math>\sigma_{ij}</math> a través de la fórmula::
<math>
\sigma_{ij}=\lambda \nabla \cdot \vec u \delta_{ij} + 2\mu \epsilon_{ij},
</math>
Siendo <math>\nabla \cdot \vec u = \frac{1}{\sqrt{g}}\cdot \ \frac{\partial(\sqrt{g}u^{i})}{\partial x^{i}}=\frac{1}{\rho}\cdot (\frac{\partial \rho - 2\rho ^{2}+ \rho ^{3}}{\partial \rho }+\frac{\partial 0}{\partial \theta })=\frac{1}{\rho }-4+3\rho
</math>

y calculando el gradiente de <math>\vec u</math> como:: <math>\nabla\vec u=u^i_{.j}\vec g_i \otimes \vec g^j</math> donde

<math>u^i_{.j}= \vec g^i \nabla\vec u \vec g_j=\frac{\partial u^i}{\partial x^j}+\Gamma^i_{kj}u^k.= \left(\begin{array}{ccc}-2+2\rho& 0 \\0& \ 0\end{array}\right)+ (1-2\rho +\rho ^{2}).\left(\begin{array}{ccc}0& 0 \\0& \ 1\rho\end{array}\right)
</math>

Podemos calcular el tensor gradiente como 1-contravariante 1-covariante:
:<math>\nabla\vec u=u^i_{.j}=\left(\begin{array}{ccc}-2+2\rho& 0 \\0& \ 1/\rho -2 +\rho\end{array}\right)
</math>

Las coordenadas 1-covariantes 1-contravariante coinciden por lo que la transposicion seá directa y el tensor quedara como <math>(\epsilon^i_{.j})=(\nabla \vec u + \nabla \vec u^t)/2</math> , que corrresponde con la matriz mostrada anteriormente

Sustituyendo ahora los valores en la fórmula del tensor de tensiones,y tomando \(\lambda=\mu=1\), siendo \(\lambda\) y \(\mu\) los conocidos como coeficientes de Lamé, obtenemos:

:<math>\left(\begin{array}{ccc}1/\rho-8+7\rho& 0 \\0& \ 3/\rho -8 +5\rho\end{array}\right)</math>

Las tensiones normales de ambos vectores de la base nos quedaran como:

Las tensiones normales en la dirección \(\vec g_{\rho}\) serían:\[\vec g_{\rho} \cdot \sigma \cdot \vec g_\rho =1/\rho-8+7\rho\]

Las tensiones normales en la dirección \(\vec g_\theta/\rho\) serían:\[\vec g_\theta/\rho \cdot \sigma \cdot \vec g_\theta/\rho = 3/\rho -8 +5\rho\]

A continuación se muestra el código Matlab empleado para dibujar estas tensiones normales en ambas direcciones:
[[Archivo:tensionesnormales.jpg|500px|derecha|Tensiones normales en la placa]]
{{matlab|codigo=
h=0.1; % Intervalo de separación

u=[1:0.1:2]; % Intervalo de rho [1,2]

v=[0:h:2*pi+h]; % Intervalo de theta [0,2*pi]

[uu,vv]=meshgrid(u,v); % Matrices de rho and theta

figure(1)

subplot(1,2,1)

xx=uu.*cos(vv); % Parametrización

yy=uu.*sin(vv);

f=(1./uu)-8+7.*uu; % Campo escalar

surf(xx,yy,f) % Dibujar el mallado

axis([-2,2,-2,2]) % Región seleccionada

axis equal

colorbar

title('Tension normal en dirección g sub rho')

view(2)

subplot(1,2,2)

f=(3./uu)-8+5.*uu; % Campo escalar

surf(xx,yy,f) % Dibujar el mallado

axis([-2,2,-2,2]) % Región del dibujo

axis equal

colorbar

title('Tension normal en dirección g sub theta')

view(2)

}}

Vemos como la tensiones serán mayores a medida que nos alejamos del centro de masas de nuestra placa y simétricas respecto de este ya que dependerá de la distancia de cada uno de los puntos al origen. Además se aprecia que en los bordes interiores de las placas dicha tensión será nula.

===Tensiones tangenciales===
Al realizar el calculo de las tensiones tangenciales de nuestra placa para las direcciones <math>\vec g_\rho</math> y a <math>\vec g_\theta/\rho</math>::
<math>|\sigma \cdot \vec g_\rho-(\vec g_\rho \cdot \sigma \cdot \vec g_\rho) \vec g_\rho| = 0
</math><math>|\sigma \cdot \vec g_\theta/\rho-(\vec g_\theta/\rho \cdot \sigma \cdot \vec g_\theta/\rho) \vec g_\theta/\rho|= 0 </math>

Vemos que ambas tensiones son nulas, algebraicamente hablando esto es debido a que nuestro campo de desplazamientos depende únicamente del valor de la distancia al origen. Esto implica que al estar definido dicho tensor por la divergencia y el tensor unitario, será nula la tensión tangencial siempre que el campo no dependa del angulo, ya que sera igual en todas las direcciones.

==Tensión de Von Mises==
Se define por la fórmula: <math> \sigma _{VM}=\sqrt{\frac{(\sigma 1-\sigma 2)^{2}+(\sigma 2-\sigma 3)^{2}+(\sigma 3-\sigma 1)^{2}}{2}}</math>
En nuestro caso el campo de tensiones de Von Misses viene dado por la fórmula <math> \sigma _{vm}=\frac{1}{\sqrt{2}}\cdot \left ( \frac{-2}{\rho }+2\cdot \rho \right ) </math>.

La representaremos mediante el siguiente programa en MATLAB:
[[Archivo:tensionesvon.jpg|500px|derecha|Tensiones Von Misses sobre la placa]]
{{matlab|codigo=
h=(0.1);
ro=1:h:2;
teta=0:h:2*pi+h;
[roro,tetateta]=meshgrid(ro,teta);
figure(1)
xx=roro.*cos(tetateta);
yy=roro.*sin(tetateta)
f=((-2./roro)-4+2.*roro)/sqrt(2); % El campo de tensiones
surf(xx,yy,f) % Dibujamos la malla
axis([-2,2,-2,2])
colorbar
title('Tensiones de Von Misses')
view(2) }}

==Masa total de la placa==
Si la densidad de la placa es <math> d(x,y,z)=xye^{-1/x^{2}} </math> y el área del anillo es <math> \pi.R^{2}-\pi.r^{2} </math> ,siendo R y r los radios exteriores e interiores respectivamente, la masa total de la placa será el producto del área y la densidad. Al tener la funcíon de la densidad, necesitamos integrarla para calcularla en todo el dominio.
Para una mejor resolución, efectuaremos el cambio a coordenadas polares:
:<math> x=\rho \cdot \cos \theta ; y=\rho \cdot \sin \theta </math>
Tendremos entonces que resolver la integral:
<math>\int_{0 }^{2\pi} \int_{1 }^{2}\rho ^{3}\cos \theta . \sin \theta .\log_{10}(\rho .\cos \theta +2)d\rho d\theta </math>
En matlab sería:
{{matlab|codigo=
h=1/100
u=1:h:2
v=0:h:2*pi
[uu,vv]=meshgrid(u,v)
f=(uu.^3).*sin(vv).*cos(vv).*log(uu.*cos(vv)+2)
a=(h^2)*f
masa=sum(sum(a))
}}
La masa de la placa es <math>m=5,0709.10^{-5}</math>

[[Categoría:Teoría de Campos]]
[[Categoría:TC14/15]]

Visualización de campos escalares y vectoriales en elasticidad. Grupo C4

2014-12-05T12:26:23Z

FernandoSancha:

{{ TrabajoED |Visualización de campos escalares y vectoriales en elasticidad. Grupo C4 | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC14/15|2013-14]] | Alberto Rojas Rivero, Eduardo Bonet García, Fernando Sancha Domínguez, Carlos Fernández Bermejo, Ricardo Mazón Cabrera, Emilio Valero Muñoz-Rojas }}

Consideramos una placa circular comprendida entre las circunferencias de radio 1 y radio 2, centradas en el origen.
En ella vamos a suponer que tenemos definidas dos cantidades físicas: la temperatura <math>T(\rho,\theta,t)</math> ,que depende de las dos variables espaciales <math>(x,y)</math> , y el tiempo <math>t</math>, y los desplazamientos <math>\vec u(\rho,\theta)</math>.De esta forma, si definimos <math>r_0(\rho,\theta)</math> el vector de posición de los puntos de la placa antes de la deformación , la posición de cada punto <math>(\rho,\theta)</math> de la placa después de la deformación viene dada por:
<math>
\vec r (\rho,\theta)= \vec r_{0}(\rho,\theta)+\vec u(\rho,\theta).
</math>
Vamos a suponer que la fuerza aplicada sobre la placa ha provocado una desplazamiento de manera que los desplazamientos en un instante <math>t_0</math> vienen dados por:
<math>\vec u(\rho,\theta)=(1-\rho)^2 \vec g_{\theta}.
</math>

== Mallado del sólido ==
Nuestro mallado consiste en la representación de los puntos interiores de una placa solida, en forma de corona circular. Tomando como paso de muestreo h=1/10, los comandos en Matlab de la visualizacion de nuestra placa es la siguiente:
[[Archivo:Sin título.jpg|350px|miniaturadeimagen|derecha|Mallado de la corona circular.]]
{{matlab|codigo=
h=0.1 % Paso de muestreo
u=1:h:2; % Intervalo [1,2]
v=0:h:2*pi+h; % Intervalo [0,2*pi]
[uu,vv]=meshgrid(u,v); % Matrices de las coordenadas u y v
figure(1)
xx=uu.*cos(vv); % parametrizacion
yy=uu.*sin(vv);
mesh(xx,yy,0*xx) % dibuja la malla
axis([-3,3,-3,3]) % selecciona la región
view(2) }}

== Temperatura del sólido ==
La temperatura es un campo escalar representado por <math>T(x,y)=e^{-y}</math>,vamos a tratar de ver como varía la temperatura en la placa que estamos estudiando . Nuestra función depende únicamente de la variable Y, esto quiere decir que la temperatura se mantendrá constante para todo valor de X. Observamos que cuanto mayor es el valor de la Y ,menor es la temperatura , luego los puntos de máxima temperatura son los más alejados del foco( el origen) y con menor valor de la variable Y, como podemos observar en la figura de la derecha.
[[Archivo:2.jpg|325px|derecha|Distribución de la temperatura en la placa]]
{{matlab|codigo=
h=0.1 % paso de muestreo
u=1:h:2; % Intervalo[1,2]
v=0:h:2*pi+h; % Intervalo [0,2*pi]
[uu,vv]=meshgrid(u,v); % matrices
figure(1)
xx=uu.*cos(vv); % parametrización
yy=uu.*sin(vv);
f=exp(-yy); % Campo escalar
surf(xx,yy,f) % visualización
axis([-3,3,-3,3]) % región de visualización
view(2)}}

==Cálculo y representación del <math>∇T </math>==

Para estudiar como varía la temperatura al movernos por el plano calculamos el gradiente de T.
<math>∇T= -e^{-y} j </math>
En la figura de la derecha podemos observar la representación del campo vectorial <math>∇T </math> y que este es perpendicular a las curvas de nivel del campo escalar T.
[[Archivo:Anillo doblado.jpg|400px|derecha| Curvas de nivel de T y gradiente de T]]
{{matlab|codigo=
g=-exp(-yy);
subplot(1,2,1);
hold on
contour(xx,yy,f,20)
surf(xx,yy,f);
quiver3(xx,yy,f,0*xx,g,0*f);
hold off
subplot(1,2,2)
hold on
quiver(xx,yy,0*xx,Gradiente); colorbar;
contour(xx,yy,f,20)
hold off }}

==Campo de desplazamientos==
Supongamos la aplicación de una fuerza que provoca una vibración sobre la placa, de manera que, en un t aleatorio los desplazamientos vienen dados por:
<math>\vec u(\rho,\theta)=(1-\rho)^2 \vec g_{\theta}.
</math>
Para distinguir el campo de vectores en los puntos del mallado del sólido utilizaremos el código matlab siguiente:
[[Archivo:campodesplazam.jpg|500px|derecha|Distribución de la temperatura en la placCampo de desplazamientos de la placa]]
{{matlab|codigo=
u=1:h:2; % definir intervalo rho [1,2]

v=0:h:2*pi+h; % definir intervalo theta [0,2*pi]

[uu,vv]=meshgrid(u,v); % definir mallado

xx=uu.*cos(vv); % parametrizacion

yy=uu.*sin(vv);

ux=(xx.*(1-uu).^2);

uy=(yy.*(1-uu).^2);

subplot(1,2,1);

mesh(xx,yy,0.*xx)

view(2)

axis([-3,3,-3,3])

hold on

quiver(xx,yy,ux,uy)

title('campo con vectores')

subplot(1,2,2)

plot3(ux,uy,0.*xx)

title('campo de deformaciones')

}}

==Sólido antes y despues de desplazamiento==
El desplazamiento realizado por la placa viene determinado por una serie de movimientos transversales en direcciones distintas. En las imagenes adjuntas calculadas mediante Matlab, se aprecia que el desplazamiento es practicamente nulo.
[[Archivo:antesydespues.jpg|500px|derecha|Placa antes y después del desplazamiento]]
{{matlab|codigo=
h=0.1;
ro=1:h:2; % Valor de ro según el muestreo del intervalo
teta=0:h:2*pi+h; % Valor de teta según el muestreo del intervalo
[roro,tetateta]=meshgrid(ro,teta); %Matrices de la placa en coordenadas polares
figure(1)
xx=roro.*cos(tetateta); % parametrización
yy=roro.*sin(tetateta);
subplot(1,3,1) % figura solido momento inicial
mesh(xx,yy,0*xx) % mallado
axis([-4,4,-4,4]); % ejes
view(2); %visualización
desx=xx+((1-roro).^2).*(-yy); % movimiento parametrizado
desy=yy+((1-roro).^2).*(xx);
subplot(1,3,2) % figura solido en movimiento
mesh(desx,desy,yy.*0) % mallado
axis([-4,4,-4,4]);
view(2);
subplot(1,3,3); % figura solido en movimiento y campo vectorial
hold on
mesh(desx,desy,yy.*0)
axis([-4,4,-4,4]); % ejes
view(2); %
quiver(xx,yy,fx,fy) % campo vectorial del gradiente de temperaturas
axis([-4,4,-4,4]) % ejes
view(2)
hold off

}}

==Divergencia de u==
La divergencia mide el cambio de volumen local, en este caso en nuestra placa, debido al desplazamiento.
Con el grafico adjunto observamos que todos los puntos tienen el mismo valor, cero.
Como demostramos en el apartado anterior, nuestra placa esta sometida a desplazamientos transversales sin variacion de volumen, por lo que no hay cambio del mismo.
[[Archivo:divergence.jpg|350px|derecha|Divergencia en la placa]]
{{matlab|codigo=
h=(0.1);
ro=1:h:2;
teta=0:h:2*pi+h;
[roro,tetateta]=meshgrid(ro,teta);
figure(1)
xx=roro.*cos(tetateta);
yy=roro.*sin(tetateta)
f=(1./roro)-4+3.*roro; % El campo divergencia
surf(xx,yy,f) % Dibujmos la malla
axis([-2,2,-2,2])
colorbar
title('Divergencia del campo del desplazamiento')
view(2)
}}

==Rotacional de u==
A traves del programa de Matlab, que podemos observar a continuacion, hemos calculado el rotacional del campo desplazamiento u en los puntos del solido y los hemos representado, observandose en la grafica cuales son los puntos que sufren mayor rotacional.
[[Archivo:rotacionalc4.jpg|400px|derecha|Rotacional en la placa]]
{{matlab|codigo=
h=(0.1);
ro=1:h:2;
teta=0:h:2*pi+h;
[roro,tetateta]=meshgrid(ro,teta);
figure(1)
xx=roro.*cos(tetateta);
yy=roro.*sin(tetateta);
f=(2+4.*roro.^2-6.*roro);
surf(xx,yy,f)
axis([-2,2,-2,2])
colorbar
title('Rotacional del campo del desplazamiento')
view(2)

}}

== Estudio del tensor de tensiones ==
Definamos <math>\epsilon(\vec u)=(\nabla \vec u + \nabla \vec u^t)/2</math>, la parte simétrica del tensor gradiente de <math>\vec u</math>, que se denomina tensor de deformaciones. En un medio elástico lineal, isótropo y homogéneo los desplazamientos permiten escribir el tensor de tensiones <math>\sigma_{ij}</math> a través de la fórmula::
<math>
\sigma_{ij}=\lambda \nabla \cdot \vec u \delta_{ij} + 2\mu \epsilon_{ij},
</math>
Siendo <math>\nabla \cdot \vec u = \frac{1}{\sqrt{g}}\cdot \ \frac{\partial(\sqrt{g}u^{i})}{\partial x^{i}}=\frac{1}{\rho}\cdot (\frac{\partial \rho - 2\rho ^{2}+ \rho ^{3}}{\partial \rho }+\frac{\partial 0}{\partial \theta })=\frac{1}{\rho }-4+3\rho
</math>

y calculando el gradiente de <math>\vec u</math> como:: <math>\nabla\vec u=u^i_{.j}\vec g_i \otimes \vec g^j</math> donde

<math>u^i_{.j}= \vec g^i \nabla\vec u \vec g_j=\frac{\partial u^i}{\partial x^j}+\Gamma^i_{kj}u^k.= \left(\begin{array}{ccc}-2+2\rho& 0 \\0& \ 0\end{array}\right)+ (1-2\rho +\rho ^{2}).\left(\begin{array}{ccc}0& 0 \\0& \ 1\rho\end{array}\right)
</math>

Podemos calcular el tensor gradiente como 1-contravariante 1-covariante:
:<math>\nabla\vec u=u^i_{.j}=\left(\begin{array}{ccc}-2+2\rho& 0 \\0& \ 1/\rho -2 +\rho\end{array}\right)
</math>

Las coordenadas 1-covariantes 1-contravariante coinciden por lo que la transposicion seá directa y el tensor quedara como <math>(\epsilon^i_{.j})=(\nabla \vec u + \nabla \vec u^t)/2</math> , que corrresponde con la matriz mostrada anteriormente

Sustituyendo ahora los valores en la fórmula del tensor de tensiones,y tomando \(\lambda=\mu=1\), siendo \(\lambda\) y \(\mu\) los conocidos como coeficientes de Lamé, obtenemos:

:<math>\left(\begin{array}{ccc}1/\rho-8+7\rho& 0 \\0& \ 3/\rho -8 +5\rho\end{array}\right)</math>

Las tensiones normales de ambos vectores de la base nos quedaran como:

Las tensiones normales en la dirección \(\vec g_{\rho}\) serían:\[\vec g_{\rho} \cdot \sigma \cdot \vec g_\rho =1/\rho-8+7\rho\]

Las tensiones normales en la dirección \(\vec g_\theta/\rho\) serían:\[\vec g_\theta/\rho \cdot \sigma \cdot \vec g_\theta/\rho = 3/\rho -8 +5\rho\]

A continuación se muestra el código Matlab empleado para dibujar estas tensiones normales en ambas direcciones:
[[Archivo:tensionesnormales.jpg|500px|derecha|Tensiones normales en la placa]]
{{matlab|codigo=
h=0.1; % Intervalo de separación

u=[1:0.1:2]; % Intervalo de rho [1,2]

v=[0:h:2*pi+h]; % Intervalo de theta [0,2*pi]

[uu,vv]=meshgrid(u,v); % Matrices de rho and theta

figure(1)

subplot(1,2,1)

xx=uu.*cos(vv); % Parametrización

yy=uu.*sin(vv);

f=(1./uu)-8+7.*uu; % Campo escalar

surf(xx,yy,f) % Dibujar el mallado

axis([-2,2,-2,2]) % Región seleccionada

axis equal

colorbar

title('Tension normal en dirección g sub rho')

view(2)

subplot(1,2,2)

f=(3./uu)-8+5.*uu; % Campo escalar

surf(xx,yy,f) % Dibujar el mallado

axis([-2,2,-2,2]) % Región del dibujo

axis equal

colorbar

title('Tension normal en dirección g sub theta')

view(2)

}}

Vemos como la tensiones serán mayores a medida que nos alejamos del centro de masas de nuestra placa y simétricas respecto de este ya que dependerá de la distancia de cada uno de los puntos al origen. Además se aprecia que en los bordes interiores de las placas dicha tensión será nula.

===Tensiones tangenciales===
Al realizar el calculo de las tensiones tangenciales de nuestra placa para las direcciones <math>\vec g_\rho</math> y a <math>\vec g_\theta/\rho</math>::
<math>|\sigma \cdot \vec g_\rho-(\vec g_\rho \cdot \sigma \cdot \vec g_\rho) \vec g_\rho| = 0
</math><math>|\sigma \cdot \vec g_\theta/\rho-(\vec g_\theta/\rho \cdot \sigma \cdot \vec g_\theta/\rho) \vec g_\theta/\rho|= 0 </math>

Vemos que ambas tensiones son nulas, algebraicamente hablando esto es debido a que nuestro campo de desplazamientos depende únicamente del valor de la distancia al origen. Esto implica que al estar definido dicho tensor por la divergencia y el tensor unitario, será nula la tensión tangencial siempre que el campo no dependa del angulo, ya que sera igual en todas las direcciones.

==Tensión de Von Mises==
Se define por la fórmula: <math> \sigma _{VM}=\sqrt{\frac{(\sigma 1-\sigma 2)^{2}+(\sigma 2-\sigma 3)^{2}+(\sigma 3-\sigma 1)^{2}}{2}}</math>
En nuestro caso el campo de tensiones de Von Misses viene dado por la fórmula <math> \sigma _{vm}=\frac{1}{\sqrt{2}}\cdot \left ( \frac{-2}{\rho }+2\cdot \rho \right ) </math>.

La representaremos mediante el siguiente programa en MATLAB:
[[Archivo:tensionesvon.jpg|500px|derecha|Distribución de la temperatura en la placa]]
{{matlab|codigo=
h=(0.1);
ro=1:h:2;
teta=0:h:2*pi+h;
[roro,tetateta]=meshgrid(ro,teta);
figure(1)
xx=roro.*cos(tetateta);
yy=roro.*sin(tetateta)
f=((-2./roro)-4+2.*roro)/sqrt(2); % El campo de tensiones
surf(xx,yy,f) % Dibujamos la malla
axis([-2,2,-2,2])
colorbar
title('Tensiones de Von Misses')
view(2) }}

==Masa total de la placa==
Si la densidad de la placa es <math> d(x,y,z)=xye^{-1/x^{2}} </math> y el área del anillo es <math> \pi.R^{2}-\pi.r^{2} </math> ,siendo R y r los radios exteriores e interiores respectivamente, la masa total de la placa será el producto del área y la densidad. Al tener la funcíon de la densidad, necesitamos integrarla para calcularla en todo el dominio.
Para una mejor resolución, efectuaremos el cambio a coordenadas polares:
:<math> x=\rho \cdot \cos \theta ; y=\rho \cdot \sin \theta </math>
Tendremos entonces que resolver la integral:
<math>\int_{0 }^{2\pi} \int_{1 }^{2}\rho ^{3}\cos \theta . \sin \theta .\log_{10}(\rho .\cos \theta +2)d\rho d\theta </math>
En matlab sería:
{{matlab|codigo=
h=1/100
u=1:h:2
v=0:h:2*pi
[uu,vv]=meshgrid(u,v)
f=(uu.^3).*sin(vv).*cos(vv).*log(uu.*cos(vv)+2)
a=(h^2)*f
masa=sum(sum(a))
}}
La masa de la placa es <math>m=5,0709.10^{-5}</math>

[[Categoría:Teoría de Campos]]
[[Categoría:TC14/15]]

Archivo:Tensionesvon.jpg

2014-12-05T12:25:08Z

FernandoSancha:

Visualización de campos escalares y vectoriales en elasticidad. Grupo C4

2014-12-05T12:24:26Z

FernandoSancha:

{{ TrabajoED |Visualización de campos escalares y vectoriales en elasticidad. Grupo C4 | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC14/15|2013-14]] | Alberto Rojas Rivero, Eduardo Bonet García, Fernando Sancha Domínguez, Carlos Fernández Bermejo, Ricardo Mazón Cabrera, Emilio Valero Muñoz-Rojas }}

Consideramos una placa circular comprendida entre las circunferencias de radio 1 y radio 2, centradas en el origen.
En ella vamos a suponer que tenemos definidas dos cantidades físicas: la temperatura <math>T(\rho,\theta,t)</math> ,que depende de las dos variables espaciales <math>(x,y)</math> , y el tiempo <math>t</math>, y los desplazamientos <math>\vec u(\rho,\theta)</math>.De esta forma, si definimos <math>r_0(\rho,\theta)</math> el vector de posición de los puntos de la placa antes de la deformación , la posición de cada punto <math>(\rho,\theta)</math> de la placa después de la deformación viene dada por:
<math>
\vec r (\rho,\theta)= \vec r_{0}(\rho,\theta)+\vec u(\rho,\theta).
</math>
Vamos a suponer que la fuerza aplicada sobre la placa ha provocado una desplazamiento de manera que los desplazamientos en un instante <math>t_0</math> vienen dados por:
<math>\vec u(\rho,\theta)=(1-\rho)^2 \vec g_{\theta}.
</math>

== Mallado del sólido ==
Nuestro mallado consiste en la representación de los puntos interiores de una placa solida, en forma de corona circular. Tomando como paso de muestreo h=1/10, los comandos en Matlab de la visualizacion de nuestra placa es la siguiente:
[[Archivo:Sin título.jpg|350px|miniaturadeimagen|derecha|Mallado de la corona circular.]]
{{matlab|codigo=
h=0.1 % Paso de muestreo
u=1:h:2; % Intervalo [1,2]
v=0:h:2*pi+h; % Intervalo [0,2*pi]
[uu,vv]=meshgrid(u,v); % Matrices de las coordenadas u y v
figure(1)
xx=uu.*cos(vv); % parametrizacion
yy=uu.*sin(vv);
mesh(xx,yy,0*xx) % dibuja la malla
axis([-3,3,-3,3]) % selecciona la región
view(2) }}

== Temperatura del sólido ==
La temperatura es un campo escalar representado por <math>T(x,y)=e^{-y}</math>,vamos a tratar de ver como varía la temperatura en la placa que estamos estudiando . Nuestra función depende únicamente de la variable Y, esto quiere decir que la temperatura se mantendrá constante para todo valor de X. Observamos que cuanto mayor es el valor de la Y ,menor es la temperatura , luego los puntos de máxima temperatura son los más alejados del foco( el origen) y con menor valor de la variable Y, como podemos observar en la figura de la derecha.
[[Archivo:2.jpg|325px|derecha|Distribución de la temperatura en la placa]]
{{matlab|codigo=
h=0.1 % paso de muestreo
u=1:h:2; % Intervalo[1,2]
v=0:h:2*pi+h; % Intervalo [0,2*pi]
[uu,vv]=meshgrid(u,v); % matrices
figure(1)
xx=uu.*cos(vv); % parametrización
yy=uu.*sin(vv);
f=exp(-yy); % Campo escalar
surf(xx,yy,f) % visualización
axis([-3,3,-3,3]) % región de visualización
view(2)}}

==Cálculo y representación del <math>∇T </math>==

Para estudiar como varía la temperatura al movernos por el plano calculamos el gradiente de T.
<math>∇T= -e^{-y} j </math>
En la figura de la derecha podemos observar la representación del campo vectorial <math>∇T </math> y que este es perpendicular a las curvas de nivel del campo escalar T.
[[Archivo:Anillo doblado.jpg|400px|derecha| Curvas de nivel de T y gradiente de T]]
{{matlab|codigo=
g=-exp(-yy);
subplot(1,2,1);
hold on
contour(xx,yy,f,20)
surf(xx,yy,f);
quiver3(xx,yy,f,0*xx,g,0*f);
hold off
subplot(1,2,2)
hold on
quiver(xx,yy,0*xx,Gradiente); colorbar;
contour(xx,yy,f,20)
hold off }}

==Campo de desplazamientos==
Supongamos la aplicación de una fuerza que provoca una vibración sobre la placa, de manera que, en un t aleatorio los desplazamientos vienen dados por:
<math>\vec u(\rho,\theta)=(1-\rho)^2 \vec g_{\theta}.
</math>
Para distinguir el campo de vectores en los puntos del mallado del sólido utilizaremos el código matlab siguiente:
[[Archivo:campodesplazam.jpg|500px|derecha|Distribución de la temperatura en la placCampo de desplazamientos de la placa]]
{{matlab|codigo=
u=1:h:2; % definir intervalo rho [1,2]

v=0:h:2*pi+h; % definir intervalo theta [0,2*pi]

[uu,vv]=meshgrid(u,v); % definir mallado

xx=uu.*cos(vv); % parametrizacion

yy=uu.*sin(vv);

ux=(xx.*(1-uu).^2);

uy=(yy.*(1-uu).^2);

subplot(1,2,1);

mesh(xx,yy,0.*xx)

view(2)

axis([-3,3,-3,3])

hold on

quiver(xx,yy,ux,uy)

title('campo con vectores')

subplot(1,2,2)

plot3(ux,uy,0.*xx)

title('campo de deformaciones')

}}

==Sólido antes y despues de desplazamiento==
El desplazamiento realizado por la placa viene determinado por una serie de movimientos transversales en direcciones distintas. En las imagenes adjuntas calculadas mediante Matlab, se aprecia que el desplazamiento es practicamente nulo.
[[Archivo:antesydespues.jpg|500px|derecha|Placa antes y después del desplazamiento]]
{{matlab|codigo=
h=0.1;
ro=1:h:2; % Valor de ro según el muestreo del intervalo
teta=0:h:2*pi+h; % Valor de teta según el muestreo del intervalo
[roro,tetateta]=meshgrid(ro,teta); %Matrices de la placa en coordenadas polares
figure(1)
xx=roro.*cos(tetateta); % parametrización
yy=roro.*sin(tetateta);
subplot(1,3,1) % figura solido momento inicial
mesh(xx,yy,0*xx) % mallado
axis([-4,4,-4,4]); % ejes
view(2); %visualización
desx=xx+((1-roro).^2).*(-yy); % movimiento parametrizado
desy=yy+((1-roro).^2).*(xx);
subplot(1,3,2) % figura solido en movimiento
mesh(desx,desy,yy.*0) % mallado
axis([-4,4,-4,4]);
view(2);
subplot(1,3,3); % figura solido en movimiento y campo vectorial
hold on
mesh(desx,desy,yy.*0)
axis([-4,4,-4,4]); % ejes
view(2); %
quiver(xx,yy,fx,fy) % campo vectorial del gradiente de temperaturas
axis([-4,4,-4,4]) % ejes
view(2)
hold off

}}

==Divergencia de u==
La divergencia mide el cambio de volumen local, en este caso en nuestra placa, debido al desplazamiento.
Con el grafico adjunto observamos que todos los puntos tienen el mismo valor, cero.
Como demostramos en el apartado anterior, nuestra placa esta sometida a desplazamientos transversales sin variacion de volumen, por lo que no hay cambio del mismo.
[[Archivo:divergence.jpg|350px|derecha|Divergencia en la placa]]
{{matlab|codigo=
h=(0.1);
ro=1:h:2;
teta=0:h:2*pi+h;
[roro,tetateta]=meshgrid(ro,teta);
figure(1)
xx=roro.*cos(tetateta);
yy=roro.*sin(tetateta)
f=(1./roro)-4+3.*roro; % El campo divergencia
surf(xx,yy,f) % Dibujmos la malla
axis([-2,2,-2,2])
colorbar
title('Divergencia del campo del desplazamiento')
view(2)
}}

==Rotacional de u==
A traves del programa de Matlab, que podemos observar a continuacion, hemos calculado el rotacional del campo desplazamiento u en los puntos del solido y los hemos representado, observandose en la grafica cuales son los puntos que sufren mayor rotacional.
[[Archivo:rotacionalc4.jpg|400px|derecha|Rotacional en la placa]]
{{matlab|codigo=
h=(0.1);
ro=1:h:2;
teta=0:h:2*pi+h;
[roro,tetateta]=meshgrid(ro,teta);
figure(1)
xx=roro.*cos(tetateta);
yy=roro.*sin(tetateta);
f=(2+4.*roro.^2-6.*roro);
surf(xx,yy,f)
axis([-2,2,-2,2])
colorbar
title('Rotacional del campo del desplazamiento')
view(2)

}}

== Estudio del tensor de tensiones ==
Definamos <math>\epsilon(\vec u)=(\nabla \vec u + \nabla \vec u^t)/2</math>, la parte simétrica del tensor gradiente de <math>\vec u</math>, que se denomina tensor de deformaciones. En un medio elástico lineal, isótropo y homogéneo los desplazamientos permiten escribir el tensor de tensiones <math>\sigma_{ij}</math> a través de la fórmula::
<math>
\sigma_{ij}=\lambda \nabla \cdot \vec u \delta_{ij} + 2\mu \epsilon_{ij},
</math>
Siendo <math>\nabla \cdot \vec u = \frac{1}{\sqrt{g}}\cdot \ \frac{\partial(\sqrt{g}u^{i})}{\partial x^{i}}=\frac{1}{\rho}\cdot (\frac{\partial \rho - 2\rho ^{2}+ \rho ^{3}}{\partial \rho }+\frac{\partial 0}{\partial \theta })=\frac{1}{\rho }-4+3\rho
</math>

y calculando el gradiente de <math>\vec u</math> como:: <math>\nabla\vec u=u^i_{.j}\vec g_i \otimes \vec g^j</math> donde

<math>u^i_{.j}= \vec g^i \nabla\vec u \vec g_j=\frac{\partial u^i}{\partial x^j}+\Gamma^i_{kj}u^k.= \left(\begin{array}{ccc}-2+2\rho& 0 \\0& \ 0\end{array}\right)+ (1-2\rho +\rho ^{2}).\left(\begin{array}{ccc}0& 0 \\0& \ 1\rho\end{array}\right)
</math>

Podemos calcular el tensor gradiente como 1-contravariante 1-covariante:
:<math>\nabla\vec u=u^i_{.j}=\left(\begin{array}{ccc}-2+2\rho& 0 \\0& \ 1/\rho -2 +\rho\end{array}\right)
</math>

Las coordenadas 1-covariantes 1-contravariante coinciden por lo que la transposicion seá directa y el tensor quedara como <math>(\epsilon^i_{.j})=(\nabla \vec u + \nabla \vec u^t)/2</math> , que corrresponde con la matriz mostrada anteriormente

Sustituyendo ahora los valores en la fórmula del tensor de tensiones,y tomando \(\lambda=\mu=1\), siendo \(\lambda\) y \(\mu\) los conocidos como coeficientes de Lamé, obtenemos:

:<math>\left(\begin{array}{ccc}1/\rho-8+7\rho& 0 \\0& \ 3/\rho -8 +5\rho\end{array}\right)</math>

Las tensiones normales de ambos vectores de la base nos quedaran como:

Las tensiones normales en la dirección \(\vec g_{\rho}\) serían:\[\vec g_{\rho} \cdot \sigma \cdot \vec g_\rho =1/\rho-8+7\rho\]

Las tensiones normales en la dirección \(\vec g_\theta/\rho\) serían:\[\vec g_\theta/\rho \cdot \sigma \cdot \vec g_\theta/\rho = 3/\rho -8 +5\rho\]

A continuación se muestra el código Matlab empleado para dibujar estas tensiones normales en ambas direcciones:
[[Archivo:tensionesnormales.jpg|500px|derecha|Tensiones normales en la placa]]
{{matlab|codigo=
h=0.1; % Intervalo de separación

u=[1:0.1:2]; % Intervalo de rho [1,2]

v=[0:h:2*pi+h]; % Intervalo de theta [0,2*pi]

[uu,vv]=meshgrid(u,v); % Matrices de rho and theta

figure(1)

subplot(1,2,1)

xx=uu.*cos(vv); % Parametrización

yy=uu.*sin(vv);

f=(1./uu)-8+7.*uu; % Campo escalar

surf(xx,yy,f) % Dibujar el mallado

axis([-2,2,-2,2]) % Región seleccionada

axis equal

colorbar

title('Tension normal en dirección g sub rho')

view(2)

subplot(1,2,2)

f=(3./uu)-8+5.*uu; % Campo escalar

surf(xx,yy,f) % Dibujar el mallado

axis([-2,2,-2,2]) % Región del dibujo

axis equal

colorbar

title('Tension normal en dirección g sub theta')

view(2)

}}

Vemos como la tensiones serán mayores a medida que nos alejamos del centro de masas de nuestra placa y simétricas respecto de este ya que dependerá de la distancia de cada uno de los puntos al origen. Además se aprecia que en los bordes interiores de las placas dicha tensión será nula.

===Tensiones tangenciales===
Al realizar el calculo de las tensiones tangenciales de nuestra placa para las direcciones <math>\vec g_\rho</math> y a <math>\vec g_\theta/\rho</math>::
<math>|\sigma \cdot \vec g_\rho-(\vec g_\rho \cdot \sigma \cdot \vec g_\rho) \vec g_\rho| = 0
</math><math>|\sigma \cdot \vec g_\theta/\rho-(\vec g_\theta/\rho \cdot \sigma \cdot \vec g_\theta/\rho) \vec g_\theta/\rho|= 0 </math>

Vemos que ambas tensiones son nulas, algebraicamente hablando esto es debido a que nuestro campo de desplazamientos depende únicamente del valor de la distancia al origen. Esto implica que al estar definido dicho tensor por la divergencia y el tensor unitario, será nula la tensión tangencial siempre que el campo no dependa del angulo, ya que sera igual en todas las direcciones.

==Tensión de Von Mises==
Se define por la fórmula: <math> \sigma _{VM}=\sqrt{\frac{(\sigma 1-\sigma 2)^{2}+(\sigma 2-\sigma 3)^{2}+(\sigma 3-\sigma 1)^{2}}{2}}</math>
En nuestro caso el campo de tensiones de Von Misses viene dado por la fórmula <math> \sigma _{vm}=\frac{1}{\sqrt{2}}\cdot \left ( \frac{-2}{\rho }+2\cdot \rho \right ) </math>.
La representaremos mediante el siguiente programa en MATLAB:
{{matlab|codigo=
h=(0.1);
ro=1:h:2;
teta=0:h:2*pi+h;
[roro,tetateta]=meshgrid(ro,teta);
figure(1)
xx=roro.*cos(tetateta);
yy=roro.*sin(tetateta)
f=((-2./roro)-4+2.*roro)/sqrt(2); % El campo de tensiones
surf(xx,yy,f) % Dibujamos la malla
axis([-2,2,-2,2])
colorbar
title('Tensiones de Von Misses')
view(2) }}

==Masa total de la placa==
Si la densidad de la placa es <math> d(x,y,z)=xye^{-1/x^{2}} </math> y el área del anillo es <math> \pi.R^{2}-\pi.r^{2} </math> ,siendo R y r los radios exteriores e interiores respectivamente, la masa total de la placa será el producto del área y la densidad. Al tener la funcíon de la densidad, necesitamos integrarla para calcularla en todo el dominio.
Para una mejor resolución, efectuaremos el cambio a coordenadas polares:
:<math> x=\rho \cdot \cos \theta ; y=\rho \cdot \sin \theta </math>
Tendremos entonces que resolver la integral:
<math>\int_{0 }^{2\pi} \int_{1 }^{2}\rho ^{3}\cos \theta . \sin \theta .\log_{10}(\rho .\cos \theta +2)d\rho d\theta </math>
En matlab sería:
{{matlab|codigo=
h=1/100
u=1:h:2
v=0:h:2*pi
[uu,vv]=meshgrid(u,v)
f=(uu.^3).*sin(vv).*cos(vv).*log(uu.*cos(vv)+2)
a=(h^2)*f
masa=sum(sum(a))
}}
La masa de la placa es <math>m=5,0709.10^{-5}</math>

[[Categoría:Teoría de Campos]]
[[Categoría:TC14/15]]

Visualización de campos escalares y vectoriales en elasticidad. Grupo C4

2014-12-05T12:21:22Z

FernandoSancha: /* Tensión de Von Mises */

{{ TrabajoED |Visualización de campos escalares y vectoriales en elasticidad. Grupo C4 | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC14/15|2013-14]] | Alberto Rojas Rivero, Eduardo Bonet García, Fernando Sancha Domínguez, Carlos Fernández Bermejo, Ricardo Mazón Cabrera, Emilio Valero Muñoz-Rojas }}

Consideramos una placa circular comprendida entre las circunferencias de radio 1 y radio 2, centradas en el origen.
En ella vamos a suponer que tenemos definidas dos cantidades físicas: la temperatura <math>T(\rho,\theta,t)</math> ,que depende de las dos variables espaciales <math>(x,y)</math> , y el tiempo <math>t</math>, y los desplazamientos <math>\vec u(\rho,\theta)</math>.De esta forma, si definimos <math>r_0(\rho,\theta)</math> el vector de posición de los puntos de la placa antes de la deformación , la posición de cada punto <math>(\rho,\theta)</math> de la placa después de la deformación viene dada por:
<math>
\vec r (\rho,\theta)= \vec r_{0}(\rho,\theta)+\vec u(\rho,\theta).
</math>
Vamos a suponer que la fuerza aplicada sobre la placa ha provocado una desplazamiento de manera que los desplazamientos en un instante <math>t_0</math> vienen dados por:
<math>\vec u(\rho,\theta)=(1-\rho)^2 \vec g_{\theta}.
</math>

== Mallado del sólido ==
Nuestro mallado consiste en la representación de los puntos interiores de una placa solida, en forma de corona circular. Tomando como paso de muestreo h=1/10, los comandos en Matlab de la visualizacion de nuestra placa es la siguiente:
[[Archivo:Sin título.jpg|350px|miniaturadeimagen|derecha|Mallado de la corona circular.]]
{{matlab|codigo=
h=0.1 % Paso de muestreo
u=1:h:2; % Intervalo [1,2]
v=0:h:2*pi+h; % Intervalo [0,2*pi]
[uu,vv]=meshgrid(u,v); % Matrices de las coordenadas u y v
figure(1)
xx=uu.*cos(vv); % parametrizacion
yy=uu.*sin(vv);
mesh(xx,yy,0*xx) % dibuja la malla
axis([-3,3,-3,3]) % selecciona la región
view(2) }}

== Temperatura del sólido ==
La temperatura es un campo escalar representado por <math>T(x,y)=e^{-y}</math>,vamos a tratar de ver como varía la temperatura en la placa que estamos estudiando . Nuestra función depende únicamente de la variable Y, esto quiere decir que la temperatura se mantendrá constante para todo valor de X. Observamos que cuanto mayor es el valor de la Y ,menor es la temperatura , luego los puntos de máxima temperatura son los más alejados del foco( el origen) y con menor valor de la variable Y, como podemos observar en la figura de la derecha.
[[Archivo:2.jpg|325px|derecha|Distribución de la temperatura en la placa]]
{{matlab|codigo=
h=0.1 % paso de muestreo
u=1:h:2; % Intervalo[1,2]
v=0:h:2*pi+h; % Intervalo [0,2*pi]
[uu,vv]=meshgrid(u,v); % matrices
figure(1)
xx=uu.*cos(vv); % parametrización
yy=uu.*sin(vv);
f=exp(-yy); % Campo escalar
surf(xx,yy,f) % visualización
axis([-3,3,-3,3]) % región de visualización
view(2)}}

==Cálculo y representación del <math>∇T </math>==

Para estudiar como varía la temperatura al movernos por el plano calculamos el gradiente de T.
<math>∇T= -e^{-y} j </math>
En la figura de la derecha podemos observar la representación del campo vectorial <math>∇T </math> y que este es perpendicular a las curvas de nivel del campo escalar T.
[[Archivo:Anillo doblado.jpg|400px|derecha| Curvas de nivel de T y gradiente de T]]
{{matlab|codigo=
g=-exp(-yy);
subplot(1,2,1);
hold on
contour(xx,yy,f,20)
surf(xx,yy,f);
quiver3(xx,yy,f,0*xx,g,0*f);
hold off
subplot(1,2,2)
hold on
quiver(xx,yy,0*xx,Gradiente); colorbar;
contour(xx,yy,f,20)
hold off }}

==Campo de desplazamientos==
Supongamos la aplicación de una fuerza que provoca una vibración sobre la placa, de manera que, en un t aleatorio los desplazamientos vienen dados por:
<math>\vec u(\rho,\theta)=(1-\rho)^2 \vec g_{\theta}.
</math>
Para distinguir el campo de vectores en los puntos del mallado del sólido utilizaremos el código matlab siguiente:
[[Archivo:campodesplazam.jpg|500px|derecha|Distribución de la temperatura en la placCampo de desplazamientos de la placa]]
{{matlab|codigo=
u=1:h:2; % definir intervalo rho [1,2]

v=0:h:2*pi+h; % definir intervalo theta [0,2*pi]

[uu,vv]=meshgrid(u,v); % definir mallado

xx=uu.*cos(vv); % parametrizacion

yy=uu.*sin(vv);

ux=(xx.*(1-uu).^2);

uy=(yy.*(1-uu).^2);

subplot(1,2,1);

mesh(xx,yy,0.*xx)

view(2)

axis([-3,3,-3,3])

hold on

quiver(xx,yy,ux,uy)

title('campo con vectores')

subplot(1,2,2)

plot3(ux,uy,0.*xx)

title('campo de deformaciones')

}}

==Sólido antes y despues de desplazamiento==
El desplazamiento realizado por la placa viene determinado por una serie de movimientos transversales en direcciones distintas. En las imagenes adjuntas calculadas mediante Matlab, se aprecia que el desplazamiento es practicamente nulo.
[[Archivo:antesydespues.jpg|500px|derecha|Placa antes y después del desplazamiento]]
{{matlab|codigo=
h=0.1;
ro=1:h:2; % Valor de ro según el muestreo del intervalo
teta=0:h:2*pi+h; % Valor de teta según el muestreo del intervalo
[roro,tetateta]=meshgrid(ro,teta); %Matrices de la placa en coordenadas polares
figure(1)
xx=roro.*cos(tetateta); % parametrización
yy=roro.*sin(tetateta);
subplot(1,3,1) % figura solido momento inicial
mesh(xx,yy,0*xx) % mallado
axis([-4,4,-4,4]); % ejes
view(2); %visualización
desx=xx+((1-roro).^2).*(-yy); % movimiento parametrizado
desy=yy+((1-roro).^2).*(xx);
subplot(1,3,2) % figura solido en movimiento
mesh(desx,desy,yy.*0) % mallado
axis([-4,4,-4,4]);
view(2);
subplot(1,3,3); % figura solido en movimiento y campo vectorial
hold on
mesh(desx,desy,yy.*0)
axis([-4,4,-4,4]); % ejes
view(2); %
quiver(xx,yy,fx,fy) % campo vectorial del gradiente de temperaturas
axis([-4,4,-4,4]) % ejes
view(2)
hold off

}}

==Divergencia de u==
La divergencia mide el cambio de volumen local, en este caso en nuestra placa, debido al desplazamiento.
Con el grafico adjunto observamos que todos los puntos tienen el mismo valor, cero.
Como demostramos en el apartado anterior, nuestra placa esta sometida a desplazamientos transversales sin variacion de volumen, por lo que no hay cambio del mismo.
[[Archivo:divergence.jpg|350px|derecha|Divergencia en la placa]]
{{matlab|codigo=
h=(0.1);
ro=1:h:2;
teta=0:h:2*pi+h;
[roro,tetateta]=meshgrid(ro,teta);
figure(1)
xx=roro.*cos(tetateta);
yy=roro.*sin(tetateta)
f=(1./roro)-4+3.*roro; % El campo divergencia
surf(xx,yy,f) % Dibujmos la malla
axis([-2,2,-2,2])
colorbar
title('Divergencia del campo del desplazamiento')
view(2)
}}

==Rotacional de u==
A traves del programa de Matlab, que podemos observar a continuacion, hemos calculado el rotacional del campo desplazamiento u en los puntos del solido y los hemos representado, observandose en la grafica cuales son los puntos que sufren mayor rotacional.
[[Archivo:rotacionalc4.jpg|400px|derecha|Rotacional en la placa]]
{{matlab|codigo=
h=(0.1);
ro=1:h:2;
teta=0:h:2*pi+h;
[roro,tetateta]=meshgrid(ro,teta);
figure(1)
xx=roro.*cos(tetateta);
yy=roro.*sin(tetateta);
f=(2+4.*roro.^2-6.*roro);
surf(xx,yy,f)
axis([-2,2,-2,2])
colorbar
title('Rotacional del campo del desplazamiento')
view(2)

}}

== Estudio del tensor de tensiones ==
Definamos <math>\epsilon(\vec u)=(\nabla \vec u + \nabla \vec u^t)/2</math>, la parte simétrica del tensor gradiente de <math>\vec u</math>, que se denomina tensor de deformaciones. En un medio elástico lineal, isótropo y homogéneo los desplazamientos permiten escribir el tensor de tensiones <math>\sigma_{ij}</math> a través de la fórmula::
<math>
\sigma_{ij}=\lambda \nabla \cdot \vec u \delta_{ij} + 2\mu \epsilon_{ij},
</math>
Siendo <math>\nabla \cdot \vec u = \frac{1}{\sqrt{g}}\cdot \ \frac{\partial(\sqrt{g}u^{i})}{\partial x^{i}}=\frac{1}{\rho}\cdot (\frac{\partial \rho - 2\rho ^{2}+ \rho ^{3}}{\partial \rho }+\frac{\partial 0}{\partial \theta })=\frac{1}{\rho }-4+3\rho
</math>

y calculando el gradiente de <math>\vec u</math> como:: <math>\nabla\vec u=u^i_{.j}\vec g_i \otimes \vec g^j</math> donde

<math>u^i_{.j}= \vec g^i \nabla\vec u \vec g_j=\frac{\partial u^i}{\partial x^j}+\Gamma^i_{kj}u^k.= \left(\begin{array}{ccc}-2+2\rho& 0 \\0& \ 0\end{array}\right)+ (1-2\rho +\rho ^{2}).\left(\begin{array}{ccc}0& 0 \\0& \ 1\rho\end{array}\right)
</math>

Podemos calcular el tensor gradiente como 1-contravariante 1-covariante:
:<math>\nabla\vec u=u^i_{.j}=\left(\begin{array}{ccc}-2+2\rho& 0 \\0& \ 1/\rho -2 +\rho\end{array}\right)
</math>

Las coordenadas 1-covariantes 1-contravariante coinciden por lo que la transposicion seá directa y el tensor quedara como <math>(\epsilon^i_{.j})=(\nabla \vec u + \nabla \vec u^t)/2</math> , que corrresponde con la matriz mostrada anteriormente

Sustituyendo ahora los valores en la fórmula del tensor de tensiones,y tomando \(\lambda=\mu=1\), siendo \(\lambda\) y \(\mu\) los conocidos como coeficientes de Lamé, obtenemos:

:<math>\left(\begin{array}{ccc}1/\rho-8+7\rho& 0 \\0& \ 3/\rho -8 +5\rho\end{array}\right)</math>

Las tensiones normales de ambos vectores de la base nos quedaran como:

Las tensiones normales en la dirección \(\vec g_{\rho}\) serían:\[\vec g_{\rho} \cdot \sigma \cdot \vec g_\rho =1/\rho-8+7\rho\]

Las tensiones normales en la dirección \(\vec g_\theta/\rho\) serían:\[\vec g_\theta/\rho \cdot \sigma \cdot \vec g_\theta/\rho = 3/\rho -8 +5\rho\]

A continuación se muestra el código Matlab empleado para dibujar estas tensiones normales en ambas direcciones:
[[Archivo:tensionesnormales.jpg|500px|derecha|Tensiones normales en la placa]]
{{matlab|codigo=
h=0.1; % Intervalo de separación

u=[1:0.1:2]; % Intervalo de rho [1,2]

v=[0:h:2*pi+h]; % Intervalo de theta [0,2*pi]

[uu,vv]=meshgrid(u,v); % Matrices de rho and theta

figure(1)

subplot(1,2,1)

xx=uu.*cos(vv); % Parametrización

yy=uu.*sin(vv);

f=(1./uu)-8+7.*uu; % Campo escalar

surf(xx,yy,f) % Dibujar el mallado

axis([-2,2,-2,2]) % Región seleccionada

axis equal

colorbar

title('Tension normal en dirección g sub rho')

view(2)

subplot(1,2,2)

f=(3./uu)-8+5.*uu; % Campo escalar

surf(xx,yy,f) % Dibujar el mallado

axis([-2,2,-2,2]) % Región del dibujo

axis equal

colorbar

title('Tension normal en dirección g sub theta')

view(2)

}}

Vemos como la tensiones serán mayores a medida que nos alejamos del centro de masas de nuestra placa y simétricas respecto de este ya que dependerá de la distancia de cada uno de los puntos al origen. Además se aprecia que en los bordes interiores de las placas dicha tensión será nula.

===Tensiones tangenciales===
Al realizar el calculo de las tensiones tangenciales de nuestra placa para las direcciones <math>\vec g_\rho</math> y a <math>\vec g_\theta/\rho</math>::
<math>|\sigma \cdot \vec g_\rho-(\vec g_\rho \cdot \sigma \cdot \vec g_\rho) \vec g_\rho| = 0
</math><math>|\sigma \cdot \vec g_\theta/\rho-(\vec g_\theta/\rho \cdot \sigma \cdot \vec g_\theta/\rho) \vec g_\theta/\rho|= 0 </math>

Vemos que ambas tensiones son nulas, algebraicamente hablando esto es debido a que nuestro campo de desplazamientos depende únicamente del valor de la distancia al origen. Esto implica que al estar definido dicho tensor por la divergencia y el tensor unitario, será nula la tensión tangencial siempre que el campo no dependa del angulo, ya que sera igual en todas las direcciones.

==Tensión de Von Mises==
Se define por la fórmula: <math> \sigma _{VM}=\sqrt{\frac{(\sigma 1-\sigma 2)^{2}+(\sigma 2-\sigma 3)^{2}+(\sigma 3-\sigma 1)^{2}}{2}}</math>
En nuestro caso el campo de tensiones de Von Misses viene dado por la fórmula <math> \sigma _{vm}=\frac{1}{\sqrt{2}}\cdot \left ( \frac{-2}{\rho }+2\cdot \rho \right ) </math>
La representaremos mediante el siguiente programa en MATLAB:

==Masa total de la placa==
Si la densidad de la placa es <math> d(x,y,z)=xye^{-1/x^{2}} </math> y el área del anillo es <math> \pi.R^{2}-\pi.r^{2} </math> ,siendo R y r los radios exteriores e interiores respectivamente, la masa total de la placa será el producto del área y la densidad. Al tener la funcíon de la densidad, necesitamos integrarla para calcularla en todo el dominio.
Para una mejor resolución, efectuaremos el cambio a coordenadas polares:
:<math> x=\rho \cdot \cos \theta ; y=\rho \cdot \sin \theta </math>
Tendremos entonces que resolver la integral:
<math>\int_{0 }^{2\pi} \int_{1 }^{2}\rho ^{3}\cos \theta . \sin \theta .\log_{10}(\rho .\cos \theta +2)d\rho d\theta </math>
En matlab sería:
{{matlab|codigo=
h=1/100
u=1:h:2
v=0:h:2*pi
[uu,vv]=meshgrid(u,v)
f=(uu.^3).*sin(vv).*cos(vv).*log(uu.*cos(vv)+2)
a=(h^2)*f
masa=sum(sum(a))
}}
La masa de la placa es <math>m=5,0709.10^{-5}</math>

[[Categoría:Teoría de Campos]]
[[Categoría:TC14/15]]

Visualización de campos escalares y vectoriales en elasticidad. Grupo C4

2014-12-05T11:54:52Z

FernandoSancha: /* Divergencia de u */

{{ TrabajoED |Visualización de campos escalares y vectoriales en elasticidad. Grupo C4 | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC14/15|2013-14]] | Alberto Rojas Rivero, Eduardo Bonet García, Fernando Sancha Domínguez, Carlos Fernández Bermejo, Ricardo Mazón Cabrera, Emilio Valero Muñoz-Rojas }}

Consideramos una placa circular comprendida entre las circunferencias de radio 1 y radio 2, centradas en el origen.
En ella vamos a suponer que tenemos definidas dos cantidades físicas: la temperatura <math>T(\rho,\theta,t)</math> ,que depende de las dos variables espaciales <math>(x,y)</math> , y el tiempo <math>t</math>, y los desplazamientos <math>\vec u(\rho,\theta)</math>.De esta forma, si definimos <math>r_0(\rho,\theta)</math> el vector de posición de los puntos de la placa antes de la deformación , la posición de cada punto <math>(\rho,\theta)</math> de la placa después de la deformación viene dada por:
<math>
\vec r (\rho,\theta)= \vec r_{0}(\rho,\theta)+\vec u(\rho,\theta).
</math>
Vamos a suponer que la fuerza aplicada sobre la placa ha provocado una desplazamiento de manera que los desplazamientos en un instante <math>t_0</math> vienen dados por:
<math>\vec u(\rho,\theta)=(1-\rho)^2 \vec g_{\theta}.
</math>

== Mallado del sólido ==
Nuestro mallado consiste en la representación de los puntos interiores de una placa solida, en forma de corona circular. Tomando como paso de muestreo h=1/10, los comandos en Matlab de la visualizacion de nuestra placa es la siguiente:
[[Archivo:Sin título.jpg|350px|miniaturadeimagen|derecha|Mallado de la corona circular.]]
{{matlab|codigo=
h=0.1 % Paso de muestreo
u=1:h:2; % Intervalo [1,2]
v=0:h:2*pi+h; % Intervalo [0,2*pi]
[uu,vv]=meshgrid(u,v); % Matrices de las coordenadas u y v
figure(1)
xx=uu.*cos(vv); % parametrizacion
yy=uu.*sin(vv);
mesh(xx,yy,0*xx) % dibuja la malla
axis([-3,3,-3,3]) % selecciona la región
view(2) }}

== Temperatura del sólido ==
La temperatura es un campo escalar representado por <math>T(x,y)=e^{-y}</math>,vamos a tratar de ver como varía la temperatura en la placa que estamos estudiando . Nuestra función depende únicamente de la variable Y, esto quiere decir que la temperatura se mantendrá constante para todo valor de X. Observamos que cuanto mayor es el valor de la Y ,menor es la temperatura , luego los puntos de máxima temperatura son los más alejados del foco( el origen) y con menor valor de la variable Y, como podemos observar en la figura de la derecha.
[[Archivo:2.jpg|325px|derecha|Distribución de la temperatura en la placa]]
{{matlab|codigo=
h=0.1 % paso de muestreo
u=1:h:2; % Intervalo[1,2]
v=0:h:2*pi+h; % Intervalo [0,2*pi]
[uu,vv]=meshgrid(u,v); % matrices
figure(1)
xx=uu.*cos(vv); % parametrización
yy=uu.*sin(vv);
f=exp(-yy); % Campo escalar
surf(xx,yy,f) % visualización
axis([-3,3,-3,3]) % región de visualización
view(2)}}

==Cálculo y representación del <math>∇T </math>==

Para estudiar como varía la temperatura al movernos por el plano calculamos el gradiente de T.
<math>∇T= -e^{-y} j </math>
En la figura de la derecha podemos observar la representación del campo vectorial <math>∇T </math> y que este es perpendicular a las curvas de nivel del campo escalar T.
[[Archivo:Anillo doblado.jpg|400px|derecha| Curvas de nivel de T y gradiente de T]]
{{matlab|codigo=
g=-exp(-yy);
subplot(1,2,1);
hold on
contour(xx,yy,f,20)
surf(xx,yy,f);
quiver3(xx,yy,f,0*xx,g,0*f);
hold off
subplot(1,2,2)
hold on
quiver(xx,yy,0*xx,Gradiente); colorbar;
contour(xx,yy,f,20)
hold off }}

==Campo de desplazamientos==
Supongamos la aplicación de una fuerza que provoca una vibración sobre la placa, de manera que, en un t aleatorio los desplazamientos vienen dados por:
<math>\vec u(\rho,\theta)=(1-\rho)^2 \vec g_{\theta}.
</math>
Para distinguir el campo de vectores en los puntos del mallado del sólido utilizaremos el código matlab siguiente:
[[Archivo:campodesplazam.jpg|500px|derecha|Distribución de la temperatura en la placCampo de desplazamientos de la placa]]
{{matlab|codigo=
u=1:h:2; % definir intervalo rho [1,2]

v=0:h:2*pi+h; % definir intervalo theta [0,2*pi]

[uu,vv]=meshgrid(u,v); % definir mallado

xx=uu.*cos(vv); % parametrizacion

yy=uu.*sin(vv);

ux=(xx.*(1-uu).^2);

uy=(yy.*(1-uu).^2);

subplot(1,2,1);

mesh(xx,yy,0.*xx)

view(2)

axis([-3,3,-3,3])

hold on

quiver(xx,yy,ux,uy)

title('campo con vectores')

subplot(1,2,2)

plot3(ux,uy,0.*xx)

title('campo de deformaciones')

}}

==Sólido antes y despues de desplazamiento==
El desplazamiento realizado por la placa viene determinado por una serie de movimientos transversales en direcciones distintas. En las imagenes adjuntas calculadas mediante Matlab, se aprecia que el desplazamiento es practicamente nulo.
[[Archivo:antesydespues.jpg|500px|derecha|Placa antes y después del desplazamiento]]
{{matlab|codigo=
h=0.1;
ro=1:h:2; % Valor de ro según el muestreo del intervalo
teta=0:h:2*pi+h; % Valor de teta según el muestreo del intervalo
[roro,tetateta]=meshgrid(ro,teta); %Matrices de la placa en coordenadas polares
figure(1)
xx=roro.*cos(tetateta); % parametrización
yy=roro.*sin(tetateta);
subplot(1,3,1) % figura solido momento inicial
mesh(xx,yy,0*xx) % mallado
axis([-4,4,-4,4]); % ejes
view(2); %visualización
desx=xx+((1-roro).^2).*(-yy); % movimiento parametrizado
desy=yy+((1-roro).^2).*(xx);
subplot(1,3,2) % figura solido en movimiento
mesh(desx,desy,yy.*0) % mallado
axis([-4,4,-4,4]);
view(2);
subplot(1,3,3); % figura solido en movimiento y campo vectorial
hold on
mesh(desx,desy,yy.*0)
axis([-4,4,-4,4]); % ejes
view(2); %
quiver(xx,yy,fx,fy) % campo vectorial del gradiente de temperaturas
axis([-4,4,-4,4]) % ejes
view(2)
hold off

}}

==Divergencia de u==
La divergencia mide el cambio de volumen local, en este caso en nuestra placa, debido al desplazamiento.
Con el grafico adjunto observamos que todos los puntos tienen el mismo valor, cero.
Como demostramos en el apartado anterior, nuestra placa esta sometida a desplazamientos transversales sin variacion de volumen, por lo que no hay cambio del mismo.
[[Archivo:divergence.jpg|350px|derecha|Divergencia en la placa]]
{{matlab|codigo=
h=(0.1);
ro=1:h:2;
teta=0:h:2*pi+h;
[roro,tetateta]=meshgrid(ro,teta);
figure(1)
xx=roro.*cos(tetateta);
yy=roro.*sin(tetateta)
f=(1./roro)-4+3.*roro; % El campo divergencia
surf(xx,yy,f) % Dibujmos la malla
axis([-2,2,-2,2])
colorbar
title('Divergencia del campo del desplazamiento')
view(2)
}}

==Rotacional de u==
A traves del programa de Matlab, que podemos observar a continuacion, hemos calculado el rotacional del campo desplazamiento u en los puntos del solido y los hemos representado, observandose en la grafica cuales son los puntos que sufren mayor rotacional.
[[Archivo:rotacionalc4.jpg|400px|derecha|Rotacional en la placa]]
{{matlab|codigo=
h=(0.1);
ro=1:h:2;
teta=0:h:2*pi+h;
[roro,tetateta]=meshgrid(ro,teta);
figure(1)
xx=roro.*cos(tetateta);
yy=roro.*sin(tetateta);
f=(2+4.*roro.^2-6.*roro);
surf(xx,yy,f)
axis([-2,2,-2,2])
colorbar
title('Rotacional del campo del desplazamiento')
view(2)

}}

== Estudio del tensor de tensiones ==
Definamos <math>\epsilon(\vec u)=(\nabla \vec u + \nabla \vec u^t)/2</math>, la parte simétrica del tensor gradiente de <math>\vec u</math>, que se denomina tensor de deformaciones. En un medio elástico lineal, isótropo y homogéneo los desplazamientos permiten escribir el tensor de tensiones <math>\sigma_{ij}</math> a través de la fórmula::
<math>
\sigma_{ij}=\lambda \nabla \cdot \vec u \delta_{ij} + 2\mu \epsilon_{ij},
</math>
Siendo <math>\nabla \cdot \vec u = \frac{1}{\sqrt{g}}\cdot \ \frac{\partial(\sqrt{g}u^{i})}{\partial x^{i}}=\frac{1}{\rho}\cdot (\frac{\partial \rho - 2\rho ^{2}+ \rho ^{3}}{\partial \rho }+\frac{\partial 0}{\partial \theta })=\frac{1}{\rho }-4+3\rho
</math>

y calculando el gradiente de <math>\vec u</math> como:: <math>\nabla\vec u=u^i_{.j}\vec g_i \otimes \vec g^j</math> donde

<math>u^i_{.j}= \vec g^i \nabla\vec u \vec g_j=\frac{\partial u^i}{\partial x^j}+\Gamma^i_{kj}u^k.= \left(\begin{array}{ccc}-2+2\rho& 0 \\0& \ 0\end{array}\right)+ (1-2\rho +\rho ^{2}).\left(\begin{array}{ccc}0& 0 \\0& \ 1\rho\end{array}\right)
</math>

Podemos calcular el tensor gradiente como 1-contravariante 1-covariante:
:<math>\nabla\vec u=u^i_{.j}=\left(\begin{array}{ccc}-2+2\rho& 0 \\0& \ 1/\rho -2 +\rho\end{array}\right)
</math>

Las coordenadas 1-covariantes 1-contravariante coinciden por lo que la transposicion seá directa y el tensor quedara como <math>(\epsilon^i_{.j})=(\nabla \vec u + \nabla \vec u^t)/2</math> , que corrresponde con la matriz mostrada anteriormente

Sustituyendo ahora los valores en la fórmula del tensor de tensiones,y tomando \(\lambda=\mu=1\), siendo \(\lambda\) y \(\mu\) los conocidos como coeficientes de Lamé, obtenemos:

:<math>\left(\begin{array}{ccc}1/\rho-8+7\rho& 0 \\0& \ 3/\rho -8 +5\rho\end{array}\right)</math>

Las tensiones normales de ambos vectores de la base nos quedaran como:

Las tensiones normales en la dirección \(\vec g_{\rho}\) serían:\[\vec g_{\rho} \cdot \sigma \cdot \vec g_\rho =1/\rho-8+7\rho\]

Las tensiones normales en la dirección \(\vec g_\theta/\rho\) serían:\[\vec g_\theta/\rho \cdot \sigma \cdot \vec g_\theta/\rho = 3/\rho -8 +5\rho\]

A continuación se muestra el código Matlab empleado para dibujar estas tensiones normales en ambas direcciones:
[[Archivo:tensionesnormales.jpg|500px|derecha|Tensiones normales en la placa]]
{{matlab|codigo=
h=0.1; % Intervalo de separación

u=[1:0.1:2]; % Intervalo de rho [1,2]

v=[0:h:2*pi+h]; % Intervalo de theta [0,2*pi]

[uu,vv]=meshgrid(u,v); % Matrices de rho and theta

figure(1)

subplot(1,2,1)

xx=uu.*cos(vv); % Parametrización

yy=uu.*sin(vv);

f=(1./uu)-8+7.*uu; % Campo escalar

surf(xx,yy,f) % Dibujar el mallado

axis([-2,2,-2,2]) % Región seleccionada

axis equal

colorbar

title('Tension normal en dirección g sub rho')

view(2)

subplot(1,2,2)

f=(3./uu)-8+5.*uu; % Campo escalar

surf(xx,yy,f) % Dibujar el mallado

axis([-2,2,-2,2]) % Región del dibujo

axis equal

colorbar

title('Tension normal en dirección g sub theta')

view(2)

}}

Vemos como la tensiones serán mayores a medida que nos alejamos del centro de masas de nuestra placa y simétricas respecto de este ya que dependerá de la distancia de cada uno de los puntos al origen. Además se aprecia que en los bordes interiores de las placas dicha tensión será nula.

===Tensiones tangenciales===
Al realizar el calculo de las tensiones tangenciales de nuestra placa para las direcciones <math>\vec g_\rho</math> y a <math>\vec g_\theta/\rho</math>::
<math>|\sigma \cdot \vec g_\rho-(\vec g_\rho \cdot \sigma \cdot \vec g_\rho) \vec g_\rho| = 0
</math><math>|\sigma \cdot \vec g_\theta/\rho-(\vec g_\theta/\rho \cdot \sigma \cdot \vec g_\theta/\rho) \vec g_\theta/\rho|= 0 </math>

Vemos que ambas tensiones son nulas, algebraicamente hablando esto es debido a que nuestro campo de desplazamientos depende únicamente del valor de la distancia al origen. Esto implica que al estar definido dicho tensor por la divergencia y el tensor unitario, será nula la tensión tangencial siempre que el campo no dependa del angulo, ya que sera igual en todas las direcciones.

==Tensión de Von Mises==
Se define por la fórmula: <math> \sigma _{VM}=\sqrt{\frac{(\sigma 1-\sigma 2)^{2}+(\sigma 2-\sigma 3)^{2}+(\sigma 3-\sigma 1)^{2}}{2}}</math>

==Masa total de la placa==
Si la densidad de la placa es <math> d(x,y,z)=xye^{-1/x^{2}} </math> y el área del anillo es <math> \pi.R^{2}-\pi.r^{2} </math> ,siendo R y r los radios exteriores e interiores respectivamente, la masa total de la placa será el producto del área y la densidad. Al tener la funcíon de la densidad, necesitamos integrarla para calcularla en todo el dominio.
Para una mejor resolución, efectuaremos el cambio a coordenadas polares:
:<math> x=\rho \cdot \cos \theta ; y=\rho \cdot \sin \theta </math>
Tendremos entonces que resolver la integral:
<math>\int_{0 }^{2\pi} \int_{1 }^{2}\rho ^{3}\cos \theta . \sin \theta .\log_{10}(\rho .\cos \theta +2)d\rho d\theta </math>
En matlab sería:
{{matlab|codigo=
h=1/100
u=1:h:2
v=0:h:2*pi
[uu,vv]=meshgrid(u,v)
f=(uu.^3).*sin(vv).*cos(vv).*log(uu.*cos(vv)+2)
a=(h^2)*f
masa=sum(sum(a))
}}
La masa de la placa es <math>m=5,0709.10^{-5}</math>

[[Categoría:Teoría de Campos]]
[[Categoría:TC14/15]]

Archivo:Divergence.jpg

2014-12-05T11:52:30Z

FernandoSancha:

Visualización de campos escalares y vectoriales en elasticidad. Grupo C4

2014-12-04T14:08:05Z

FernandoSancha: /* Campo de desplazamientos */

{{ TrabajoED |Visualización de campos escalares y vectoriales en elasticidad. Grupo C4 | [[:Categoría:Teoría de Campos|Teoría de Campos]]|[[:Categoría:TC14/15|2013-14]] | Alberto Rojas Rivero, Eduardo Bonet García, Fernando Sancha Domínguez, Carlos Fernandes Bermejo, Ricardo Mazón Cabrera, Emilio Valero Muñoz-Rojas }}

Consideramos una placa circular comprendida entre las circunferencias de radio 1 y radio 2, centradas en el origen.
En ella vamos a suponer que tenemos definidas dos cantidades físicas: la temperatura <math>T(\rho,\theta,t)</math> ,que depende de las dos variables espaciales <math>(x,y)</math> , y el tiempo <math>t</math>, y los desplazamientos <math>\vec u(\rho,\theta)</math>.De esta forma, si definimos <math>r_0(\rho,\theta)</math> el vector de posición de los puntos de la placa antes de la deformación , la posición de cada punto <math>(\rho,\theta)</math> de la placa después de la deformación viene dada por:
<math>
\vec r (\rho,\theta)= \vec r_{0}(\rho,\theta)+\vec u(\rho,\theta).
</math>
Vamos a suponer que la fuerza aplicada sobre la placa ha provocado una desplazamiento de manera que los desplazamientos en un instante <math>t_0</math> vienen dados por:
<math>\vec u(\rho,\theta)=(1-\rho)^2 \vec g_{\theta}.
</math>

== Mallado del sólido ==
Nuestro mallado consiste en la representación de los puntos interiores de una placa solida, en forma de corona circular. Tomando como paso de muestreo h=1/10, los comandos en Matlab de la visualizacion de nuestra placa es la siguiente:
[[Archivo:Sin título.jpg|350px|miniaturadeimagen|derecha|Mallado de la corona circular.]]
{{matlab|codigo=
h=0.1 % Paso de muestreo
u=1:h:2; % Intervalo [1,2]
v=0:h:2*pi+h; % Intervalo [0,2*pi]
[uu,vv]=meshgrid(u,v); % Matrices de las coordenadas u y v
figure(1)
xx=uu.*cos(vv); % parametrizacion
yy=uu.*sin(vv);
mesh(xx,yy,0*xx) % dibuja la malla
axis([-3,3,-3,3]) % selecciona la región
view(2) }}

== Temperatura del sólido ==
La temperatura es un campo escalar representado por <math>T(x,y)=e^{-y}</math>,vamos a tratar de ver como varía la temperatura en la placa que estamos estudiando . Nuestra función depende únicamente de la variable Y, esto quiere decir que la temperatura se mantendrá constante para todo valor de X. Observamos que cuanto mayor es el valor de la Y ,menor es la temperatura , luego los puntos de máxima temperatura son los más alejados del foco( el origen) y con menor valor de la variable Y, como podemos observar en la figura de la derecha.
[[Archivo:2.jpg|325px|derecha|Distribución de la temperatura en la placa]]
{{matlab|codigo=
h=0.1 % paso de muestreo
u=1:h:2; % Intervalo[1,2]
v=0:h:2*pi+h; % Intervalo [0,2*pi]
[uu,vv]=meshgrid(u,v); % matrices
figure(1)
xx=uu.*cos(vv); % parametrización
yy=uu.*sin(vv);
f=exp(-yy); % Campo escalar
surf(xx,yy,f) % visualización
axis([-3,3,-3,3]) % región de visualización
view(2)}}

==Cálculo y representación del <math>∇T </math>==

Para estudiar como varía la temperatura al movernos por el plano calculamos el gradiente de T.
<math>∇T= -e^{-y} j </math>
En la figura de la derecha podemos observar la representación del campo vectorial <math>∇T </math> y que este es perpendicular a las curvas de nivel del campo escalar T.
[[Archivo:Anillo doblado.jpg|400px|derecha| Curvas de nivel de T y gradiente de T]]
{{matlab|codigo=
g=-exp(-yy); % Gradiente
subplot(1,2,1); % Muestra varias imágenes
hold on
contour(xx,yy,f,20)
surf(xx,yy,f);
quiver3(xx,yy,f,0*xx,g,0*f);
hold off
subplot(1,2,2)
hold on
quiver(xx,yy,0*xx,Gradiente); colorbar;
contour(xx,yy,f,20)
hold off }}

==Campo de desplazamientos==
Supongamos la aplicación de una fuerza que provoca una vibración sobre la placa, de manera que, en un t aleatorio los desplazamientos vienen dados por:
<math>\vec u(\rho,\theta)=(1-\rho)^2 \vec g_{\theta}.
</math>
Para distinguir el campo de vectores en los puntos del mallado del sólido utilizaremos el código matlab siguiente:
[[Archivo:campodesplazam.jpg|500px|derecha|Distribución de la temperatura en la placCampo de desplazamientos de la placa]]
{{matlab|codigo=
u=1:h:2; % definir intervalo rho [1,2]

v=0:h:2*pi+h; % definir intervalo theta [0,2*pi]

[uu,vv]=meshgrid(u,v); % definir mallado

xx=uu.*cos(vv); % parametrizacion

yy=uu.*sin(vv);

ux=(xx.*(1-uu).^2);

uy=(yy.*(1-uu).^2);

subplot(1,2,1);

mesh(xx,yy,0.*xx)

view(2)

axis([-3,3,-3,3])

hold on

quiver(xx,yy,ux,uy)

title('campo con vectores')

subplot(1,2,2)

plot3(ux,uy,0.*xx)

title('campo de deformaciones')

}}

==Sólido antes y despues de desplazamiento==
El desplazamiento realizado por la placa viene determinado por una serie de movimientos transversales en direcciones distintas. En las imagenes adjuntas calculadas mediante Matlab, se aprecia que el desplazamiento es practicamente nulo.
[[Archivo:antesydespues.jpg|500px|derecha|Placa antes y después del desplazamiento]]
{{matlab|codigo=
h=0.1;
ro=1:h:2; % Valor de ro según el muestreo del intervalo
teta=0:h:2*pi+h; % Valor de teta según el muestreo del intervalo
[roro,tetateta]=meshgrid(ro,teta); %Matrices de la placa en coordenadas polares
figure(1)
xx=roro.*cos(tetateta); % parametrización
yy=roro.*sin(tetateta);
subplot(1,3,1) % figura solido momento inicial
mesh(xx,yy,0*xx) % mallado
axis([-4,4,-4,4]); % ejes
view(2); %visualización
desx=xx+((1-roro).^2).*(-yy); % movimiento parametrizado
desy=yy+((1-roro).^2).*(xx);
subplot(1,3,2) % figura solido en movimiento
mesh(desx,desy,yy.*0) % mallado
axis([-4,4,-4,4]);
view(2);
subplot(1,3,3); % figura solido en movimiento y campo vectorial
hold on
mesh(desx,desy,yy.*0)
axis([-4,4,-4,4]); % ejes
view(2); %
quiver(xx,yy,fx,fy) % campo vectorial del gradiente de temperaturas
axis([-4,4,-4,4]) % ejes
view(2)
hold off

}}

==Divergencia de u==
La divergencia mide el cambio de volumen local, en este caso en nuestra placa, debido al desplazamiento.
Con el grafico adjunto observamos que todos los puntos tienen el mismo valor, cero.
Como demostramos en el apartado anterior, nuestra placa esta sometida a desplazamientos transversales sin variacion de volumen, por lo que no hay cambio del mismo.
[[Archivo:divergenciac4.jpg|350px|derecha|Divergencia en la placa]]
{{matlab|codigo=
h=(0.1);
ro=1:h:2;
teta=0:h:2*pi+h;
[roro,tetateta]=meshgrid(ro,teta);
figure(1)
xx=roro.*cos(tetateta);
yy=roro.*sin(tetateta)
f=0*yy; % The scalar field
surf(xx,yy,f) % Draw the mesh
axis([-2,2,-2,2])
colorbar
title('Divergencia del campo del desplazamiento')
view(2)

}}

==Rotacional de u==
A traves del programa de Matlab, que podemos observar a continuacion, hemos calculado el rotacional del campo desplazamiento u en los puntos del solido y los hemos representado, observandose en la grafica cuales son los puntos que sufren mayor rotacional.
[[Archivo:rotacionalc4.jpg|400px|derecha|Rotacional en la placa]]
{{matlab|codigo=
h=(0.1);
ro=1:h:2;
teta=0:h:2*pi+h;
[roro,tetateta]=meshgrid(ro,teta);
figure(1)
xx=roro.*cos(tetateta);
yy=roro.*sin(tetateta);
f=(2+4.*roro.^2-6.*roro);
surf(xx,yy,f)
axis([-2,2,-2,2])
colorbar
title('Rotacional del campo del desplazamiento')
view(2)

}}

== Estudio del tensor de tensiones ==
Definamos <math>\epsilon(\vec u)=(\nabla \vec u + \nabla \vec u^t)/2</math>, la parte simétrica del tensor gradiente de <math>\vec u</math>, que se denomina tensor de deformaciones. En un medio elástico lineal, isótropo y homogéneo los desplazamientos permiten escribir el tensor de tensiones <math>\sigma_{ij}</math> a través de la fórmula::
<math>
\sigma_{ij}=\lambda \nabla \cdot \vec u \delta_{ij} + 2\mu \epsilon_{ij},
</math>
Siendo <math>\nabla \cdot \vec u = \frac{1}{\sqrt{g}}\cdot \ \frac{\partial(\sqrt{g}u^{i})}{\partial x^{i}}=\frac{1}{\rho}\cdot (\frac{\partial \rho - 2\rho ^{2}+ \rho ^{3}}{\partial \rho }+\frac{\partial 0}{\partial \theta })=\frac{1}{\rho }-4+3\rho
</math>

y calculando el gradiente de <math>\vec u</math> como:: <math>\nabla\vec u=u^i_{.j}\vec g_i \otimes \vec g^j</math> donde

<math>u^i_{.j}= \vec g^i \nabla\vec u \vec g_j=\frac{\partial u^i}{\partial x^j}+\Gamma^i_{kj}u^k.= \left(\begin{array}{ccc}-2+2\rho& 0 \\0& \ 0\end{array}\right)+ (1-2\rho +\rho ^{2}).\left(\begin{array}{ccc}0& 0 \\0& \ 1\rho\end{array}\right)
</math>

Podemos calcular el tensor gradiente como 1-contravariante 1-covariante:
:<math>\nabla\vec u=u^i_{.j}=\left(\begin{array}{ccc}-2+2\rho& 0 \\0& \ 1/\rho -2 +\rho\end{array}\right)
</math>

Las coordenadas 1-covariantes 1-contravariante coinciden por lo que la transposicion seá directa y el tensor quedara como <math>(\epsilon^i_{.j})=(\nabla \vec u + \nabla \vec u^t)/2</math> , que corrresponde con la matriz mostrada anteriormente

Sustituyendo ahora los valores en la fórmula del tensor de tensiones,y tomando \(\lambda=\mu=1\), siendo \(\lambda\) y \(\mu\) los conocidos como coeficientes de Lamé, obtenemos:

:<math>\left(\begin{array}{ccc}1/\rho-8+7\rho& 0 \\0& \ 3/\rho -8 +5\rho\end{array}\right)</math>

Las tensiones normales de ambos vectores de la base nos quedaran como:

Las tensiones normales en la dirección \(\vec g_{\rho}\) serían:\[\vec g_{\rho} \cdot \sigma \cdot \vec g_\rho =1/\rho-8+7\rho\]

Las tensiones normales en la dirección \(\vec g_\theta/\rho\) serían:\[\vec g_\theta/\rho \cdot \sigma \cdot \vec g_\theta/\rho = 3/\rho -8 +5\rho\]

A continuación se muestra el código Matlab empleado para dibujar estas tensiones normales en ambas direcciones:
[[Archivo:tensionesnormales.jpg|500px|derecha|Tensiones normales en la placa]]
{{matlab|codigo=
h=0.1; % Intervalo de separación

u=[1:0.1:2]; % Intervalo de rho [1,2]

v=[0:h:2*pi+h]; % Intervalo de theta [0,2*pi]

[uu,vv]=meshgrid(u,v); % Matrices de rho and theta

figure(1)

subplot(1,2,1)

xx=uu.*cos(vv); % Parametrización

yy=uu.*sin(vv);

f=(1./uu)-8+7.*uu; % Campo escalar

surf(xx,yy,f) % Dibujar el mallado

axis([-2,2,-2,2]) % Región seleccionada

axis equal

colorbar

title('Tension normal en dirección g sub rho')

view(2)

subplot(1,2,2)

f=(3./uu)-8+5.*uu; % Campo escalar

surf(xx,yy,f) % Dibujar el mallado

axis([-2,2,-2,2]) % Región del dibujo

axis equal

colorbar

title('Tension normal en dirección g sub theta')

view(2)

}}

Vemos como la tensiones serán mayores a medida que nos alejamos del centro de masas de nuestra placa y simétricas respecto de este ya que dependerá de la distancia de cada uno de los puntos al origen. Además se aprecia que en los bordes interiores de las placas dicha tensión será nula.

===Tensiones tangenciales===
Al realizar el calculo de las tensiones tangenciales de nuestra placa para las direcciones <math>\vec g_\rho</math> y a <math>\vec g_\theta/\rho</math>::
<math>|\sigma \cdot \vec g_\rho-(\vec g_\rho \cdot \sigma \cdot \vec g_\rho) \vec g_\rho| = 0
</math><math>|\sigma \cdot \vec g_\theta/\rho-(\vec g_\theta/\rho \cdot \sigma \cdot \vec g_\theta/\rho) \vec g_\theta/\rho|= 0 </math>

Vemos que ambas tensiones son nulas, algebraicamente hablando esto es debido a que nuestro campo de desplazamientos depende únicamente del valor de la distancia al origen. Esto implica que al estar definido dicho tensor por la divergencia y el tensor unitario, será nula la tensión tangencial siempre que el campo no dependa del angulo, ya que sera igual en todas las direcciones.

==Tensión de Von Mises==
Se define por la fórmula: <math> \sigma _{VM}=\sqrt{\frac{(\sigma 1-\sigma 2)^{2}+(\sigma 2-\sigma 3)^{2}+(\sigma 3-\sigma 1)^{2}}{2}}</math>

==Masa total de la placa==
Si la densidad de la placa es <math> d(x,y,z)=xye^{-1/x^{2}} </math> y el área del anillo es <math> \pi.R^{2}-\pi.r^{2} </math> ,siendo R y r los radios exteriores e interiores respectivamente, la masa total de la placa será el producto del área y la densidad. Al tener la funcíon de la densidad, necesitamos integrarla para calcularla en todo el dominio.
Para una mejor resolución, efectuaremos el cambio a coordenadas polares:
:<math> x=\rho \cdot \cos \theta ; y=\rho \cdot \sin \theta </math>
Tendremos entonces que resolver la integral:
<math>\int_{0 }^{2\pi} \int_{1 }^{2}\rho ^{3}\cos \theta . \sin \theta .\log_{10}(\rho .\cos \theta +2)d\rho d\theta </math>
En matlab sería:
{{matlab|codigo=
h=1/100
u=1:h:2
v=0:h:2*pi
[uu,vv]=meshgrid(u,v)
f=(uu.^3).*sin(vv).*cos(vv).*log(uu.*cos(vv)+2)
a=(h^2)*f
masa=sum(sum(a))
}}
La masa de la placa es <math>m=5,0709.10^{-5}</math>

[[Categoría:Teoría de Campos]]
[[Categoría:TC14/15]]