Studying the Kinetic Energy of a Flow

clc
clear
L=40;
M=50;
N=50;
LX=1;
LY=1;
LZ=1;
DX=LX/L;
DY=LY/M;
DZ=LZ/N;
g=9.81;
for i=1:L;
    CD(i)=2*10^(-3);
    R(i)=1.29;
    TV(i)=298+0.01*rand(1,1);
    TR(i)=0;
    FHSCS(i)=-0.02*randn(1,1);
    MW(i)=0.7+0.0001*rand(1,1);
    US(i)=(CD(i)*MW(i))^(0.5);
    UP(i)=0.01*randn(1,1);
    VP(i)=0.01*randn(1,1);
    WP(i)=0.01*randn(1,1);
    UA(i)=0.1*randn(1,1);
    VA(i)=0.1*randn(1,1);
    WA(i)=0.1*randn(1,1);
    TKE(i)=0.5*((UP(i))^(2)+(VP(i))^(2)+(WP(i))^(2));
    LEPS(i)=50;
    EPS(i)=((TKE(i))^(1.5))/LEPS(i);
    B(i)=(g/TV(i))*FHSCS(i);
end
for i=1:L-1;
    DTKE(i)=TKE(i+1)-TKE(i);
    A(i)=-(UA(i)*DTKE(i))/DX-(VA(i)*DTKE(i))/DY-(WA(i)*DTKE(i))/DZ;
    S(i)=(MW(i+1)-MW(i))*((US(i))^(2))/DZ;
    RF(i)=-B(i)/S(i);
    DTKE_DT(i)=A(i)+S(i)+B(i)+TR(i)-EPS(i);
end
NSAMPLES=(1:1:L-1);

scrsz = get(0,'ScreenSize');
figure('Position',[10 10 1500 800])
subplot 231
plot(UP,'r')
Xlabel('Number of Samples')
ylabel('UP')
grid on
subplot 232
plot(VP)
Xlabel('Number of Samples')
ylabel('VP')
grid on
subplot 233
plot(WP)
Xlabel('Number of Samples')
ylabel('WP')
grid on
subplot 234
plot(UA)
Xlabel('Number of Samples')
ylabel('UA')
grid on
subplot 235
plot(VA)
Xlabel('Number of Samples')
ylabel('VA')
grid on
subplot 236
plot(WA)
xlabel('Number of Samples')
ylabel('WA')
grid on
pause(10)
close all

scrsz = get(0,'ScreenSize');
figure('Position',[10 10 1500 800])
subplot 231
plot(RF,'r')
Xlabel('Number of Samples')
ylabel('RF')
grid on
subplot 232
plot(S)
Xlabel('Number of Samples')
ylabel('S')
grid on
subplot 233
plot(A)
Xlabel('Number of Samples')
ylabel('A')
grid on
subplot 234
plot(DTKE)
Xlabel('Number of Samples')
ylabel('DTKE')
grid on
subplot 235
plot(DTKE_DT)
Xlabel('Number of Samples')
ylabel('DTKE-DT')
grid on
subplot 236
plot(B(1:39),S(1:39),'*')
xlabel('B')
ylabel('S')
grid on
pause(10)
close all

for i=1:L-1;
    RF(i)=-B(i)/abs(S(i));
    theta(i)=acot(RF(i));
end

for i=1:L-1;
    x(i)=cos(theta(i));
    y(i)=abs(sin(theta(i)));
    xx(1)=0;
    yy(1)=0;
    v1=[ x(i) xx(1) ]
    v2=[ y(i) yy(1) ]
    f=line(v1,v2);
    hold on
end
dd=max(x)
ddd=max(y)
set(gca,'XLim',[-dd dd],'YLim',[-(ddd+0.9*ddd) (ddd+0.9*ddd)])
xlabel('B')
ylabel('s')
hold on
plot(x,y,'*')
axis equal
grid on

Unless otherwise noted, all content on this site is @Copyright by Ahmed Al Makky 2012-2013 - http://cfd2012.com