- #1

nitish

- 1

- 0

I have one serious problem while solving rayleigh equation using blasius profile, in which so as to remove the singularity the intergration contour is defined in a complex plane.

4th order runge kutta is used but if the step size( h) is is in complex, it is giving some error. would anyone please help me solving it?

here is the MATLAB code. thank you...reply as soon as possible.

nstep=100;

pi=3.14;

error=1e-4;

ycr=1.35;

yci=1.35;

ymax=8;

nstep2=5*nstep-1;

ncirc=nstep;

ncirc2=3*ncirc+1;

for j=1:ncirc2

theta=pi*(j-1)/(ncirc2-1);

yr=ycr*(1-cos(theta));

yi=-yci*sin(theta);

Y(j)= (yr+i*yi);

end

for j=ncirc2:nstep2

yr=ycr*2+ymax*(j-ncirc2)/(nstep2-ncirc2);

Y(j)=(yr+i*0);

end

%y=zeros(1,61);

%z=zeros(1,61);

%p=zeros(1,61);

%u=zeros(1,61);

F(1)=0+0*i;

U(1)=0+0*i;

UP(1)=1.660287+0*i;

UPP(1)=0+0*i;

[y,z,p,u]=blasius1(nstep2,Y,F,U,UP,UPP);

pnew=p(1);

znew=z(nstep2);

UP(1)=1.05*UP(1);

[y,z,p,u]=blasius1(nstep2,Y,F,U,UP,UPP);

pold=p(1);

zold=z(nstep2);

residue=1e30;

while(residue>error)

U(1),UPP(1),F(1),Y,nstep2;

UP(1)=pnew-(((znew)-1)*((pold)-(pnew)))/((zold)-(znew));

[y,z,p,u]=blasius1(nstep2,Y,F,U,UP,UPP);

zold=znew;

pold=pnew;

znew=z(nstep2);

pnew=p(1);

residue=abs((z(nstep2))-1.0);

end

%plot(real(Y),real(z));

hold on

plot(real(Y),abs(p));

hold on

%plot(real(Y),real(y));

hold on

%plot(real(Y),real(u));

================================================== ===============

blasius.m(separate file):

function[y,z,p,u]=blasius1(num,Y,F,U,UP,UPP)

y(1)=F(1);

z(1)=U(1);

p(1)=UP(1);

u(1)=UPP(1);

for j=1num-1)

h=((Y(j+1))-(Y(j)))*0.5;

h=real(h);

A11=(h)*(z(j));

A12=(h)*(p(j));

A13=(-1*(h)*((y(j))*(p(j))));

A14=(-(1.4806)*(h)*(((y(j))*(u(j))+(z(j))*(p(j)))));

A21=(h)*((z(j))+(A12)*0.5);

A22=(h)*((p(j))+(A13)*0.5);

A23=(-1*(h)*(((y(j))+(A11)*0.5)*((p(j))+(A13)*0.5)));

A24=(-(1.4806)*(h)*((((y(j))+(A11)*0.5)*((u(j))+(A14)*0. 5)+((z(j))+(A12)*0.5)*((p(j))+(A13)*0.5))));

A31=(h)*((z(j))+(A22)*0.5);

A32=(h)*((p(j))+(A23)*0.5);

A33=(-1*(h)*(((y(j))+(A21)*0.5)*((p(j))+(A23)*0.5)));

A34=(-(1.4806)*(h)*((((y(j))+(A21)*0.5)*((u(j))+(A24)*0. 5)+((z(j))+(A22)*0.5)*((p(j))+(A23)*0.5))));

A41=(h)*((z(j))+(A32));

A42=(h)*((p(j))+(A33));

A43=(-1*(h)*(((y(j))+(A31))*((p(j))+(A33))));

A44=(-(1.4806)*(h)*((((y(j))+(A31))*((u(j))+(A34))+((z(j ))+(A32))*((p(j))+(A33)))));

y(j+1)=((y(j))+(((A11)+2*(A21)+2*(A31)+(A41))/6));

z(j+1)=((z(j))+(((A12)+2*(A22)+2*(A32)+(A42))/6));

p(j+1)=((p(j))+(((A13)+2*(A23)+2*(A33)+(A43))/6));

u(j+1)=((u(j))+(((A14)+2*(A24)+2*(A34)+(A44))/6));

end