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

Writing a CFD code

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

Like Tree1Likes
  • 1 Post By vatavuk

Reply
 
LinkBack Thread Tools Display Modes
Old   May 26, 2011, 06:06
Default Writing a CFD code
  #1
New Member
 
Dave Smith
Join Date: Jul 2010
Posts: 27
Rep Power: 7
davesmith_01 is on a distinguished road
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
davesmith_01 is offline   Reply With Quote

Old   May 27, 2011, 04:39
Default
  #2
New Member
 
Join Date: May 2011
Posts: 6
Rep Power: 6
Raambi is on a distinguished road
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
Raambi is offline   Reply With Quote

Old   May 27, 2011, 06:46
Default
  #3
New Member
 
Dave Smith
Join Date: Jul 2010
Posts: 27
Rep Power: 7
davesmith_01 is on a distinguished road
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?
davesmith_01 is offline   Reply With Quote

Old   May 27, 2011, 09:26
Default
  #4
Senior Member
 
Paulo Vatavuk
Join Date: Mar 2009
Location: Campinas, Brasil
Posts: 111
Rep Power: 8
vatavuk is on a distinguished road
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

Last edited by vatavuk; May 29, 2011 at 15:34. Reason: correct a small error
vatavuk is offline   Reply With Quote

Old   June 2, 2011, 14:18
Default
  #5
New Member
 
Anonymous
Join Date: Jul 2010
Posts: 3
Rep Power: 7
LESter is on a distinguished road
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.
LESter is offline   Reply With Quote

Old   June 3, 2011, 00:21
Default
  #6
New Member
 
Felipe Hernandez
Join Date: Apr 2011
Posts: 13
Rep Power: 6
Petrov is on a distinguished road
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!
Petrov is offline   Reply With Quote

Old   September 14, 2012, 14:13
Default How to run the peric's code
  #7
New Member
 
Join Date: Dec 2011
Posts: 12
Rep Power: 5
adampaya is on a distinguished road
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.
adampaya is offline   Reply With Quote

Old   September 15, 2012, 07:35
Default
  #8
Senior Member
 
Paulo Vatavuk
Join Date: Mar 2009
Location: Campinas, Brasil
Posts: 111
Rep Power: 8
vatavuk is on a distinguished road
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 likes this.
vatavuk is offline   Reply With Quote

Old   September 16, 2012, 03:54
Default
  #9
New Member
 
Join Date: Dec 2011
Posts: 12
Rep Power: 5
adampaya is on a distinguished road
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???
adampaya is offline   Reply With Quote

Old   September 17, 2012, 16:48
Default
  #10
Senior Member
 
Paulo Vatavuk
Join Date: Mar 2009
Location: Campinas, Brasil
Posts: 111
Rep Power: 8
vatavuk is on a distinguished road
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
vatavuk is offline   Reply With Quote

Old   September 22, 2012, 11:46
Default
  #11
New Member
 
Join Date: Dec 2011
Posts: 12
Rep Power: 5
adampaya is on a distinguished road
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.
adampaya is offline   Reply With Quote

Old   September 23, 2012, 01:41
Default
  #12
Senior Member
 
Cean
Join Date: Feb 2010
Posts: 128
Rep Power: 7
shirazbj is on a distinguished road
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.
shirazbj is offline   Reply With Quote

Old   September 26, 2012, 10:31
Default
  #13
New Member
 
Join Date: Dec 2011
Posts: 12
Rep Power: 5
adampaya is on a distinguished road
The error numbers are LNK2001 and LNK1120.
I don't think that the problem would be related to the version of Fortran.
adampaya is offline   Reply With Quote

Old   October 6, 2012, 04:25
Default
  #14
New Member
 
nas
Join Date: Apr 2012
Posts: 2
Rep Power: 0
yeganeh is on a distinguished road
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
yeganeh 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


Similar Threads
Thread Thread Starter Forum Replies Last Post
writing a cfd code kelvin Main CFD Forum 13 October 13, 2005 05:41
A simple CFD code for teaching basic CFD? Christoph Lund Main CFD Forum 13 September 14, 2005 04:36
CFD JOBS and Expected Salary.... Noel Harrison Main CFD Forum 11 November 22, 2000 08:15
CFD Salary CFD Main CFD Forum 15 September 4, 1999 14:04
Commercial CFD code Hanson G. He Main CFD Forum 1 October 15, 1998 08:49


All times are GMT -4. The time now is 13:58.