CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Programming & Development

How to interface a Fortran thermodynamic tool with OpenFOAM?

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

Like Tree4Likes

Reply
 
LinkBack Thread Tools Display Modes
Old   October 29, 2010, 09:18
Default How to interface a Fortran thermodynamic tool with OpenFOAM?
  #1
Cyp
Senior Member
 
Cyprien
Join Date: Feb 2010
Location: Stanford University
Posts: 240
Rep Power: 9
Cyp is on a distinguished road
Hi everybody!

I am looking to interface a specific thermodynamic tool with OpenFOAM. This program is written in Fortran. From a Pressure, a temperature and a composition it can provide a density : rho=f(P,T,x)

I wonder how to bound these different programs. I would like to keep the OpenFOAM structure and evaluate the density at each time step with a command as : thermo.rho();

Is anybody already performed such an interfacing ? Can you give some hints to start ? from which thermophysicalModels can I be inspired?

Regards,
Cyp
Cyp is offline   Reply With Quote

Old   November 1, 2010, 06:13
Default
  #2
Senior Member
 
Christian Lucas
Join Date: Aug 2009
Location: Braunschweig, Germany
Posts: 199
Rep Power: 8
Chris Lucas is on a distinguished road
Hi,

compile the fortran programm as a dynamic library. Then, use that library in OpenFoam.

Is the programm you want the use in OpenFoam Refprop (NIST)?

Regards,
Christian
Chris Lucas is offline   Reply With Quote

Old   November 1, 2010, 11:44
Default
  #3
Cyp
Senior Member
 
Cyprien
Join Date: Feb 2010
Location: Stanford University
Posts: 240
Rep Power: 9
Cyp is on a distinguished road
Hi Chris!

Thank you for your answer!! Indeed it is exactly the NIST REFPROP I want to use with OpenFOAM! Did you already experienced to call REFPROP from OpenFOAM ?


For the moment, I compiled the Fortran subroutine as a shared object (.so). How can I call it from OF ??

Best regards,
Cyp
Cyp is offline   Reply With Quote

Old   November 1, 2010, 11:53
Default
  #4
Senior Member
 
Christian Lucas
Join Date: Aug 2009
Location: Braunschweig, Germany
Posts: 199
Rep Power: 8
Chris Lucas is on a distinguished road
Hi,

yes, I have implemented Refprop in OpenFoam. However, I noticed that Refprop is far too slow to be used in OpenFOAM. The simulation time will increase in order of many magnitudes.

Regards,
Christian
Chris Lucas is offline   Reply With Quote

Old   November 1, 2010, 12:03
Default
  #5
Member
 
kshitij neroorkar
Join Date: Mar 2009
Location: Michigan, USA
Posts: 32
Rep Power: 8
kdneroorkar is on a distinguished road
Hi
another way to do this (we use this with REFPROP) is to create a lookup table using REFPROP with rho as a function of (P,T,x). Then write a class called, say, thermoprops. when you call thermo.rho(P,T,x) (where thermo is of type thermoprops) , it would interpolate in this table and return the required rho.
Hope this helps
Kshitij
kdneroorkar is offline   Reply With Quote

Old   November 2, 2010, 05:18
Default
  #6
Cyp
Senior Member
 
Cyprien
Join Date: Feb 2010
Location: Stanford University
Posts: 240
Rep Power: 9
Cyp is on a distinguished road
Hi!

Quote:
Hi,

yes, I have implemented Refprop in OpenFoam. However, I noticed that Refprop is far too slow to be used in OpenFOAM. The simulation time will increase in order of many magnitudes.

Regards,
Christian
It is not really a good news....

Quote:
Originally Posted by kdneroorkar View Post
Hi
another way to do this (we use this with REFPROP) is to create a lookup table using REFPROP with rho as a function of (P,T,x). Then write a class called, say, thermoprops. when you call thermo.rho(P,T,x) (where thermo is of type thermoprops) , it would interpolate in this table and return the required rho.
Hope this helps
Kshitij
Do you mean I have to generate a table from REFPROP and then pick up the right value of rho in this table from P, T and x ?

I think my main problem is either I chose a direct call of REFPROP or an interpolation from a table, is the creation of the thermoprops class... From which existed class can I start ??

Regards,
Cyp
Cyp is offline   Reply With Quote

Old   November 2, 2010, 09:47
Default
  #7
Member
 
kshitij neroorkar
Join Date: Mar 2009
Location: Michigan, USA
Posts: 32
Rep Power: 8
kdneroorkar is on a distinguished road
yes. so you would have something like
#P T x rho
1000 273 0.01 900
1200 273 0.01 890
.....
and you would interpolate to get rho=f(P,T,x). another idea you could try ( I dont know if this is feasible in your case), instead of using P, T and x, you could use P,H. This way the interpolation would be simplified.

I dont know of any existing class that does this interpolation. We had to write our own.
kdneroorkar is offline   Reply With Quote

Old   November 2, 2010, 10:37
Default
  #8
Cyp
Senior Member
 
Cyprien
Join Date: Feb 2010
Location: Stanford University
Posts: 240
Rep Power: 9
Cyp is on a distinguished road
Thank you for your answers! For the moment, I would like to test to implement a direct call from REFPROP instead of a reading in a table.

As a first test, and to become familiar with the thermodynamic class programmation in OpenFOAM, I try to create my own class from basicRhoThermo : rho=f(T,P). However, when I read the .C and .H files, I can't where rho is calculated...

Do you have any idea ?

Regards,
Cyp
Cyp is offline   Reply With Quote

Old   November 2, 2010, 11:38
Default
  #9
Member
 
kshitij neroorkar
Join Date: Mar 2009
Location: Michigan, USA
Posts: 32
Rep Power: 8
kdneroorkar is on a distinguished road
basicRhoThermo is a base class. since rho() is a virtual function, it would be calculated in some class that is derived from basicRhoThermo. If you look at some class that is derived from basicRhoThermo, you might get an idea.
I am still using 1.5-dev so I dont have basicRhoThermo. Sorry that I cant be of more help
kdneroorkar is offline   Reply With Quote

Old   November 2, 2010, 16:55
Default
  #10
Senior Member
 
Christian Lucas
Join Date: Aug 2009
Location: Braunschweig, Germany
Posts: 199
Rep Power: 8
Chris Lucas is on a distinguished road
Hi,

might I ask, what exactly you want to simulated?

Ok, here is what you have to do to implement refprop in OpenFoam. This is not the best way to do it, but it works.

Start from the file hPsiThermo (or hRhoThermo). In this file, you have to call the Refprop functions. In the constuctor, you have to create the Refprop object. You also need a second "rho" field. This is needed because rho not directly connected to psi anymore (as it is for perfect gases). In this class, all the functions have to be changed so that the Refprop functions are called.

Now comes the tricky part. Firstly, you have to change the enthalpy BC's. The "normal" BC are programmed for perfect gases, so we have h(T) instead of h(T,p). You also need to change the functions in basicThermo called by the enthalpy BC.



As I said, Refprop is very slow and you have to solve the fundamental equations for each grid point (without changing refprop even more than once).

About using tables, it is a possibility. But one problem is that you get an error and I don't mean the "normal" interpolation error. The problem is that thermo physical properties are connected to each other e.g. c_p=dh/dT. So, the interpolated value of c_p might not correspond correctly with you enthalpy.


Regards,
Christian
mgg likes this.

Last edited by Chris Lucas; November 3, 2010 at 04:42.
Chris Lucas is offline   Reply With Quote

Old   November 3, 2010, 06:12
Default
  #11
Cyp
Senior Member
 
Cyprien
Join Date: Feb 2010
Location: Stanford University
Posts: 240
Rep Power: 9
Cyp is on a distinguished road
Hi Christian!

Thank you for your extended answer! In my investigation fields, the final step will be to simulate a two-phase multicomponent flow. I already have developed models for the mass fractions in each phase and I need to obtain the density of each phase from REFPROP. For each phase, density depends on Temperature, Pressure and Composition. It is my final objectives. I will consider a rough mesh. So I hope that the use of REFPROP will not have a consequent impact on my CPU time...

For the moment, to become familiar with OpenFoam/REFPROP interfacing, I tried the implementation on a simpler model : just one phase is considered. There is no question of enthalpy, Cp and other thermodynamic properties. Just the density (I call the TPRHO() Refprop subroutine).

As you advised me, I started from hPsiThermo. I saw enthalpy, Cp and Cv variables. Since I not consider such variables, I suppose I can rid of them ? And so I can rid of the part where I should change the enthalpy BC, can't I ?

I am not sure to well understand the part with the second density. Do you mean I have to create a rho2 variable in hPsiThermo.H file ?

Moreover, if I understand well, it is in the hPsiThermo.C file that I call the REFPROP function ?

Thank you very much for your help.
Cyp is offline   Reply With Quote

Old   November 3, 2010, 06:48
Default
  #12
Senior Member
 
Christian Lucas
Join Date: Aug 2009
Location: Braunschweig, Germany
Posts: 199
Rep Power: 8
Chris Lucas is on a distinguished road
Hi,

I'm also interested in multiphase flow. So you have phase change in your application? How did you get your pressure equation?

But back to your question:

"As you advised me, I started from hPsiThermo. I saw enthalpy, Cp and Cv variable. Since I not consider such variables, I suppose I can rid of them ?"

For now, yes. But for a real compressible flow, no, you need this functions like
"template<class MixtureType>
Foam::tmp<Foam::scalarField> Foam::hPsiThermo<MixtureType>::Cp"
even so the variables given to the function have to be changed. Within, this function, you must call the refprop routines.

"I am not sure to well understand the part with the second density. Do you mean I have to create a rho2 variable in hPsiThermo.H file?"

About the second density field. If you look at the rho function in basicPsiThermo.C
( virtual tmp<volScalarField> rho() const
{
return p_*psi();
})

You see that rho is connected to psi for perfect gases, so that no rho field has to be defined to save the results for rho. For real fluids, there is no relationship as the one above between rho and psi, so you have the save the results for rho separately.

"And so I can rid of the part where I should change the enthalpy BC, can't I ?"

If you don't have a energy equation, yes, you bypass the BC problem. However this is not correct for compressible flows
This BC functions are in different classes (within a subsubfolder of the basic folder). This classes call functions in basicThermo, which are later overwriten by the functions in your new "hPsiThermo".

"(I call the TPRHO() Refprop subroutine)".

Unless you rewrite the energy equations (I would not use energy equations based on temperature for real fluids), you need a function like rho(p,h) and rho(p,T) (for BC's).

"Moreover, if I understand well, it is in the hPsiThermo.C file that I call the REFPROP function?"

Yes, you can.
Regards,
Christian
Chris Lucas is offline   Reply With Quote

Old   November 3, 2010, 10:38
Default
  #13
Cyp
Senior Member
 
Cyprien
Join Date: Feb 2010
Location: Stanford University
Posts: 240
Rep Power: 9
Cyp is on a distinguished road
Yes I have phase change in my application. For the pressure equation, since I work in a porous medium, I developed an IMPES method (Implicit Pressure, Explicit Saturation). Pressures within each phase are related by a capillary pressure which depends on the saturation.

I was thinking about the creation of a second rho to rid of the perfect gas law. I tried to define a new base class which is a copy of basicPsiThermo (I called it testThermo1). In the equivalent of the makeBasicPsiThermo.H file, I could precised only two function : the Cthermo and the Mixture. So I could have only one rho variable. What do you think about that ?

In a second step, I defined a new class deriving from the previous class and MixtureType. (exactly what have been done for the hPsiThermo class, I called it testThermo2). Then I modify the equivalent of the hPsiThermos.C file in order to define only Cthermo and pureMixture.... but it doesn't compile.....

I get :

Code:
cyp@cyp-laptop:~/OpenFOAM/cyp-1.7.1/applications/lib$ wmake libso
Making dependency list for source file testThermo2/testThermos2.C
SOURCE=testThermo2/testThermos2.C  ;  g++ -m64 -Dlinux64 -DWM_DP -Wall -Wextra -Wno-unused-parameter  -Wold-style-cast -Wnon-virtual-dtor -O3  -DNoRepository  -ftemplate-depth-40 -I/opt/openfoam171/src/finiteVolume/lnInclude     -I/opt/openfoam171/src/thermophysicalModels/specie/lnInclude     -I/opt/openfoam171/src/thermophysicalModels/basic/lnInclude -IlnInclude -I. -I/opt/openfoam171/src/OpenFOAM/lnInclude -I/opt/openfoam171/src/OSspecific/POSIX/lnInclude   -fPIC -c $SOURCE -o Make/linux64GccDPOpt/testThermos2.o
testThermo2/testThermos2.C:49: error: type/value mismatch at argument 1  in template parameter list for ‘template<class MixtureType> class  Foam::testThermo2’
testThermo2/testThermos2.C:49: error:   expected a type, got ‘pureMixture’
testThermo2/testThermos2.C:49: error: invalid type in declaration before ‘;’ token
testThermo2/testThermos2.C:49: error: ‘testThermo2pureMixture’ is not a class or namespace
testThermo2/testThermos2.C:49: error: ‘testThermo2pureMixture’ is not a class or namespace
testThermo2/testThermos2.C:49: warning: unused variable ‘Foam::debug’
testThermo2/testThermos2.C: In constructor ‘Foam::testThermo1::addfvMeshConstructorToTable<testThermo1Type>::addfvMeshConstructorToTable(const Foam::word&) [with testThermo1Type = int]’:
testThermo2/testThermos2.C:56: error: ‘typeName’ is not a member of ‘int’
In file included from lnInclude/makeTestThermo1.H:35,
                 from testThermo2/testThermos2.C:28:
lnInclude/testThermo1.H: In static member function ‘static Foam::autoPtr<Foam::testThermo1> Foam::testThermo1::addfvMeshConstructorToTable<testThermo1Type>::New(const Foam::fvMesh&) [with testThermo1Type = int]’:
lnInclude/testThermo1.H:78:   instantiated from ‘Foam::testThermo1::addfvMeshConstructorToTable<testThermo1Type>::addfvMeshConstructorToTable(const Foam::word&) [with testThermo1Type = int]’
testThermo2/testThermos2.C:56:   instantiated from here
lnInclude/testThermo1.H:71: error: cannot convert ‘const Foam::fvMesh’ to ‘int’ in initialization
I do not understand why testThermo2pureMixture is written a single word...

I think it should come from my makeTestThermo.H file :
Quote:
#ifndef makeTestThermo1_H
#define makeTestThermo1_H

#include "testThermo1.H"
#include "addToRunTimeSelectionTable.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#define makeTestThermo1(Cthermo,Mixture) \
\
typedef Cthermo<Mixture> \
Cthermo##Mixture; \
\
defineTemplateTypeNameAndDebugWithName \
( \
Cthermo##Mixture, \
#Cthermo \
"<"#Mixture">", \
0 \
); \
\
addToRunTimeSelectionTable \
( \
testThermo1, \
Cthermo##Mixture, \
fvMesh \
)


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#endif
I can't see my mistake...

Best Regards,
Cyp

Last edited by Cyp; November 3, 2010 at 11:36.
Cyp is offline   Reply With Quote

Old   November 16, 2010, 21:40
Default
  #14
Senior Member
 
Nima Samkhaniani
Join Date: Sep 2009
Location: Tehran, Iran
Posts: 1,166
Blog Entries: 1
Rep Power: 15
nimasam is on a distinguished road
how can i make dynamic file from a fortan code and use it in openFoam ?
nimasam is offline   Reply With Quote

Old   January 20, 2011, 04:49
Default Basicthermo.H
  #15
New Member
 
karthik
Join Date: Dec 2010
Location: munich
Posts: 16
Rep Power: 7
karthik1414 is on a distinguished road
hello everyone,
I am new to the OpenFoam software.
Can anyone tell me how I can make use of Basicthermo.H while using Interfoam as solver??
karthik1414 is offline   Reply With Quote

Old   January 20, 2011, 05:23
Default
  #16
Senior Member
 
Nima Samkhaniani
Join Date: Sep 2009
Location: Tehran, Iran
Posts: 1,166
Blog Entries: 1
Rep Power: 15
nimasam is on a distinguished road
may be its good to open new thread for ur issue, however you should look at user guide chapter 3 but a brief description can be like that:
in interfoam.C u should include ur new class:
# include "Basicthermo.H"
then in case directory use "wmake" to compile it.
nimasam is offline   Reply With Quote

Old   January 24, 2011, 09:43
Default
  #17
New Member
 
karthik
Join Date: Dec 2010
Location: munich
Posts: 16
Rep Power: 7
karthik1414 is on a distinguished road
hey nimasam,
Thank you for your reply, i tried that.
Any idea how to vary density with pressure change??
karthik1414 is offline   Reply With Quote

Old   January 24, 2011, 09:56
Default
  #18
Senior Member
 
Nima Samkhaniani
Join Date: Sep 2009
Location: Tehran, Iran
Posts: 1,166
Blog Entries: 1
Rep Power: 15
nimasam is on a distinguished road
look in compressible solver for example "compressibleInterFoam" and compare it with interFoam
nimasam is offline   Reply With Quote

Old   March 7, 2011, 15:15
Default
  #19
New Member
 
Evren
Join Date: Mar 2010
Posts: 20
Rep Power: 7
pbe_cfd is on a distinguished road
Hi,

I have couple of fortran.90 subroutines which I would like to call from OpenFOAM1.6. What I know direct call of fortran subroutine from c++ is not good. Thus, one has to have C interface to routines and some more technical difficulties as, fortran stores rowwise C stores columnwise Whatever, these are general problems of C++ and Fortran mixing. My questions are more OpenFOAM specific:

1) What should be the compiled fortran source codes, dynamic library, shared objects (.so), objects (.o) ... ? If more than one is good, which one performs better?

2) How should be this "executables" (compiled files) linked to OpenFOAM, let's say scalarFoam? How should the <file> and <options> in Make folder be modified?

Please no suggestion like "why don't you rewrite your subroutines in C++" ;-)

All the best,

@ Cyp: You see, I am chasing you ;-)
pbe_cfd is offline   Reply With Quote

Old   March 9, 2011, 06:15
Default
  #20
Cyp
Senior Member
 
Cyprien
Join Date: Feb 2010
Location: Stanford University
Posts: 240
Rep Power: 9
Cyp is on a distinguished road
Hi!

To link a library with OpenFOAM you can follow this 2 steps :

- First ,you create a *.so file using "wmake libso"
- Then you link this *.so file to OpenFOAM by including the following command in your controlDict file :
Code:
libs
(
    "*.so "
);
Best,
Cyp
Cyp 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
Call Fortran program from OpenFoam tsjb00 OpenFOAM Running, Solving & CFD 6 April 22, 2013 15:26
2D Mesh Generation Tutorial for GMSH aeroslacker Open Source Meshers: Gmsh, Netgen, CGNS, ... 12 January 19, 2012 04:52
Native OpenFOAM interface in Pointwise Chris Sideroff Main CFD Forum 0 January 16, 2009 13:37
Convective Heat Transfer - Heat Exchanger Mark CFX 6 November 15, 2004 16:55
Replace periodic by inlet-outlet pair lego CFX 3 November 5, 2002 21:09


All times are GMT -4. The time now is 05:37.