|
[Sponsors] |
May 26, 2011, 07:06 |
Writing a CFD code
|
#1 |
New Member
Dave Smith
Join Date: Jul 2010
Posts: 27
Rep Power: 16 |
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 |
|
May 27, 2011, 05:39 |
|
#2 |
New Member
Join Date: May 2011
Posts: 6
Rep Power: 15 |
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 |
|
May 27, 2011, 07:46 |
|
#3 |
New Member
Dave Smith
Join Date: Jul 2010
Posts: 27
Rep Power: 16 |
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? |
|
May 27, 2011, 10:26 |
|
#4 |
Senior Member
Paulo Vatavuk
Join Date: Mar 2009
Location: Campinas, Brasil
Posts: 200
Rep Power: 18 |
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 16:34. Reason: correct a small error |
|
June 2, 2011, 15:18 |
|
#5 |
New Member
Anonymous
Join Date: Jul 2010
Posts: 3
Rep Power: 16 |
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.
|
|
June 3, 2011, 01:21 |
|
#6 |
New Member
Felipe Hernandez
Join Date: Apr 2011
Posts: 13
Rep Power: 15 |
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! |
|
September 14, 2012, 15:13 |
How to run the peric's code
|
#7 |
New Member
Join Date: Dec 2011
Posts: 12
Rep Power: 14 |
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. |
|
September 15, 2012, 08:35 |
|
#8 |
Senior Member
Paulo Vatavuk
Join Date: Mar 2009
Location: Campinas, Brasil
Posts: 200
Rep Power: 18 |
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 |
|
September 16, 2012, 04:54 |
|
#9 |
New Member
Join Date: Dec 2011
Posts: 12
Rep Power: 14 |
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??? |
|
September 17, 2012, 17:48 |
|
#10 |
Senior Member
Paulo Vatavuk
Join Date: Mar 2009
Location: Campinas, Brasil
Posts: 200
Rep Power: 18 |
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 |
|
September 22, 2012, 12:46 |
|
#11 |
New Member
Join Date: Dec 2011
Posts: 12
Rep Power: 14 |
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. |
|
September 23, 2012, 02:41 |
|
#12 |
Senior Member
Cean
Join Date: Feb 2010
Posts: 128
Rep Power: 16 |
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. |
|
September 26, 2012, 11:31 |
|
#13 |
New Member
Join Date: Dec 2011
Posts: 12
Rep Power: 14 |
The error numbers are LNK2001 and LNK1120.
I don't think that the problem would be related to the version of Fortran. |
|
October 6, 2012, 05:25 |
|
#14 |
New Member
nas
Join Date: Apr 2012
Posts: 2
Rep Power: 0 |
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 |
|
|
|
Similar Threads | ||||
Thread | Thread Starter | Forum | Replies | Last Post |
CFD Salary | CFD | Main CFD Forum | 17 | January 3, 2017 18:09 |
writing a cfd code | kelvin | Main CFD Forum | 13 | October 13, 2005 06:41 |
A simple CFD code for teaching basic CFD? | Christoph Lund | Main CFD Forum | 13 | September 14, 2005 05:36 |
CFD JOBS and Expected Salary.... | Noel Harrison | Main CFD Forum | 11 | November 22, 2000 08:15 |
Commercial CFD code | Hanson G. He | Main CFD Forum | 1 | October 15, 1998 09:49 |