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

Solver for FSI strongly coupled

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

Like Tree4Likes

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   February 17, 2011, 05:25
Post Solver for FSI strongly coupled
  #1
New Member
 
charbel's Avatar
 
Charbel Habchi
Join Date: Jan 2011
Location: Lebanon
Posts: 27
Rep Power: 15
charbel is on a distinguished road
Hi,
I am new user of OpenFOAM (OF-1.6) and I want to study a flow past a flexible, elastic, tab fixed on a channel wall. The first simulations I made were for non-flexible tab using pisoFoam and everything worked fine. Now I want to study the FSI between the flow and the flexible tab.

My questions are below:
1) Can I use icoFsiFoam for my case?
2) Is icoFsiFoam adapted for weak or strong dynamical coupling?
3) I tried to compile icoFsiFoam on OF-1.6 but several erros have occurred. Do I need to have the older version OF-1.4 or there is some solution to compile icoFsiFoam on OF-1.6?
4) Can icoFsiFoam be used for turbulent flows?

Thank you very much.
Charbel.
hua1015 and Komeil like this.
charbel is offline   Reply With Quote

Old   February 23, 2011, 10:43
Default icoFsiFoam with OpenFOAM-1.6-ext
  #2
New Member
 
charbel's Avatar
 
Charbel Habchi
Join Date: Jan 2011
Location: Lebanon
Posts: 27
Rep Power: 15
charbel is on a distinguished road
Hi,
Thanks Mathieu.

icoFsiFoam works with OpenFOAM-1.6-ext:
http://sourceforge.net/projects/openfoam-extend/develop

I have successfully compiled OpenFoam-1.6-ext and icoFsiFoam tutorial is working well.

Best regards,
Charbel.
charbel is offline   Reply With Quote

Old   April 13, 2011, 12:35
Exclamation FSI solver based on pimpleDyMFoam
  #3
New Member
 
charbel's Avatar
 
Charbel Habchi
Join Date: Jan 2011
Location: Lebanon
Posts: 27
Rep Power: 15
charbel is on a distinguished road
Hi,
enclosed please find a solver and a test case that I am trying to made for fluid/structure interaction basing on pimpleDyMFoam ...
I am not sure from the outer loop iterations I made (using Aitken or constant relaxation). Can any one give me some hints please!!
In fact when using the solver for the enclosed test case, the mesh blows after some time steps I do not know why!! The test case is the same as flappingConsole with some modifications.
Can you test try this solver for some other cases and tell me if it works?!!
Best regards,
Charbel.
Attached Files
File Type: gz vorflexFoam.tar.gz (8.1 KB, 139 views)
File Type: gz testVorflexFoam.tar.gz (63.0 KB, 137 views)
charbel is offline   Reply With Quote

Old   April 28, 2011, 05:46
Default
  #4
Member
 
Timo K.
Join Date: Feb 2010
Location: University of Stuttgart
Posts: 66
Rep Power: 16
timo_IHS is on a distinguished road
Hello Charbel,

your approach looks really nice.
something is wrong with the connection between solid and fluid, even for weak coupling! you should concentrate on that! I would first test it with solver of icoFsiFoam...

Best regards!
timo_IHS is offline   Reply With Quote

Old   May 1, 2011, 10:25
Default
  #5
New Member
 
charbel's Avatar
 
Charbel Habchi
Join Date: Jan 2011
Location: Lebanon
Posts: 27
Rep Power: 15
charbel is on a distinguished road
Hello Timo,

Thank you for your post.

For weak coupling I am using the same connection between Fluid & Solid as in icoFsiFoam. I do not understand where you see the error. Can you please precise?

Best regards,
Charbel.
charbel is offline   Reply With Quote

Old   May 2, 2011, 12:28
Default
  #6
Member
 
Timo K.
Join Date: Feb 2010
Location: University of Stuttgart
Posts: 66
Rep Power: 16
timo_IHS is on a distinguished road
When I start your solver with weakCoupling, the mesh blows.
With you it doesn't!?
timo_IHS is offline   Reply With Quote

Old   May 3, 2011, 06:26
Default
  #7
Member
 
Timo K.
Join Date: Feb 2010
Location: University of Stuttgart
Posts: 66
Rep Power: 16
timo_IHS is on a distinguished road
Okay I found the problem (for weak coupling):
you have to initialize the UsolidPrev like this:

volVectorField UsolidPrev
(
IOobject
(
"UsolidPrev",
runTime.timeName(),
stressMesh,
IOobject::NO_READ
),
Usolid
);
UsolidPrev = Usolid.oldTime();
timo_IHS is offline   Reply With Quote

Old   May 3, 2011, 09:14
Default
  #8
New Member
 
charbel's Avatar
 
Charbel Habchi
Join Date: Jan 2011
Location: Lebanon
Posts: 27
Rep Power: 15
charbel is on a distinguished road
Hello Timo,
Indeed I have not pay attention to this error, thank you a lot.
Please correct me if I am wrong: I can do the same thing if I modify the weakCoupling.H to:

{

# include "solveFluid.H"
# include "setPressure.H"
# include "solveSolid.H"
{
// Setting mesh motion
pointVectorField solidPointsDispl = cpi.interpolate(Usolid - Usolid.oldTime());

vectorField newPoints = stressMesh.points() + solidPointsDispl.internalField();

stressMesh.movePoints(newPoints);

interpolatorSolidFluid.movePoints();

vectorField fluidPatchPointsDispl =
interpolatorSolidFluid.pointInterpolate
(
solidPointsDispl.boundaryField()[solidPatchID].
patchInternalField()
);

motionUFluidPatch ==
tppi.pointToPointInterpolate
(
fluidPatchPointsDispl/runTime.deltaT().value()
);

mesh.update();

# include "volContinuity.H"

Info << "Motion magnitude: mean = "
<< average(mag(Usolid.boundaryField()[solidPatchID]))
<< " max = "
<< max(mag(Usolid.boundaryField()[solidPatchID])) << endl;
}

}

Best regards,
Charbel.
charbel is offline   Reply With Quote

Old   May 3, 2011, 10:03
Default
  #9
Member
 
Timo K.
Join Date: Feb 2010
Location: University of Stuttgart
Posts: 66
Rep Power: 16
timo_IHS is on a distinguished road
You're right, you can also use the original implementation.
But you can't replace
Usolid.oldTime()
by
volVectorField UsolidPrevious

you have to use the IOobject, which is related to "stressMesh".


BTW: what's your application for FSI?
timo_IHS is offline   Reply With Quote

Old   May 3, 2011, 11:11
Default
  #10
New Member
 
charbel's Avatar
 
Charbel Habchi
Join Date: Jan 2011
Location: Lebanon
Posts: 27
Rep Power: 15
charbel is on a distinguished road
must I also apply this when I define
volVectorField UsolidPrevious = 0*Usolid; //which is used in setMotion.H
volVectorField DeltaU1 = 0*Usolid;
volVectorField DeltaU2 = 0*Usolid;
volVectorField rr1 = 0*Usolid;
for the strong coupling?
In fact I want to apply the FSI for mixing enhancement.
charbel is offline   Reply With Quote

Old   May 6, 2011, 11:40
Exclamation
  #11
New Member
 
charbel's Avatar
 
Charbel Habchi
Join Date: Jan 2011
Location: Lebanon
Posts: 27
Rep Power: 15
charbel is on a distinguished road
I modified weakCoupling.H to:

Code:
{
# include "solveFluid.H"
# include "setPressure.H"
# include "solveSolid.H"
{
// Setting mesh motion
    pointVectorField solidPointsDispl = cpi.interpolate(Usolid - Usolid.oldTime());
    vectorField newPoints =
        stressMesh.points() + solidPointsDispl.internalField();
    stressMesh.movePoints(newPoints);
    interpolatorSolidFluid.movePoints();
     vectorField fluidPatchPointsDispl =
        interpolatorSolidFluid.pointInterpolate
        (
            solidPointsDispl.boundaryField()[solidPatchID].
            patchInternalField()
        );
    motionUFluidPatch ==
        tppi.pointToPointInterpolate
        (
            fluidPatchPointsDispl/runTime.deltaT().value()
        );

    mesh.update();
#   include "volContinuity.H"
    Info << "Motion magnitude: mean = " 
        << average(mag(Usolid.boundaryField()[solidPatchID]))
        << " max = " 
        << max(mag(Usolid.boundaryField()[solidPatchID])) << endl;
}
}
is it good now?
Best regards,
Charbel.
charbel is offline   Reply With Quote

Old   July 7, 2011, 08:54
Default
  #12
New Member
 
charbel's Avatar
 
Charbel Habchi
Join Date: Jan 2011
Location: Lebanon
Posts: 27
Rep Power: 15
charbel is on a distinguished road
Hi everyone,

You can find here a new solver that I am trying to make for fluid-structure interaction.
I am working on version 1.6-ext.

This solver included:
1/ pimpleDyMFoam in which we have ALE and turbulence models can be used
2/ outer fixed point adaptive relaxation by using Aitken's method for strong coupling
3/ solve transport equations for heat and mass transfer
4/ solid solver for large deformattion (2nd Piola-Krichhoff stress tensor...)
5/ computing traction in setPressure.H but it is not taken into account in tractionDisplacementFvPatchVectorField.C can you give some hint please

Can you try this solver on the attached testCase and give me remarks because I am sure there will be some mistakes especially for:
1/ Aitken's method
2/ solid solver and its part in tractionDisplacementFvPatchVectorField.C
3/ setting traction in tractionDisplacementFvPatchVectorField.C


Thank you in advance.
Charbel.
Attached Files
File Type: gz solverAndTestCaseFSI.tar.gz (75.4 KB, 203 views)
charbel is offline   Reply With Quote

Old   July 7, 2011, 12:02
Default
  #13
Senior Member
 
Elvis
Join Date: Mar 2009
Location: Sindelfingen, Germany
Posts: 620
Blog Entries: 6
Rep Power: 24
elvis will become famous soon enough
Hello,

on the last OpenFoam workshop at Penn State FSI was one topic of the trainings.
http://www.openfoamworkshop.org/6th_...m/training.htm
http://www.openfoamworkshop.org/6th_...ven_slides.pdf http://www.openfoamworkshop.org/6th_...erTraining.tgz

*this might be best *http://www.openfoamworkshop.org/6th_...ell_slides.pdf http://www.openfoamworkshop.org/6th_...mpbell_fsi.tgz

Numerical simulation of sailing boats dynamics FSI and shape optimization RANS and LES of Turbulent Combustion based on Flamelet Generated Manifolds
Matteo Lombardi (École Polytechnique Fédérale de Lausanne)

more FSI stuff see other presentations
etc.

http://www.openfoamworkshop.org/6th_...rdi_slides.pdf
elvis is offline   Reply With Quote

Old   July 11, 2011, 06:57
Default
  #14
New Member
 
charbel's Avatar
 
Charbel Habchi
Join Date: Jan 2011
Location: Lebanon
Posts: 27
Rep Power: 15
charbel is on a distinguished road
Thank you Elvis, I ll see these presentations
__________________
--
Charbel
charbel is offline   Reply With Quote

Old   July 20, 2011, 12:21
Default fluid structure interaction problem
  #15
New Member
 
charbel's Avatar
 
Charbel Habchi
Join Date: Jan 2011
Location: Lebanon
Posts: 27
Rep Power: 15
charbel is on a distinguished road
I have small deformation values relative to weak coupling.
It seems that the problem comes from the first predicted interface displacement (see Kuttler and Wall 2008 page 64 section 3.1 before equation (22)).
Can you please give me a hint on how can I make it. Actually I am doing as follows (first 4 lines) but it is not working:
Code:
//******************************** Initialize for new time step ******************************//
# include "solveFluid.H"
# include "setPressure.H"
# include "solveSolid.H"
# include "setMotionOldTime.H"
//here solidPointsDispl = cpi.interpolate(Usolid - Usolid.oldTime())
//*******************************************************************************************//
int myLoopCount = 0;
do
{
  if (myLoopCount == 0)
   {
//********************* Computing Aitken's relaxation parameter for k = 0 ******************//
     minGammaAitken = min(GammaAitken.boundaryField()[solidPatchID]);
     minGammaAitken = std::sqrt(sqr(minGammaAitken));
     maxGammaAitken = max(GammaAitken.boundaryField()[solidPatchID]);
     maxGammaAitken = std::sqrt(sqr(maxGammaAitken));
     if (minGammaAitken > 1)
      {
        minGammaAitken = 0;
      }
     if (maxGammaAitken > 1)
      {
        maxGammaAitken = 0.9;
      }
     avgGammaAitken = (minGammaAitken + maxGammaAitken)/2;
     GammaAitken = min(avgGammaAitken, initialRelaxation) + 0*GammaAitken;
     avgGammaAitken = average(GammaAitken.boundaryField()[solidPatchID]);
     Info<< "*** myLoopCount = " << myLoopCount << " \n" << endl;
     Info<< "*** Aitken relaxation parameter (k=0) = " 
         << avgGammaAitken << " \n" << endl;
//*******************************************************************************************//
//******************************* Solving Fluid -> Solving Solid ****************************//
//#    include "setMotion.H"
#    include "solveFluid.H"
#    include "setPressure.H"
     UsolidPrevious = 1*Usolid;
#    include "solveSolid.H"
//*******************************************************************************************//
//*************************** Computing Aitken's residuales for k = 0 ***********************//
     rr2.boundaryField()[solidPatchID] = Usolid.boundaryField()[solidPatchID] - UsolidPrevious.boundaryField()[solidPatchID];
     rr2Size = rr2.boundaryField()[solidPatchID].size();
     outerError = std::sqrt(max(magSqr(rr2.boundaryField()[solidPatchID])/rr2Size));
     Info<< "*** Aitken loop error (k=0) = " << outerError << " \n" << endl;
     outerError = 1;
//*******************************************************************************************//
//********************************** Relaxe Aitken for k = 0 ********************************//
     Usolid.boundaryField()[solidPatchID] = UsolidPrevious.boundaryField()[solidPatchID]
                                            + GammaAitken.boundaryField()[solidPatchID]*rr2.boundaryField()[solidPatchID];
//*******************************************************************************************//
   }
  else
   {
   if (myLoopCount > 1)
    {
//********************* Computing Aitken's relaxation parameter for k > 0 ******************//
     GammaAitken.boundaryField()[solidPatchID] =
              (-1*GammaAitken.boundaryField()[solidPatchID]*
                (
                 (rr2.boundaryField()[solidPatchID] - rr1.boundaryField()[solidPatchID])
                 &rr1.boundaryField()[solidPatchID]
                )
              )
              /magSqr(rr2.boundaryField()[solidPatchID] - rr1.boundaryField()[solidPatchID]);
     minGammaAitken = min(GammaAitken.boundaryField()[solidPatchID]);
     minGammaAitken = std::sqrt(sqr(minGammaAitken));
     maxGammaAitken = max(GammaAitken.boundaryField()[solidPatchID]);
     maxGammaAitken = std::sqrt(sqr(maxGammaAitken));
     if (minGammaAitken > 1)
      {
        minGammaAitken = 0;
      }
     if (maxGammaAitken > 1)
      {
        maxGammaAitken = 0.9;
      }
     avgGammaAitken = (minGammaAitken + maxGammaAitken)/2;
     GammaAitken = avgGammaAitken + 0*GammaAitken;
     Info<< "*** myLoopCount = " << myLoopCount << " \n" << endl;
     Info<< "*** Aitken relaxation parameter = " 
         << avgGammaAitken << " \n" << endl;
//*******************************************************************************************//
//********************************** Relaxe Aitken for k > 0 ********************************//
     Usolid.boundaryField()[solidPatchID] = UsolidPrevious.boundaryField()[solidPatchID]
                                            + GammaAitken.boundaryField()[solidPatchID]*rr2.boundaryField()[solidPatchID];
//*******************************************************************************************//
    }
//********************** Solving Mesh -> Solving Fluid -> Solving Solid *********************//
#    include "setMotion.H"
//here solidPointsDispl = cpi.interpolate(Usolid - UsolidPrevious)
#    include "solveFluid.H"
#    include "setPressure.H"
     UsolidPrevious = 1*Usolid;
#    include "solveSolid.H"
//*******************************************************************************************//
//*************************** Computing Aitken's residuales for k > 0 ***********************//
     rr1.boundaryField()[solidPatchID] = rr2.boundaryField()[solidPatchID];
     rr2.boundaryField()[solidPatchID] = Usolid.boundaryField()[solidPatchID] - UsolidPrevious.boundaryField()[solidPatchID];
     rr2Size = rr2.boundaryField()[solidPatchID].size();
     outerError = std::sqrt(max(magSqr(rr2.boundaryField()[solidPatchID])/rr2Size));
     Info<< "*** Aitken loop error = " << outerError << "  \n" << endl;
     if (myLoopCount == 1)
      {
        outerError = 1;
      }
//*******************************************************************************************//
   }
  ++myLoopCount;
} while ((outerError > outerLoopConvergence) && (myLoopCount < nLoops));
Best regards.
__________________
--
Charbel
charbel is offline   Reply With Quote

Old   July 22, 2011, 09:50
Default traction displacement
  #16
New Member
 
charbel's Avatar
 
Charbel Habchi
Join Date: Jan 2011
Location: Lebanon
Posts: 27
Rep Power: 15
charbel is on a distinguished road
Hello,
I want to compute the traction field so it can be transferred, with the pressure, from fluid mesh to solid in the setPressure.H.
So I make as follows:

Code:
In createFields.H

    volSymmTensorField sigmaFluid
    (
        IOobject
        (
            "sigmaFluid",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
        2*nu*dev(symm(fvc::grad(U)))
    );

In setPressure.H

    sigmaFluid = 2*nu*dev(symm(fvc::grad(U)));
    
    vectorField faceNormal = (stressMesh.Sf().boundaryField()[solidPatchID])
                             /(stressMesh.magSf().boundaryField()[solidPatchID]);

    symmTensorField solidPatchTraction =
        interpolatorFluidSolid.faceInterpolate
        (
          (sigmaFluid.boundaryField()[fluidPatchID])
        );

    solidPatchTraction *= rhoFluid.value();

    tForce.traction() = solidPatchTraction & faceNormal;
However, it seems that the traction is not taken into account in the tractionDisplacementFvPatchVectorField.C
What most I modify in tractionDisplacementFvPatchVectorField.C .H so traction can be transferred to solid patch?


Regards.
__________________
--
Charbel
charbel is offline   Reply With Quote

Old   July 23, 2011, 00:31
Default
  #17
Member
 
Santiago
Join Date: Dec 2009
Posts: 85
Rep Power: 16
gascortado is on a distinguished road
Hi, I'm trying to test your solver. I untar the files and I placed the files in a folder. However, when I run Allrun in that folder, the command is not recognized. Do I have to place the files in a specifically maned folder? How do I run the case? Please let me know. THanks
gascortado is offline   Reply With Quote

Old   July 24, 2011, 12:29
Default
  #18
New Member
 
charbel's Avatar
 
Charbel Habchi
Join Date: Jan 2011
Location: Lebanon
Posts: 27
Rep Power: 15
charbel is on a distinguished road
Hi Santiago,
Did you first compile the solver without problems?
Once you compile it you go in the testCase directory and then ./Allrun
__________________
--
Charbel
charbel is offline   Reply With Quote

Old   July 24, 2011, 21:33
Default
  #19
Member
 
Santiago
Join Date: Dec 2009
Posts: 85
Rep Power: 16
gascortado is on a distinguished road
I execute "wmake libso" in "testVorflexFoam" and I get the following error:

testVorflexFoam.C:41: fatal error: tetFemMatrices.H: No such file or directory
compilation terminated.
make: *** [Make/linuxGccDPOpt/testVorflexFoam.o] Error 1

This also happened to me while trying to compile "icoFsiFoam". I tried with both OF-1.5-dev and OF-1.6-ext with the same results. Any idea of what I need to do now?
gascortado is offline   Reply With Quote

Old   July 25, 2011, 05:13
Default
  #20
New Member
 
charbel's Avatar
 
Charbel Habchi
Join Date: Jan 2011
Location: Lebanon
Posts: 27
Rep Power: 15
charbel is on a distinguished road
To compile solvers or utilities you moste execute wmake. The same for icoFsiFoam. In fact, you use wmake libso when you want to compile a library.

Regards.
__________________
--
Charbel
charbel is offline   Reply With Quote

Reply

Tags
elastic, icofsifoam, strongly weakly coupling

Thread Tools Search this Thread
Search this Thread:

Advanced Search
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 Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
NASA Rotor 37 & Pressure-Based Coupled Solver Razvan FLUENT 6 June 11, 2015 04:58
coupled solver wont work in star ccm+ richie Siemens 5 November 4, 2008 05:51
coupled solver in FLUENT ztdep FLUENT 1 September 14, 2006 23:34
The seperated solver and coupled one summer FLUENT 0 July 5, 2006 16:06
can use coupled solver to calculate combustion cfdfans FLUENT 0 November 25, 2003 05:48


All times are GMT -4. The time now is 07:44.