CFD Online Discussion Forums

CFD Online Discussion Forums (
-   OpenFOAM Programming & Development (
-   -   anisotropic porous media in OpenFOAM (

hawkeye321 October 23, 2012 18:47

anisotropic porous media in OpenFOAM
Does anyone have any experience in simulating flow in an-isotropic porous media in OF? I mean simulations in which the permeability in x direction is different from the permeability in y direction.

My main problem is that when I am adding the darcy soource term as
+ ((1/Perm)*(fvm::Sp(Visc,U)))

I do not know how to tell the OF that Perm will be different based on the direction!

comments are appreciated!

nimasam October 24, 2012 04:48

hi mahdi
is in this condition permeability ("Prem") a vectorField?

if yes, you can define a volVectorField for Prem, then
assign your directional variation with setFields in 0 folder

hawkeye321 October 24, 2012 09:37

Permeability Vector
Hi nimasam
Thanks for the comment.
But, if I define the permeability as volVectorField then, when I want to multiply it by velocity what kind of vector multiplication do I have to use; definitely not inner or cross product. Right?

Hisham October 24, 2012 09:52

Hi Mahdi

In the case of non-isotropic permeability, the permeability is a symmetric tensor. For anisotropy, you should assign zero off-diagonal values and be left with kxx, kyy and kzz.

Hope this helps!


hawkeye321 October 24, 2012 09:57

multiplying permeability tensor and velocity vector
Dear Hisham
Thanks for the comment. I know that true, but, as I mentioned in my previous post, I am not sure how to multiply permeability tensor and velocity vector in the momentum equation. What kind of multiplication operator do I have to use?

nimasam October 24, 2012 10:26

dear mahdi,
if you want to divide Ux/Kx , Uy/Ky , Uz/Kz
then i guess you may want to use some thing like this :
S[celli].x() = U[celli].x()/Kx[celli].x();

Hisham October 24, 2012 10:35

I think you first would have to define the tensor as k (kxx 0 0 kyy 0 kzz), then get the inverse of the tensor as inv(k). The product should be an inner product: + (inv(k) & (fvm::Sp(Visc,U))

I would appreciate it if you report what happened!

Best regards,

hawkeye321 October 24, 2012 10:55

Dear Hisham
I have define a tensor named PermTenInv which has diagonal terms of 1/kxx,1/kyy and 1/kzz. Now when i am doing explicit calculations in form of

+ ((PermTen) & (fvc::Sp(Visc,U)))

there is no error. However, for implicit calculations in form of

((PermTen) & (fvm::Sp(Visc,U)))

I get the following error
BR1IM003.C:84:40: error: no match for ‘operator&’ in ‘PermTen & Foam::fvm::Sp(const Foam::dimensionedScalar&, Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>&) [with Type = Foam::Vector<double>, Foam::dimensionedScalar = Foam::dimensioned<double>](((Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>&)(& U)))’
/opt/openfoam171/src/OpenFOAM/lnInclude/wordI.H:167:19: note: candidates are: Foam::word Foam::operator&(const Foam::word&, const Foam::word&)
/opt/openfoam171/src/OpenFOAM/lnInclude/dimensionSet.H:273:29: note: Foam::dimensionSet Foam::operator&(const Foam::dimensionSet&, const Foam::dimensionSet&)
/opt/openfoam171/src/finiteVolume/lnInclude/readPISOControls.H:8:10: warning: unused variable ‘momentumPredictor’
/opt/openfoam171/src/finiteVolume/lnInclude/readPISOControls.H:11:10: warning: unused variable ‘transonic’
/opt/openfoam171/src/finiteVolume/lnInclude/readPISOControls.H:14:9: warning: unused variable ‘nOuterCorr’
make: *** [Make/linux64GccDPOpt/BR1IM003.o] Error 1
Unfortunately, I have to discretise this source term implicitly. Any ideas how to fix this error?

Hisham October 24, 2012 11:23

Dear Mahdi,

It is really sad that this did not work directly. You can report this as a bug (assuming a similar feature should be available). If you are highly motivated for your problem you can try to add a new fvm::Sp function to fvmSup.C/.H that can take a "dimensionedVector & sp" argument. You then have to dig deep into OpenFOAM's number crunching mechanisms and figure out which elements of fvm.diag() should you multiply by which component. Here ends my ideas. I hope you can find an easier solution!

template<class Type>
00153 Foam::tmp<Foam::fvMatrix<Type> >
00154 Foam::fvm::Sp
00155 (
00156    const dimensionedScalar& sp,
00157    const GeometricField<Type, fvPatchField, volMesh>& vf
00158 )
00159 {
00160    const fvMesh& mesh = vf.mesh();
00162    tmp<fvMatrix<Type> > tfvm
00163    (
00164        new fvMatrix<Type>
00165        (
00166            vf,
00167            dimVol*sp.dimensions()*vf.dimensions()
00168        )
00169    );
00170    fvMatrix<Type>& fvm = tfvm();
00172    fvm.diag() += mesh.V()*sp.value();
00174    return tfvm;
00175 }


hawkeye321 October 24, 2012 12:44

THnaks Hisham. I am going to discuss it in a new post to be sure that this is a bug in OF.

AnjaMiehe September 16, 2013 09:05

Hello Mhadi,

is there any progress in the subject? I am coding something similar for solidification and encounter the same problem.

Regards, Anja April 7, 2016 03:34


All times are GMT -4. The time now is 09:42.