Help calculating the thrust coefficient
12 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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
0 commentaires
Réponses (1)
TAMILAN
le 3 Juil 2024
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
0 commentaires
Voir également
Catégories
En savoir plus sur Fluid Mechanics dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!