Help calculating the thrust coefficient

Based on the code provided

My thrust coefficient is half of what is supossed to be. Can anyone point out what my mistakes are on my equation. The code starts from %%verification of conservation of momentum. I know based on my TSR my thrust coefficient is supposed to be somewhere around .7

Also my professor helped me calculate the power coefficient but i do not understand what this lines of code are doing?

What is tke1 and tke2 doing?

i know i have to do the same thing for the i,j and i,k but how would i do it?

could you do it for i,j component and comment how what is that it is doing it so that i can understand and do it fo i,k.

tke1=(u(istart,j,k)^2.0+v(istart,j,k)^2.0+w(istart,j,k)^2.0)/2;
tke2=(u(iend,j,k)^2.0+v(iend,j,k)^2.0+w(iend,j,k)^2.0)/2;

power_1=power_1+((-u(istart,j,k)))*((tke1)+pr(istart,j,k))*da;
power_2=power_2+(u(iend,j,k))*((tke2)+pr(iend,j,k))*da;

D = 126; %diamter
R = D/2;
Re = 7.3*10^7; %reynolds number
Rho = 1.225; %density
Mu = 1.785*10^-5 %viscosity of air
velocity = Re*Mu/Rho/D %velocity profile 
TSR = 7.5 %tip speed ration
Omega = TSR*velocity/R %angular velocity 
Power_in = 1/2*Rho*pi*R^2*velocity^3
figure 
utmpx(:,:)=u(:,:,n3/2);
contourf(xx,yy,utmpx',100, 'linestyle','none')
daspect([1,1,1])

figure 
[~,ind]=min(abs(xx-3));
utmpz(:,:)=u(ind,:,:);
contourf(zz,yy,utmpz,80, 'linestyle','none')
daspect([1,1,1])

%flux  going through the faces
f1=0;
f2=0;
for j=1: n2-1
    for k=1:n3-1
        da=(yy(j+1)-yy(j))*(zz(k+1)-zz(k));
        f1=f1+u(1,j,k)*da;
        f2=f2+u(n1,j,k)*da;
    end 
end
f3=0;
f4=0;
for i=1:n1-1
    for j=1:n2-1
        da=(yy(j+1)-yy(j))*(xx(i+1)-xx(i));
        f3=f3+w(i,j,1)*da;
        f4=f4+w(i,j,n3)*da;
    end
end
f5=0;
f6=0;
for i=1:n1-1
    for k=1:n3-1
        da=(xx(i+1)-xx(i))*(zz(k+1)-zz(k));
        f5=f5+v(i,n2,k)*da;
        f6=f6+v(i,1,k)*da;
    end
end
total_flux=-f1+f2+f3-f4+f5-f6

%%verification of conservation of momentum
f1=0;
f2=0;
p1=0;
p2=0;

f1=0;
f2=0;
for j=1: n2-1
    for k=1:n3-1
        da=(yy(j+1)-yy(j))*(zz(k+1)-zz(k));
        f1=f1+u(1,j,k)^2*da;
        f2=f2+u(n1,j,k)^2*da;
        p1=p1+pr(1,j,k)*da;
        p2=p2+pr(n1,j,k)*da;
    end 
end
f3=0;
f4=0;
for i=1:n1-1
    for j=1:n2-1
        da=(yy(j+1)-yy(j))*(xx(i+1)-xx(i));
        f3=f3+u(i,j,1)*w(i,j,1)*da;
        f4=f4+u(i,j,n3)*w(i,j,n3)*da;
    end
end
f5=0;
f6=0;
for i=1:n1-1
    for k=1:n3-1
        da=(xx(i+1)-xx(i))*(zz(k+1)-zz(k));
        f5=f5+u(i,n2,k)*v(i,n2,k)*da;
        f6=f6+u(i,1,k)*v(i,1,k)*da;
    end
end
trust_coefficient =-f1+f2+f3-f4+f5-f6+p1-p2

%Thrust Force 
TF = .5*trust_coefficient*1.225*velocity^2*pi*63^2

% Power extracted by the wind turbine
power_1=0;
power_2=0;
istart=75
iend=250
for j=1:n2-1
    for k=1:n3-1
        da=(zz(k+1)-zz(k))*(yy(j+1)-yy(j));
        tke1=(u(istart,j,k)^2.0+v(istart,j,k)^2.0+w(istart,j,k)^2.0)/2;
        tke2=(u(iend,j,k)^2.0+v(iend,j,k)^2.0+w(iend,j,k)^2.0)/2;

        power_1=power_1+((-u(istart,j,k)))*((tke1)+pr(istart,j,k))*da;
        power_2=power_2+(u(iend,j,k))*((tke2)+pr(iend,j,k))*da;
   end
end
power_3=0;
power_4=0;
for i=1:n1-1
    for j=1:n2-1
        da=(yy(j+1)-yy(j))*(xx(i+1)-xx(i));
        power_3=power_3+((-w(i,j,1)))*(((u(i,j,1)^2.0)/2)+pr(i,j,1))*da;
        power_4=power_4+(w(i,j,n3))*(((u(i,j,n3)^2.0)/2)+pr(i,j,n3))*da;
    end
end
power_5=0;
power_6=0;
for i=1:n1-1
    for k=1:n3-1
        da=(xx(i+1)-xx(i))*(zz(k+1)-zz(k));
        power_5=power_5+((v(i,100,k)))*(((u(i,100,k)^2.0)/2)+pr(i,100,k))*da;
        power_6=power_6+(-v(i,1,k))*(((u(i,1,k)^2.0)/2)+pr(i,1,k))*da;
    end
end
WT_power_1=4/3.14*(power_1+power_2)*2
WT_power_2=power_3+power_4
WT_power_3=power_5+power_6
c_p = -(WT_power_1+WT_power_2+WT_power_3)

%For visualisationn purpose Z and Y are transposed 
figure 
hold on
[X,Y,Z]=meshgrid(xx,zz,yy);
U=permute(u,[3,1,2]);
W=permute(v,[3,1,2]);
V=permute(w,[3,1,2]);

%Particles lines from 
[sx,sy,sz]=meshgrid(0,linspace(1.25,1.75,5),linspace(.4,.8,5));
h=streamline(X,Y,Z,U,V,W,sx,sy,sz);
set(h,'color','green');
daspect([1,1,1])
axis([min(xx),max(xx),min(zz),max(zz),min(yy),max(yy)])
z0=1.5;y0=.7031;
theta=0:.1:2*pi;
zc=.5*cos(theta)+z0;
yc=.5*sin(theta)+y0;
fill3(ones(1,length(zc))*3,zc,yc,[.6,.6,.6])
grid on
view(3)
xlabel('x')
ylabel('y')
zlabel('z')
        

% for visualization Z and Y are transposed 
figure 
hold on
[X,Y,Z]=meshgrid(xx,zz,yy);
U=permute(u,[3,1,2]);
W=permute(v,[3,1,2]);
V=permute(w,[3,1,2]);
%particles lines released from rotor plane 
[sx,sy,sz]=meshgrid(3,linspace(1.25,1.75,5),linspace(.4,.8,5));
h=streamline(X,Y,Z,U,V,W,sx,sy,sz);
set(h,'color','red');
daspect([1,1,1])
axis([min(xx),max(xx),min(zz),max(zz),min(yy),max(yy)])
z0=1.5;y0=.7031;
theta=0:.1:2*pi;
zc=.5*cos(theta)+z0;
yc=.5*sin(theta)+y0;
fill3(ones(1,length(zc))*3,zc,yc,[.6,.6,.6])
grid on
view(3)
xlabel('x')
ylabel('y')
zlabel('z')
%Pressure plot
figure 
pr_vsec(:,:)=pr(:,:,n3/2);
contourf(xx,yy,pr_vsec',180,'linestyle',' none')
daspect([1,1,1])
xlabel('x')
ylabel('y')

power = 1/2*Rho*velocity^3*pi*R^2*c_p

Answer 1

[TAMILAN]

D = 126; %diamter

R = D/2;

Re = 7.3*10^7; %reynolds number

Rho = 1.225; %density

Mu = 1.785*10^-5 %viscosity of air

velocity = Re*Mu/Rho/D %velocity profile

TSR = 7.5 %tip speed ration

Omega = TSR*velocity/R %angular velocity

Power_in = 1/2*Rho*pi*R^2*velocity^3

figure

utmpx(:,:)=u(:,:,n3/2);

contourf(xx,yy,utmpx',100, 'linestyle','none')

daspect([1,1,1])

figure

[~,ind]=min(abs(xx-3));

utmpz(:,:)=u(ind,:,:);

contourf(zz,yy,utmpz,80, 'linestyle','none')

daspect([1,1,1])

%flux going through the faces

f1=0;

f2=0;

for j=1: n2-1

for k=1:n3-1

da=(yy(j+1)-yy(j))*(zz(k+1)-zz(k));

f1=f1+u(1,j,k)*da;

f2=f2+u(n1,j,k)*da;

end

end

f3=0;

f4=0;

for i=1:n1-1

for j=1:n2-1

da=(yy(j+1)-yy(j))*(xx(i+1)-xx(i));

f3=f3+w(i,j,1)*da;

f4=f4+w(i,j,n3)*da;

end

end

f5=0;

f6=0;

for i=1:n1-1

for k=1:n3-1

da=(xx(i+1)-xx(i))*(zz(k+1)-zz(k));

f5=f5+v(i,n2,k)*da;

f6=f6+v(i,1,k)*da;

end

end

total_flux=-f1+f2+f3-f4+f5-f6

%%verification of conservation of momentum

f1=0;

f2=0;

p1=0;

p2=0;

f1=0;

f2=0;

for j=1: n2-1

for k=1:n3-1

da=(yy(j+1)-yy(j))*(zz(k+1)-zz(k));

f1=f1+u(1,j,k)^2*da;

f2=f2+u(n1,j,k)^2*da;

p1=p1+pr(1,j,k)*da;

p2=p2+pr(n1,j,k)*da;

end

end

f3=0;

f4=0;

for i=1:n1-1

for j=1:n2-1

da=(yy(j+1)-yy(j))*(xx(i+1)-xx(i));

f3=f3+u(i,j,1)*w(i,j,1)*da;

f4=f4+u(i,j,n3)*w(i,j,n3)*da;

end

end

f5=0;

f6=0;

for i=1:n1-1

for k=1:n3-1

da=(xx(i+1)-xx(i))*(zz(k+1)-zz(k));

f5=f5+u(i,n2,k)*v(i,n2,k)*da;

f6=f6+u(i,1,k)*v(i,1,k)*da;

end

end

trust_coefficient =-f1+f2+f3-f4+f5-f6+p1-p2

%Thrust Force

TF = .5*trust_coefficient*1.225*velocity^2*pi*63^2

% Power extracted by the wind turbine

power_1=0;

power_2=0;

istart=75

iend=250

for j=1:n2-1

for k=1:n3-1

da=(zz(k+1)-zz(k))*(yy(j+1)-yy(j));

tke1=(u(istart,j,k)^2.0+v(istart,j,k)^2.0+w(istart,j,k)^2.0)/2;

tke2=(u(iend,j,k)^2.0+v(iend,j,k)^2.0+w(iend,j,k)^2.0)/2;

power_1=power_1+((-u(istart,j,k)))*((tke1)+pr(istart,j,k))*da;

power_2=power_2+(u(iend,j,k))*((tke2)+pr(iend,j,k))*da;

end

end

power_3=0;

power_4=0;

for i=1:n1-1

for j=1:n2-1

da=(yy(j+1)-yy(j))*(xx(i+1)-xx(i));

power_3=power_3+((-w(i,j,1)))*(((u(i,j,1)^2.0)/2)+pr(i,j,1))*da;

power_4=power_4+(w(i,j,n3))*(((u(i,j,n3)^2.0)/2)+pr(i,j,n3))*da;

end

end

power_5=0;

power_6=0;

for i=1:n1-1

for k=1:n3-1

da=(xx(i+1)-xx(i))*(zz(k+1)-zz(k));

power_5=power_5+((v(i,100,k)))*(((u(i,100,k)^2.0)/2)+pr(i,100,k))*da;

power_6=power_6+(-v(i,1,k))*(((u(i,1,k)^2.0)/2)+pr(i,1,k))*da;

end

end

WT_power_1=4/3.14*(power_1+power_2)*2

WT_power_2=power_3+power_4

WT_power_3=power_5+power_6

c_p = -(WT_power_1+WT_power_2+WT_power_3)

%For visualisationn purpose Z and Y are transposed

figure

hold on

[X,Y,Z]=meshgrid(xx,zz,yy);

U=permute(u,[3,1,2]);

W=permute(v,[3,1,2]);

V=permute(w,[3,1,2]);

%Particles lines from

[sx,sy,sz]=meshgrid(0,linspace(1.25,1.75,5),linspace(.4,.8,5));

h=streamline(X,Y,Z,U,V,W,sx,sy,sz);

set(h,'color','green');

daspect([1,1,1])

axis([min(xx),max(xx),min(zz),max(zz),min(yy),max(yy)])

z0=1.5;y0=.7031;

theta=0:.1:2*pi;

zc=.5*cos(theta)+z0;

yc=.5*sin(theta)+y0;

fill3(ones(1,length(zc))*3,zc,yc,[.6,.6,.6])

grid on

view(3)

xlabel('x')

ylabel('y')

zlabel('z')

% for visualization Z and Y are transposed

figure

hold on

[X,Y,Z]=meshgrid(xx,zz,yy);

U=permute(u,[3,1,2]);

W=permute(v,[3,1,2]);

V=permute(w,[3,1,2]);

%particles lines released from rotor plane

[sx,sy,sz]=meshgrid(3,linspace(1.25,1.75,5),linspace(.4,.8,5));

h=streamline(X,Y,Z,U,V,W,sx,sy,sz);

set(h,'color','red');

daspect([1,1,1])

axis([min(xx),max(xx),min(zz),max(zz),min(yy),max(yy)])

z0=1.5;y0=.7031;

theta=0:.1:2*pi;

zc=.5*cos(theta)+z0;

yc=.5*sin(theta)+y0;

fill3(ones(1,length(zc))*3,zc,yc,[.6,.6,.6])

grid on

view(3)

xlabel('x')

ylabel('y')

zlabel('z')

%Pressure plot

figure

pr_vsec(:,:)=pr(:,:,n3/2);

contourf(xx,yy,pr_vsec',180,'linestyle',' none')

daspect([1,1,1])

xlabel('x')

ylabel('y')

power = 1/2*Rho*velocity^3*pi*R^2*c_p