[archeoptyrx]

Hello everyone ,

I have wrote a code for solving compressible navier stokes code over a flat plate (same as given in the Anderson CFD Book) . when i tried simulating it with free stream conditions at higher altitude for example

Plate length of 1e-5 encounters slip regime above 10km altitude , which means that slip boundary conditions should be incorporated to get proper results , and the code blows up for free stream conditions above 20km (because transition regime comes into play )

Am i correct guys ? Navier stokes equation works only for continuum regime right ?

So i am using MACCORMACK Method with the time step given in book . Its because of the time step , which brings up complex numbers and its causes the program to stop .

What should i do now ?

Any suggestion will be of great help

Thanks in advance

[flotus1]

While the NS equations are only valid in the continuum regime, a NS solver should still be able to produce a (wrong) solution in the rarefied regime.
If your NS code diverges for higher Knudsen numbers, there must be something wrong with the code.

[archeoptyrx]

Alex,

I checked the code and i find that i implemented the time step (stability criteria) wrong . This produces the complex number first which makes the solver to blow up for higher Knudsen numbers

[archeoptyrx]

Take a look at my time step , I implemented it wrong .

function delt=Tstep(u,v,T,rho,mu)
% K         -Fudge factor 0.5<K<0.8
% delx      - xinterval
% dely      - y interval
% u         - horizontal velocity at (i,j)
% v         - Vertical velocity at (i,j)
% a         - Speed of sound at (i,j)
% mu        - Dynamic viscosity at (i,j)
% Pr        - Prandtl Number
% rho       - Density at (i,j)
% gamma     - iscentropic constant
global K Pr gamma imax jmax delx dely;
% global u v T rho mu;

a=Sound(T);
vda=zeros(imax,jmax);
deltcfl=zeros(imax,jmax);
for i=1:imax
    for j=1:jmax
%                mu(i,j)=Dynvis(T(i,j),2);
        vda(i,j)=((4*mu(i,j)*gamma*mu(i,j))/(3*rho(i,j)*Pr));
    end
end
vdas=max(max(vda(:)));
for i=1:imax
    for j=1:jmax
        A=abs(u(i,j))/delx;
        B=abs(v(i,j))/dely;
        C=a(i,j)*sqrt((1/(delx^2))+(1/(dely^2)));
        D=2*vdas*((1/(delx^2))+(1/(dely^2)));
        deltcfl(i,j)=(A+B+C+D)^(-1);
    end
end

delt=real(min(min(K*deltcfl(:))));
end

[obdwinston]

Hey there, just wanted to ask if you have found a solution to your problem (I know it had been about 5 years since you posted this thread).

I am currently working on the same problem (Chapter 10 of J.D. Anderson's CFD book) and encountered the same issue with complex numbers showing up in the flow variables. A bit of digging led me to the time step being too large. I used the time step formula on page 457 and got a time step of 6.1850e-12 for the first iteration (with the initial conditions). When I reduce the time step further to 1e-15, the complex numbers after the first iteration disappeared. I am still working on finishing up the rest of the code (I have yet to run a complete simulation), but just decided I should post this in case you still read the forum.

Do you have any idea if the time step formula (page 457) is accurate? Also, may I have the approximate time step you used (just the order of magnitude is fine)?

Please see my MATLAB code below for the part that computes the time step.

Cheers!

K = 0.6; % Courant number
A = abs(u)/dx;
B = abs(v)/dy;
C = sqrt(1/dx/dx + 1/dy/dy)*sqrt(gamma*R*T);
D = ones(i_max,j_max)*2*(1/dx/dx + 1/dy/dy)*max(max((4/3)*(gamma/Pr)*mu.*mu./rho));
dt = min(min(K./(A+B+C+D)));

[obdwinston]

I found a mistake in my code. Missed out a square root to compute del (boundary layer thickness)