CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (http://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   anisotropic porous media in OpenFOAM (http://www.cfd-online.com/Forums/openfoam-programming-development/108441-anisotropic-porous-media-openfoam.html)

 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!

 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!

Regards
Hisham

 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 :
forAll(U,celli)
{
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,
Hisham

 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&)
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!
Code:

```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(); 00161 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(); 00171 00172    fvm.diag() += mesh.V()*sp.value(); 00173 00174    return tfvm; 00175 }```
Regards,
Hisham

 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