[davesmith_01]

Hi

I have always wanted to write my own 3D Navier stokes code for subsonic incompressible flow, I know this is a very complex task and takes a long time. The problem I am not sure where to start as I do not have any example codes to have a look at for a general guide and then make my own code once I have seen a structure of how the code should be written. Can anyone please help by providing some 2D and 3D N-S codes which I can use as guide

Thanks, Dave

[Raambi]

Hello!

I learned a lot from the tutorials on this page:

http://www.cfd.tu-berlin.de/index.ph...n&lang=english

U can load the pdf-files without an password. But it is in german and the standard solution codes are not available at everytime :/

If u have done the tutorials from CFD1 and CFD2 u should have an overview of the techniques.

After that u can take a look on the Code of OpenFoam and change it if u want.

regards

[davesmith_01]

Thanks friend

I will convert some of it and see if it helps.

Do you have any other good tutorials helping to solve the N-S equations and coding them?

[vatavuk]

Hi Dave,
If you know FORTRAN, I would suggest that you study the code that comes with the book of Ferziger and Peric "Computational Methods for Fluid Dynamics". The code can be downloaded at ftp://ftp.springer.de/pub/technik/peric, it’s written in a clear way and the algorithms are easy to understand because they are explained in the book.
There is also the Dolfyn code, which is more advanced and is also based in the book of Ferziger and Peric.
Paulo

[LESter]

If you are truly starting a code from scratch, start small and keep it simple. I would suggest starting with the 2D Euler equations and slowly build from there. Once you have a basic structure adding additional pieces and dimensions is much easier. If you are really new to this, you might want to start playing around with some model problems like the 1D advection equation, using the type of method you are interested in.

[Petrov]

Here's how I did it:
First, I used this book: http://www.amazon.com/Finite-Element.../dp/0486411818 to write a code for 2D heat conduction.

Then, I used this http://ta.twi.tudelft.nl/users/vuik/.../fem_notes.pdf to learn how to apply those techniques to the Navier-Stokes equations (again in 2D because I'm lazy, ok).

This is for using the FEM. I couldn't find good (and cheap) resources for the DG method.

Oh, and here's a tip that solved weird bugs twice: always make sure you take the absolute value of the jacobian.

Good luck!

[adampaya]

Hello everyone
I have some challenge in using the peric code
For example consider 2DGL, when I generate the gird for the geometry, what's the next step??? Should I use the Caffa.f or user.f ????
I read the read me file, but I don't understand.

Thanks for your help.

[vatavuk]

Hi adampaya,
If you look at the end of the caffa.f file, you'll see the command INCLUDE 'user.f', which copies the contents of user.f into the current file. The user.f file is available for the programmer to include special instructions depending on the case you are solving. You can create your own user.f file, modifying the user.gen file or use the user.cav and user.cyl files. To use one of them you must rename the file you want to user.f.
So the answer to your question is: you should copile and run the caffa file.
I hope this helps
Paulo

[adampaya]

Hi Paulo
yes you're right. But for instance consider the example of "cav30l".
the grid can be generated by Grid.f using the instructions gave in the cav30l.gin in the examples directory.
but when I'm trying to run the user.f code after generating the grid,
it gives error, and is not run. what is the problem???

[vatavuk]

Hi adampaya,
A few information could help to understand the problem:
The Caffa.f file was copiled without problems? What is the error message that you receive?
Regards,
Paulo

[adampaya]

C################################################# ########
SUBROUTINE BCIN(K)
C################################################# ########
C This routine sets inlet boundary conditions (which
C may also change in time in an unsteady problem; in
C that case, set the logical variable INIBC to 'true'.)
C It is called once at the begining of calculations on
C each grid level, unles INIBC=.TRUE. forces calls in
C each time step.
C This routine includes several options which are most
C often required and serves the test cases provided with
C the code; it may require adaptation to the problem
C being solved, especially regarding velocity or
C temperature distribution at boundaries.
C================================================= ========
INCLUDE 'param.inc'
INCLUDE 'indexc.inc'
INCLUDE 'rcont.inc'
INCLUDE 'geo.inc'
INCLUDE 'var.inc'
INCLUDE 'bound.inc'
INCLUDE 'logic.inc'
C
C.....SET INDICES, INITIALIZE MOMENTUM AND MASS FLOW RATE
C
INIBC=.FALSE.
CALL SETIND(K)
FLOMOM=0.
FLOMAS=0.
C
C.....SET HOT AND COLD WALL TEMPERATURE (WEST COLD, EAST HOT)
C
IF(LCAL(IEN)) THEN
DO IW=IWS(K)+1,IWS(K)+NWALI(K)
IJ=IJW(IW)
IF((IJW1(IW)-IJW2(IW)).LT.0) THEN
T(IJ)=TC
ELSE
T(IJ)=TH
ENDIF
END DO
ENDIF
C
C.....SET ADIABATIC WALL TEMPERATURE TO REFERENCE TEMPERATURE
C
IF(LCAL(IEN)) THEN
DO IW=IWAS(K)+1,IWAS(K)+NWALA(K)
IJ=IJW(IW)
T(IJ)=TREF
END DO
ENDIF
C
C.....SET LID VELOCITY (FOR LID-DRIVEN CAVITIES ONLY; SET ULID=0. FOR
C BUOYANCY_DRIVEN CAVITY FLOWS IN THE *.CIN FILE)
C
DO I=2,NIM
IJ=LI(I+IST)+NJ
U(IJ)=ULID
END DO
C
C.....SET RESIDUAL NORMALISATION FACTORS
C
DO L=1,NFI
RNOR(L)=1.
RESOR(L)=0.0
END DO
C
IF(FLOMAS.LT.SMALL) FLOMAS=1.
IF(FLOMOM.LT.SMALL) FLOMOM=1.
C
RNOR(IU)=1./(FLOMOM+SMALL)
RNOR(IV)=RNOR(IU)
RNOR(IP)=1./(FLOMAS+SMALL)
RNOR(IEN)=1./((FLOMAS*TREF)+SMALL)
C
RETURN
END
C
C
C################################################# #######
SUBROUTINE SOUT(K)
C################################################# #######
C This routine prints some additional quantities, as
C programmed by the user for the problem considered.
C Note that profiles can be obtained in post-processor.
C================================================= =======
INCLUDE 'param.inc'
INCLUDE 'indexc.inc'
INCLUDE 'rcont.inc'
INCLUDE 'var.inc'
INCLUDE 'geo.inc'
INCLUDE 'bound.inc'
INCLUDE 'logic.inc'
C
CALL SETIND(K)
C
C.....Calculate total shear force on a part of wall boundary,
C.....and the distribution of the shear stress. Here: bottom wall.
C
FTAUX=0.
FTAUY=0.
WRITE(2,*) ' '
WRITE(2,*) ' XC YC TAUX TAUY '
C
DO IW=IWS(K)+1,IWS(K)+NWAL(K)
IF(YC(IJW(IW)).LT.1.E-5) THEN
IJP=IJPW(IW)
IJB=IJW(IW)
COEF=VIS(IJB)*SRDW(IW)
TAUX=COEF*((U(IJP)-U(IJB))*XTW(IW)**2+
* (V(IJP)-V(IJB))*XTW(IW)*YTW(IW))
TAUY=COEF*((U(IJP)-U(IJB))*XTW(IW)*YTW(IW)+
* (V(IJP)-V(IJB))*YTW(IW)**2)
S=SQRT((X(IJW1(IW))-X(IJW2(IW)))**2+
* (Y(IJW1(IW))-Y(IJW2(IW)))**2)
IF(LAXIS) S=S*YC(IJB)
FTAUX=FTAUX+TAUX
FTAUY=FTAUY+TAUY
WRITE(2,'(1P4E14.4)') XC(IJB),YC(IJB),TAUX/S,TAUY/S
ENDIF
END DO
WRITE(2,*) ' '
WRITE(2,*) ' TOTAL SHEAR FORCE IN X-DIRECTION: ',FTAUX
WRITE(2,*) ' TOTAL SHEAR FORCE IN Y-DIRECTION: ',FTAUY
WRITE(2,*) ' '
C
C.....Calculate total pressure force on a part of wall boundary,
C.....and the distribution of the pressure. Here: bottom wall.
C.....Note: this is the force exerted by fluid onto wall.
C
FPRX=0.
FPRY=0.
WRITE(2,*) ' '
WRITE(2,*) ' XC YC PRESSURE '
C
DO IW=IWS(K)+1,IWS(K)+NWAL(K)
IF(YC(IJW(IW)).LT.1.E-5) THEN
IJB=IJW(IW)
SX=Y(IJW1(IW))-Y(IJW2(IW))
SY=X(IJW1(IW))-X(IJW2(IW))
IF(LAXIS) SX=SX*YC(IJB)
IF(LAXIS) SY=SY*YC(IJB)
FPRX=FPRX+P(IJB)*SX
FPRY=FPRY-P(IJB)*SY
WRITE(2,'(1P3E14.4)') XC(IJB),YC(IJB),P(IJB)
ENDIF
END DO
WRITE(2,*) ' '
WRITE(2,*) ' TOTAL PRESSURE FORCE IN X-DIRECTION: ',FPRX
WRITE(2,*) ' TOTAL PRESSURE FORCE IN Y-DIRECTION: ',FPRY
WRITE(2,*) ' '
C
C.....CALCULATE HEAT FLUX THROUGH AN ISOTHERMAL WALL (here HOT)
C.....HEAT(IW) is the specific heat flux (per unit area); HEATS is
C.....the total flux through the wall surface.
C
IF(LCAL(IEN)) THEN
WRITE(2,*) ' '
WRITE(2,*) ' XC YC HEAT FLUX '
HEATS=0.
DO IW=IWS(K)+1,IWS(K)+NWALI(K)
IF(T(IJW(IW)).GT.TREF) THEN
HEAT=(VISC/PRANL)*SRDW(IW)*(T(IJW(IW))-T(IJPW(IW)))
HEATS=HEATS+HEAT
S=SQRT((X(IJW1(IW))-X(IJW2(IW)))**2+
* (Y(IJW1(IW))-Y(IJW2(IW)))**2)
IF(LAXIS) S=S*YC(IJW(IW))
WRITE(2,'(1P3E14.4)') XC(IJW(IW)),YC(IJW(IW)),HEAT/S
ENDIF
END DO
WRITE(2,*) ' '
WRITE(2,*) ' TOTAL HEAT FLUX THROUGH THE HOT WALL: ',HEATS
WRITE(2,*) ' '
ENDIF
C
RETURN
END
C
C
C################################################# #######
SUBROUTINE TOUT(K)
C################################################# #######
C This routine prints some additional quantities, as
C programmed by the user for the problem considered.
C Note that profiles can be obtained in post-processor.
C================================================= =======
INCLUDE 'param.inc'
INCLUDE 'indexc.inc'
INCLUDE 'rcont.inc'
INCLUDE 'var.inc'
INCLUDE 'varold.inc'
INCLUDE 'geo.inc'
INCLUDE 'bound.inc'
INCLUDE 'logic.inc'
C
C
RETURN
END





Hi Paulo
I use this code (user.f) after I made the grid of the cavity (for example).
but it doesn't run and gives error. what changes should I do to get the code run?

Thanks.

[shirazbj]

Hi adampaya,

What compiler are you using and under what OS?

Since the code was written for old F77, it may not working for new complier like today's gfortran.

for example, this sentence in grid.f may not work:

NOT=MAX(NOT,1)

since NOT is reserved by the new complier.

[adampaya]

The error numbers are LNK2001 and LNK1120.
I don't think that the problem would be related to the version of Fortran.

[yeganeh]

hi
i am looking for fortran code for solution energy and navier stokes equations by Simple method for numerical study in triangular lid driven ?
thanks