CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Main CFD Forum

Navier stokes and maxwell slip boundary

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree1Likes

Reply
 
LinkBack Thread Tools Display Modes
Old   July 22, 2013, 15:20
Question Navier stokes and maxwell slip boundary
  #1
Member
 
sud
Join Date: May 2013
Posts: 65
Rep Power: 4
archeoptyrx is on a distinguished road
Hello ppl ,

I wrote a MATLAB code for the compressible navier stokes equation as per the Anderson CFD book . now I am trying to implement the Maxwell slip and temperature jump boundary conditions . its in first order .

Could anyone help me in implementing this to the code ?
I am no familiar with this Neumann type of boundary conditions .

I have a grid of 70*70 . with flat plate at y=0 and ranging from 0<=x<=L

Previously I just used ,

For i=2:imax
U(i,1)=0;
T(i,1)=300;
End

How should I modify it ???
Looking for someone to help me understand in implementinv this derivative boundary conditions.


Thanks
archeoptyrx is offline   Reply With Quote

Old   July 22, 2013, 16:20
Default
  #2
Senior Member
 
flotus1's Avatar
 
Alex
Join Date: Jun 2012
Location: Germany
Posts: 1,098
Rep Power: 19
flotus1 will become famous soon enoughflotus1 will become famous soon enough
You have several options for the implementation of a first order Maxwell boundary condition.

The easiest (for finite volumes) would be the way Fluent does it. They obtain the wall-normal gradients for velocity and temperature from finite differences between the values at the wall and a the first cell center.
Then they apply the value calculated from the Maxwell boundary condition to the boundary.
So you need a temporal storage to store the boundary value from the last iteration to use it in the finite differences for the next iteration.

Keep in mind that in doing so, the spatial accuracy of the whole computation is reduced to first order. To prevent this, you could for example evaluate the wall-normal gradient at the wall with a higher order finite difference scheme that uses additional points further away from the wall.

Some kind of under-relaxation might be necessary to keep this kind of BC stable, especially with higher order finite difference schemes for the gradients.
flotus1 is offline   Reply With Quote

Old   July 22, 2013, 16:45
Question
  #3
Member
 
sud
Join Date: May 2013
Posts: 65
Rep Power: 4
archeoptyrx is on a distinguished road
OK . my Matlab code is explicit mccormack method . now . I set up the free stream conditions .

Next for the first iteration , I m struck in initializing the domain as in t pic below .

At wall ,

What should I do for velocity and temperature? For t first iteration . should I include ,
Us=0
Two=300;

And for t successive iteration , I attached a pic . please take a look and help me out if you can .

Thanks

Quote:
Originally Posted by flotus1 View Post
You have several options for the implementation of a first order Maxwell boundary condition.

The easiest (for finite volumes) would be the way Fluent does it. They obtain the wall-normal gradients for velocity and temperature from finite differences between the values at the wall and a the first cell center.
Then they apply the value calculated from the Maxwell boundary condition to the boundary.
So you need a temporal storage to store the boundary value from the last iteration to use it in the finite differences for the next iteration.

Keep in mind that in doing so, the spatial accuracy of the whole computation is reduced to first order. To prevent this, you could for example evaluate the wall-normal gradient at the wall with a higher order finite difference scheme that uses additional points further away from the wall.

Some kind of under-relaxation might be necessary to keep this kind of BC stable, especially with higher order finite difference schemes for the gradients.
Attached Images
File Type: jpg 1374525702623.jpg (92.7 KB, 17 views)
File Type: jpg 1374525865000.jpg (95.2 KB, 13 views)
archeoptyrx is offline   Reply With Quote

Old   July 23, 2013, 05:21
Default
  #4
Senior Member
 
flotus1's Avatar
 
Alex
Join Date: Jun 2012
Location: Germany
Posts: 1,098
Rep Power: 19
flotus1 will become famous soon enoughflotus1 will become famous soon enough
From what I recall, the non-isothermal slip boundary conditions for flat walls read:

u_s^t-u_w^t = \frac{2-\sigma_v}{\sigma_v} \lambda_f \frac{\partial u^t}{\partial n} +\text{neglected term from tangential temperature gradient}

where superscript t denotes the wall-tangential direction and \partial n denotes the derivative in wall-normal direction

T_s-T_w = \frac{2-\sigma_t}{\sigma_t} \frac{2\gamma}{\gamma+1} \frac{\lambda_f}{Pr} \frac{\partial T}{\partial n}

The tangential velocity of the wall u_w^t appears to be 0 in your setup while the temperature of the wall T_w seems to be 300K.
Given the macroscopic gradients, this provides you with explicit expressions for the boundary values to be set in the computation.

Suitable initial conditions for slip velocity and temperature could be the macroscopic values of these quantities at the wall (0 m/s and 300K).

Maybe this document provides some guidance for the implementation of the code:
http://fluid.ippt.gov.pl/seminar/text/steffen110407.pdf
archeoptyrx likes this.
flotus1 is offline   Reply With Quote

Old   July 23, 2013, 15:27
Default
  #5
Member
 
sud
Join Date: May 2013
Posts: 65
Rep Power: 4
archeoptyrx is on a distinguished road
Thanks a lot Mr.Alex .

I will try running the code now . Let me post the results in some couple of hours
archeoptyrx is offline   Reply With Quote

Old   July 29, 2013, 15:22
Unhappy
  #6
Member
 
sud
Join Date: May 2013
Posts: 65
Rep Power: 4
archeoptyrx is on a distinguished road
Quote:
Originally Posted by flotus1 View Post
From what I recall, the non-isothermal slip boundary conditions for flat walls read:

u_s^t-u_w^t = \frac{2-\sigma_v}{\sigma_v} \lambda_f \frac{\partial u^t}{\partial n} +\text{neglected term from tangential temperature gradient}

where superscript t denotes the wall-tangential direction and \partial n denotes the derivative in wall-normal direction

T_s-T_w = \frac{2-\sigma_t}{\sigma_t} \frac{2\gamma}{\gamma+1} \frac{\lambda_f}{Pr} \frac{\partial T}{\partial n}

The tangential velocity of the wall u_w^t appears to be 0 in your setup while the temperature of the wall T_w seems to be 300K.
Given the macroscopic gradients, this provides you with explicit expressions for the boundary values to be set in the computation.

Suitable initial conditions for slip velocity and temperature could be the macroscopic values of these quantities at the wall (0 m/s and 300K).

Maybe this document provides some guidance for the implementation of the code:
http://fluid.ippt.gov.pl/seminar/text/steffen110407.pdf
HI Alex ,

I tried running in the way u said . But after 40000 iterations , it didn end . When i plottted the contour of Pressure , it showed some very weird plot. what could be the problem ?

My model is same as Anderson Flat plate of same length .

so LHORI = 1e-5;

number density =8..266e24;

should produce a knudsen number of KN=0.0156 which lies in the slip regime .

And Maxwell first order boundary condition should work in this regime . But for me it didnt even produce any reasonable results .

Any suggestions ??
archeoptyrx is offline   Reply With Quote

Old   July 29, 2013, 15:33
Default
  #7
Member
 
sud
Join Date: May 2013
Posts: 65
Rep Power: 4
archeoptyrx is on a distinguished road
Mr.Alex ,

If possible can you take a look at these files.
Attached Images
File Type: jpg Capture.jpg (30.2 KB, 14 views)
archeoptyrx is offline   Reply With Quote

Old   July 29, 2013, 16:58
Default
  #8
Senior Member
 
flotus1's Avatar
 
Alex
Join Date: Jun 2012
Location: Germany
Posts: 1,098
Rep Power: 19
flotus1 will become famous soon enoughflotus1 will become famous soon enough
I assume your code produces reasonable results without the Maxwell boundary condition?

I will have a closer look at the code tomorrow.
For now, you could try to apply some under-relaxation to the boundary values. In my implementation in a Lattice Boltzmann code, this step was crucial for the stability of the Maxwell BC. No results could be obtained without under-relaxation.

Do I need to know what LHORI means?
flotus1 is offline   Reply With Quote

Old   July 29, 2013, 17:01
Default
  #9
Member
 
sud
Join Date: May 2013
Posts: 65
Rep Power: 4
archeoptyrx is on a distinguished road
Yes . My code produces good results with no-slip BC . I even compared the plots to Anderson book . When it comes to slip , it just give me wrong values or some weird plots .
archeoptyrx is offline   Reply With Quote

Old   July 29, 2013, 17:09
Default
  #10
Senior Member
 
flotus1's Avatar
 
Alex
Join Date: Jun 2012
Location: Germany
Posts: 1,098
Rep Power: 19
flotus1 will become famous soon enoughflotus1 will become famous soon enough
Forgot to mention: In addition (or if your code is stable instead of) of applying under-relaxation, you can also limit the boundary value to the value at the nearest fluid node.
My preferred sequence for this is first apply the limitation, then the under-relaxation.

Hope your code runs better now. If you still want me to comment on the implementation of the Maxwell boundary, please post the code in text form (no screenshots).
That will make it much easier for me.

Last edited by flotus1; July 30, 2013 at 17:31.
flotus1 is offline   Reply With Quote

Old   August 1, 2013, 15:55
Default
  #11
Member
 
sud
Join Date: May 2013
Posts: 65
Rep Power: 4
archeoptyrx is on a distinguished road
Quote:
Originally Posted by flotus1 View Post
Forgot to mention: In addition (or if your code is stable instead of) of applying under-relaxation, you can also limit the boundary value to the value at the nearest fluid node.
My preferred sequence for this is first apply the limitation, then the under-relaxation.

Hope your code runs better now. If you still want me to comment on the implementation of the Maxwell boundary, please post the code in text form (no screenshots).
That will make it much easier for me.
Mr.Alex,

I am sorry for posting images. I tried uploading the .m files and it didnt allow me .

Have a look at this plot which is obtained after 200 iterations (i used second order accurate for the gradient terms and u=0,Tw=300 as initial wall conditions) . Should i run the simulation for some 1000 iterations and look for the results ..?

The Knudsen number for this simulation is Kn=0.07 . I have the results from openfoam DSMC solver for the same to validate my code results .


Also have a look at the BC. Help me out if i am wrong .
Attached Images
File Type: jpg Pressure.jpg (30.0 KB, 11 views)
Attached Files
File Type: txt BC.txt (2.0 KB, 5 views)
archeoptyrx is offline   Reply With Quote

Old   August 1, 2013, 17:29
Default
  #12
Senior Member
 
flotus1's Avatar
 
Alex
Join Date: Jun 2012
Location: Germany
Posts: 1,098
Rep Power: 19
flotus1 will become famous soon enoughflotus1 will become famous soon enough
I imagined something like this, but at least we have text now
What is Ac in your expression for the slip velocity?

Code:
% surface boundary j=1 with Maxwell Slip Boundary 1st order
for i=2:imax
% for simplicity and to eliminate possible sources of error, I would set a fixed mean free path for now
%   Meanfree(i,1)=(mu(i,1)/rho(i,1))*sqrt(3.14/(2*R*T(i,1)));
      Meanfree(i,1) = 0.07;
% I cant recommend going to second order as long as we are still facing stability issues. Lets start with first order accurate derivatives.
%   du=(-u(i,3)+4*u(i,2)-3*u(i,1))/(2*dely);
%   dT=(-T(i,3)+4*T(i,2)-3*T(i,1))/(2*dely);
du = (u(i,2) - u(i,1))/dely;
dT = (T(i,2) - T(i,1))/dely;
% still no relaxation and limiters...
u_temp = u(i,1);
     u(i,1)=(2-sigV)*Meanfree(i,1)*du/sigV;
% limiter
if (abs(u(i,1)) > abs(u(i,2)))
  u(i,1) = u(i,2);
end
% under-relaxation, typical value for urf between 0.01 and 0.1
u(i,1) = urf * u(i,1) + (1.0-urf) * u_temp;
% apply the same for the temperature
    T(i,1)=Tw+(2-sigT)*2*gamma*Meanfree(i,1)*dT/(sigT*Pr*(gamma+1));
%     u(i,1)=Us(i,1);
    v(i,1)=0;
%     T(i,1)=Tw;
    P(i,1)=2*P(i,2)-P(i,3);
    rho(i,1)=P(i,1)/(R*T(i,1));
    k(i,1)=Thermc(mu(i,1));
    e(i,1)=cv*T(i,1);
    Et(i,1)=rho(i,1)*(e(i,1)+0.5*((u(i,1)^2)+(v(i,1)^2)));
end

Last edited by flotus1; August 2, 2013 at 02:27.
flotus1 is offline   Reply With Quote

Old   August 1, 2013, 17:42
Default
  #13
Member
 
sud
Join Date: May 2013
Posts: 65
Rep Power: 4
archeoptyrx is on a distinguished road
Ac =0.8 . Like i am taking only 80% of the calculated values . If its not needed i ll remove it . I modified it . Let me try and post the results ..

Thanks
archep
archeoptyrx is offline   Reply With Quote

Old   August 1, 2013, 18:08
Default
  #14
Senior Member
 
flotus1's Avatar
 
Alex
Join Date: Jun 2012
Location: Germany
Posts: 1,098
Rep Power: 19
flotus1 will become famous soon enoughflotus1 will become famous soon enough
In this case, remove Ac from the formula.
flotus1 is offline   Reply With Quote

Old   August 1, 2013, 23:40
Default
  #15
Member
 
sud
Join Date: May 2013
Posts: 65
Rep Power: 4
archeoptyrx is on a distinguished road
Mr.Alex ,

I have deleted the Ac in the program. Now it showed me mass flow discontinuity after 5 iterations . Kindly check for the convergence criteria and Mass conservation file attached .May be bcz of that the meanfree path of 0.07 , to the value that produces slip .

should I have to change something in the files attached ?
Attached Files
File Type: txt Convergence.txt (194 Bytes, 3 views)
File Type: txt Mdot.txt (335 Bytes, 3 views)
archeoptyrx is offline   Reply With Quote

Old   August 2, 2013, 02:44
Default
  #16
Senior Member
 
flotus1's Avatar
 
Alex
Join Date: Jun 2012
Location: Germany
Posts: 1,098
Rep Power: 19
flotus1 will become famous soon enoughflotus1 will become famous soon enough
Please post your functions with [CODE] tags.

What is this supposed to do?
As soon as the absolute difference for the density is larger than 1e-8 in any element in any time step, the solution is considered diverged?
I disagree on such a criterion for divergence.

Code:
 Convergence=1;
for i=1:imax
    for j=1:jmax
        if(abs(rho(i,j)-oldrho(i,j))>10^(-8))
            % Solution is not converged
            Convergence=0;
        end
    end
end

Again, I would not subscribe to such a criterion for "discontinuity".
This function simply checks if the mass flow rate at the outlet is higher than at the inlet. If this is the case, you say it is a mass discontinuity.
During the iterative solution process, such conditions can not be avoided.

Code:
function Verf=Mdot(rho,u)
i=1;
ML=0;
global imax jmax ;
for j=1:jmax
    ML=ML+rho(i,j)*u(i,j);
end
i=imax;
MR=0;
for j=1:jmax
    MR=MR+rho(i,j)*u(i,j);
end

diff=abs(ML-MR)*100/ML;
if diff>1
%     disp('Mass flow discontinuity')
    Verf=0;
else
%     disp('Mass flow is conserved')
    Verf=1;
end
end
In my opinion, none of the above functions helps to determine how the solution is proceeding.
Did you really stop iterating when the warning messages from these functions showed up?
flotus1 is offline   Reply With Quote

Old   August 2, 2013, 03:23
Default
  #17
Member
 
sud
Join Date: May 2013
Posts: 65
Rep Power: 4
archeoptyrx is on a distinguished road
Quote:
Originally Posted by flotus1 View Post
Please post your functions with [CODE] tags.

What is this supposed to do?
As soon as the absolute difference for the density is larger than 1e-8 in any element in any time step, the solution is considered diverged?
I disagree on such a criterion for divergence.

Code:
 Convergence=1;
for i=1:imax
    for j=1:jmax
        if(abs(rho(i,j)-oldrho(i,j))>10^(-8))
            % Solution is not converged
            Convergence=0;
        end
    end
end
Mr.Alex ,

Actually this code of convergence is the one provided in book . Now the thing is it actually worked for no slip conditions producing exactly same results as given in the book . The book says "The solution is considered as converged when the density at each grid point changes no more than 10^-8 between time steps ."
So in each time step , the main program sets time step , does maccormick calculations , then calls this convergence function , and if the diff btween previous time step is less than the prescribed one , it will stop , else it will exchange the values of new rho to oldrho and goes for second time step .

Code:
 
%%MAin program (Excluding the global declarations /free stream conditions)
 if iter==1
        %Initialize the flow field variables at t=0 (ie at x=0)
        for i =1:imax
            for j=1:jmax
                u(i,j)=Vinf;
                v(i,j)=0;
                P(i,j)=Pinf;
                T(i,j)=Tinf;
            end
        end
        
        % Free stream boundary j=jmax
        for i =1:imax
            u(i,jmax)=Vinf;
            v(i,jmax)=0;
            P(i,jmax)=Pinf;
            T(i,jmax)=Tinf;
        end
        % Inflow boundary condition i==1
        for j=2:jmax-1
            u(1,j)=Vinf;
            v(1,j)=0;
            P(1,j)=Pinf;
            T(1,j)=Tinf;
        end
        u(1,1)=0;
        v(1,1)=0;
        P(1,1)=Pinf;
        T(1,1)=Tinf;
        % surface boundary j=1
        for i=2:imax
            u(i,1)=0;
            v(i,1)=0;
            T(i,1)=Tw;
            P(i,1)=2*P(i,2)-P(i,3);
        end
        % Outflow boundary
        for j=2:jmax-1
            u(imax,j)=2*u(imax-1,j)-u(imax-2,j);
            v(imax,j)=2*v(imax-1,j)-v(imax-2,j);
            P(imax,j)=2*P(imax-1,j)-P(imax-2,j);
            T(imax,j)=2*T(imax-1,j)-T(imax-2,j);
        end
    end
    iter     % display the iteration number
    delt(iter)=Tstep();                                    % Call the Time step for the iteration
    
    MAC(.15*delt(iter));       % Call the Maccormick Technique
    
    Convergence=Conver(rho,oldrho);                      % Call the Convergence Function to check convergence
    
    if Convergence==1
        break;
    end
    oldrho=rho;
    
end


Verf=Mdot(rho,u);
if Verf==0
    fprintf('\n MAss flow discontinuity !!!!! :O\n')
    break
end
I posted the convergence file and Mdot file already . Now for defining the Mdot file , the book says " the deviation between the mass flow rate at entrance and exit is less than 1%" . so i took that way . Am i wrong ? , May be i am wrong in the convergence criteria bcz for the same problem , the solution presented in book get converged for 4339 iterations and for me it is some 10000 iterations .
archeoptyrx is offline   Reply With Quote

Old   August 2, 2013, 03:33
Default
  #18
Senior Member
 
flotus1's Avatar
 
Alex
Join Date: Jun 2012
Location: Germany
Posts: 1,098
Rep Power: 19
flotus1 will become famous soon enoughflotus1 will become famous soon enough
Ok, now I get the point for the convergence criterion with the mass flow discontinuity. I actually misinterpreted the code.
But why do you break the solution process as soon as the mass imbalance is larger than 1% between inlet and outlet?
flotus1 is offline   Reply With Quote

Old   August 2, 2013, 03:38
Default
  #19
Member
 
sud
Join Date: May 2013
Posts: 65
Rep Power: 4
archeoptyrx is on a distinguished road
No i am not breaking . Till complete convergence , the program runs , after the convergence the mass imbalance is checked to say whether the mass is conserved at inlet n exit . So now I tried with what you said . I ended up with errors for fixing the mean free path . I tried to use the formula for mu + limiter ----> errors it stopped at 45 iterations

What should i do further ?
But the way i implemented the slip and jump conditions are correct right ??
archeoptyrx is offline   Reply With Quote

Old   August 2, 2013, 04:39
Default
  #20
Senior Member
 
flotus1's Avatar
 
Alex
Join Date: Jun 2012
Location: Germany
Posts: 1,098
Rep Power: 19
flotus1 will become famous soon enoughflotus1 will become famous soon enough
I am running out of ideas.
The way the velocity and temperature boundary conditions are implemented should be correct, provided your adjustments for the temperature are correct.

What happens if you set the slip velocity and temperature to constant values (but not the actual values of the wall)
flotus1 is offline   Reply With Quote

Reply

Thread Tools
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are On
Pingbacks are On
Refbacks are On



All times are GMT -4. The time now is 23:52.