
[Sponsors] 
October 9, 2007, 18:00 
Using ODE in OpenFoam

#1 
Senior Member
Senthil Kabilan
Join Date: Mar 2009
Posts: 113
Rep Power: 8 
Hi,
I am new to OpenFoam and was trying to use the ODE solver. I took a look at the ODETest but was not able to understand how the equations are cast. can some one show me how I can solve a simple system like the following? C  INITIAL VALUES X=0.0D0 Y(1)=2.0D0 Y(2)=0.66D0 C  ENDPOINT OF INTEGRATION XEND=2.0D0 XDELTA=0.2 C  SYSTEM OF EQUATIONS F(1)=Y(2) F(2)=((1Y(1)**2)*Y(2)Y(1))/2.0 ODETest \**/ #include "argList.H" #include "IOmanip.H" #include "ODE.H" #include "ODESolver.H" #include "RK.H" using namespace Foam; // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // class testODE : public ODE { public: testODE() {} label nEqns() const { return 4; } void derivatives ( const scalar x, const scalarField& y, scalarField& dydx ) const { dydx[0] = y[1]; dydx[1] = y[0]  (1.0/x)*y[1]; dydx[2] = y[1]  (2.0/x)*y[2]; dydx[3] = y[2]  (3.0/x)*y[3]; } void jacobian ( const scalar x, const scalarField& y, scalarField& dfdx, Matrix<scalar>& dfdy ) const { dfdx[0] = 0.0; dfdx[1] = (1.0/sqr(x))*y[1]; dfdx[2] = (2.0/sqr(x))*y[2]; dfdx[3] = (3.0/sqr(x))*y[3]; dfdy[0][0] = 0.0; dfdy[0][1] = 1.0; dfdy[0][2] = 0.0; dfdy[0][3] = 0.0; dfdy[1][0] = 1.0; dfdy[1][1] = 1.0/x; dfdy[1][2] = 0.0; dfdy[1][3] = 0.0; dfdy[2][0] = 0.0; dfdy[2][1] = 1.0; dfdy[2][2] = 2.0/x; dfdy[2][3] = 0.0; dfdy[3][0] = 0.0; dfdy[3][1] = 0.0; dfdy[3][2] = 1.0; dfdy[3][3] = 3.0/x; } }; // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // Main program: int main(int argc, char *argv[]) { argList::validArgs.clear(); argList::validArgs.append("ODESolver"); argList args(argc, argv); word ODESolverName(args.args()[1]); testODE ode; autoPtr<odesolver> odeSolver = ODESolver::New(ODESolverName, ode); scalar xStart = 1.0; scalarField yStart(ode.nEqns()); yStart[0] = ::Foam::j0(xStart); yStart[1] = ::Foam::j1(xStart); yStart[2] = ::Foam::jn(2, xStart); yStart[3] = ::Foam::jn(3, xStart); scalarField dyStart(ode.nEqns()); ode.derivatives(xStart, yStart, dyStart); Info<< setw(10) << "eps" << setw(12) << "hEst"; Info<< setw(13) << "hDid" << setw(14) << "hNext" << endl; Info<< setprecision(6); for (label i=0; i<15; i++) { scalar eps = ::Foam::exp(scalar(i + 1)); scalar x = xStart; scalarField y = yStart; scalarField dydx = dyStart; scalarField yScale(ode.nEqns(), 1.0); scalar hEst = 0.6; scalar hDid, hNext; odeSolver>solve(ode, x, y, dydx, eps, yScale, hEst, hDid, hNext); Info<< scientific << setw(13) << eps; Info<< fixed << setw(11) << hEst; Info<< setw(13) << hDid << setw(13) << hNext << setw(13) << y[0] << setw(13) << y[1] << setw(13) << y[2] << setw(13) << y[3] << endl; } scalar x = xStart; scalar xEnd = x + 1.0; scalarField y = yStart; scalarField yEnd(ode.nEqns()); yEnd[0] = ::Foam::j0(xEnd); yEnd[1] = ::Foam::j1(xEnd); yEnd[2] = ::Foam::jn(2, xEnd); yEnd[3] = ::Foam::jn(3, xEnd); scalar hEst = 0.5; odeSolver>solve(ode, x, xEnd, y, 1e4, hEst); Info<< nl << "Analytical: y(2.0) = " << yEnd << endl; Info << "Numerical: y(2.0) = " << y << ", hEst = " << hEst << endl; Info << "\nEnd\n" << endl; return 0; } // ************************************************** *********************** // 

October 19, 2007, 14:23 
Hi,
I figured out how to so

#2 
Senior Member
Senthil Kabilan
Join Date: Mar 2009
Posts: 113
Rep Power: 8 
Hi,
I figured out how to solve ODE's using the above driver. 1) Change the dimension of the system. abel nEqns() const { return 4; } 2) Change the following equations accordingly { dydx[0] = y[1]; dydx[1] = y[0]  (1.0/x)*y[1]; dydx[2] = y[1]  (2.0/x)*y[2]; dydx[3] = y[2]  (3.0/x)*y[3]; } 3) Specifiy the " scalar xStart " , " scalar xEnd " and you can specify the initial values for yStart's if required. Feel free to contact if not clear! 

October 19, 2007, 14:56 
Good work! Does the solution c

#3 
Senior Member
Srinath Madhavan (a.k.a pUl)
Join Date: Mar 2009
Location: Edmonton, AB, Canada
Posts: 700
Rep Power: 12 
Good work! Does the solution compare well to the analytical result?


February 19, 2008, 18:01 
Hello there
Does anyone kno

#4 
New Member
Bernardo D Flores
Join Date: Mar 2009
Posts: 21
Rep Power: 8 
Hello there
Does anyone knows how to integrate the ODE on OpenFOAM, i.e., using a solver like icoFoam and a ODE system serving as input on icoFOAM inlets to be solved at each time step. Any pointers are greatly appreciated Regards 

September 2, 2015, 19:30 

#5 
New Member
methma Rajamuni
Join Date: Jul 2015
Posts: 8
Rep Power: 2 
Hi,
Can you please explain in more details the steps of solving an ODE in OpenFOAM? Thanks, Methma 

September 14, 2015, 02:54 

#7 
New Member
methma Rajamuni
Join Date: Jul 2015
Posts: 8
Rep Power: 2 
Hi Hassan ,
I went through the tutorial you mentioned. I changed the ODETest.C file so that it contains my equations. After that I run the wmake command to complete the C file. It was compiled by making dependency list for source file. But I don't know how to get the results of the ODE solver. Can you please help me? I am new to OpenFOAM. I really need to understand this. My final goal is to work on fluid structure interface (Study on vortex induced vibration of a sphere/cylinder on free stream). There I need to solve the structure motion ODE in each time step. Thank you. Methma 

September 14, 2015, 03:12 

#8 
New Member
methma Rajamuni
Join Date: Jul 2015
Posts: 8
Rep Power: 2 
Hi Hassan ,
I went through the tutorial you mentioned. I changed the TestODE.C file so that it contains my equations. After that I run the wmake command to complete the C file. It was compiled by making dependency list for source file. But I don't know how to get the results of the ODE solver. When I compile the C file does it run automatically or do I need to run it separately? Can you please help me? I am new to OpenFOAM. I really need to understand this. My final goal is to work on fluid structure interface (Study on vortex induced vibration of a sphere/cylinder on free stream). There I need to solve the structure motion ODE in each time step. Thank you. Methma 

September 14, 2015, 03:24 

#9 
Senior Member
Hassan Kassem
Join Date: May 2010
Location: UK
Posts: 138
Rep Power: 7 
Hello,
I am working on similar project but for 3D wings. The results of the ODE system are stored in y (in TestODE) which basically the initial conditions yStart. Best wishes, Hassan 

Yesterday, 01:17 

#10 
New Member
methma Rajamuni
Join Date: Jul 2015
Posts: 8
Rep Power: 2 
Hi Hassan,
I went through the ODETest files and figured out how it works. I have implemented the ODE solver for the solid motion inside the icoFoam solver. Now, I need to take the force data to solve the solid motion solver in each time step. I went through the force.C and force.H files, but still have no idea how to use them to take the force data. I know the force.C file was written to write the force data, but I don't know how to edit those files to get the force data inside the icoFoam solver. I am not quite familiar with the library files and how OpenFoam uses the library files when we run a solver. So can you please help me on this? Thanks. Methma 

Yesterday, 06:17 

#11  
Senior Member
Hassan Kassem
Join Date: May 2010
Location: UK
Posts: 138
Rep Power: 7 
Quote:
That is a great progress, well done. Regarding the forces class, it is not only for writing. You can use it to get the forces directly [LINK]. There is a very good example for that in OpenFOAM source in sixDoFRigidBody [LINK] library. Here is a snippet of this library; Code:
// define a dictionary first dictionary forcesDict; // add the required parameters forcesDict.add("type", forces::typeName); forcesDict.add("patches", wordList(1, ptPatch.name())); forcesDict.add("rhoInf", rhoInf_); forcesDict.add("rhoName", rhoName_); forcesDict.add("CofR", motion_.centreOfRotation()); // create force object forces f("forces", db(), forcesDict); // calculate the forces f.calcForcesMoment(); // finally you can get the forces and moment vector F = f.forceEff(); vector M = f.momentEff(); Hassan 

Yesterday, 18:38 

#12 
New Member
methma Rajamuni
Join Date: Jul 2015
Posts: 8
Rep Power: 2 
Thanks a lot Hassan, I will go through the links you mentioned. I hope this is very helpful.
Best, Methma 

Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
OpenFoam vs CFX5 mass balance in OpenFoam  tangd  OpenFOAM Running, Solving & CFD  33  May 23, 2010 16:36 
CheckMesh error using a tutorial from OpenFOAM 114 with openFOAM 13  martapajon  OpenFOAM Native Meshers: blockMesh  7  January 21, 2008 13:52 
OpenFOAM users in Munich OpenFOAM benutzer in M%c3%bcnchen  jaswi  OpenFOAM  0  August 3, 2007 13:11 
A new Howto on the OpenFOAM Wiki Compiling OpenFOAM under Unix  mbeaudoin  OpenFOAM Installation  2  April 28, 2006 08:54 
Is GUI there for OpenFoam  navaladi  OpenFOAM  3  March 28, 2006 09:08 