Adding the force from IBM(immersed boundary method) to assign the vel…

…ocity to zero forcibly.
iurnus · Feb 15, 2015 · c696f1e · c696f1e
1 parent d133d77
commit c696f1e
Show file tree

Hide file tree

Showing 2 changed files with 58 additions and 0 deletions.
diff --git a/lammpsFoam/UEqns.H b/lammpsFoam/UEqns.H
@@ -28,5 +28,18 @@ betaPhib = betaf*phib;
  + fvc::average(beta)*gradP.flowDirection()*gradP.value() // Adding pressure gradient
  );
 
+ if (addIBMForce)
+ {
+ UbEqn -= fvm::Sp(-ibmIndicatorPtr()/ibmRelaxTime, Ub);
+ }
+
+ // // Set the velocity directly
+ // if (addIBMForce)
+ // {
+ // vectorField zerosVelField(ibmIndicatorList.size(), Foam::vector::zero);
+ // UbEqn.setValues(ibmIndicatorList, zerosVelField);
+ // }
+
  UbEqn.relax();
+
 }
diff --git a/lammpsFoam/createFields.H b/lammpsFoam/createFields.H
@@ -254,3 +254,48 @@
 
 # include "createTurbulence.H"
 
+ Switch addIBMForce(transportProperties.lookupOrDefault<Switch>("addIBMForce", false));
+
+ autoPtr<volScalarField> ibmIndicatorPtr;
+
+ dimensionedScalar ibmRelaxTime
+ (
+ dimensioned<scalar>::lookupOrDefault
+ (
+ "ibmRelaxTime",
+ transportProperties,
+ runTime.deltaT().value() * scalar(3),
+ dimTime
+ )
+ );
+
+ DynamicList<label> ibmIndicatorList(0);
+
+ if (addIBMForce)
+ {
+ ibmIndicatorPtr.reset
+ (
+ new volScalarField
+ (
+ IOobject
+ (
+ "ibmIndicator",
+ runTime.timeName(),
+ mesh,
+ IOobject::MUST_READ,
+ IOobject::AUTO_WRITE
+ ),
+ mesh
+ )
+ );
+
+ volScalarField& ibmIndicatorField = ibmIndicatorPtr();
+
+ forAll(ibmIndicatorField, cellI)
+ {
+ if (ibmIndicatorField[cellI] > 0)
+ {
+ ibmIndicatorList.append(cellI);
+ }
+ }
+ }