CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (http://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   interFoam: contact angle (http://www.cfd-online.com/Forums/openfoam-programming-development/68768-interfoam-contact-angle.html)

ala October 1, 2009 06:28

interFoam: contact angle
 
Dear community,

I have a question concerning interFOAM: Where/How is the contactAngle during the simulation calculated? Respectively where/how gets the contactAngle respected during the simulation?Unfortunately I cannot figure out the corresponding code part.

Any hints/help would be appreciated!
Thanks in advance,
Ala

kathrin_kissling October 2, 2009 05:26

Take a look into ../src/transportProperties/interfaceProperties/interfaceProperties.C / interfaceProperties.H the member function is called correctContactAngle. The definite calculation of the contact angle is in src/transportProperties/interfaceProperties/calculateContactAngle and in the subfolders. Just take a look. If you have any questions, feel free to ask Kathi

ala October 5, 2009 11:04

Hi Kathi,

thank you for your reply. correctContactAngle is for correcting dynamicContactAngle, after the gamma equation being solved. But how gets the constantContactAngle during the simulation respected?

Greetings, Ala

kathrin_kissling October 6, 2009 07:50

Hi ala,
it is respected in UEqn.H (in case you enable the momentum corrector) and in the pEqn.
There is one term

fvc::interpolate(interface.sigmaK())*fvc::snGrad(g amma))

which is described by Brackbill 1992.

sigmaK is a member function of interfaceProperties class. And calls calculateContact angle.

Hope this helps

Kathrin

ala September 27, 2010 11:24

Thank you : ) it had helps,

greetings, ala

duongquaphim September 28, 2010 11:41

Dear all,

I would like to know the differences between constant contact angle and dynamic contact angle in OF. I did some simulation with steady moving bubble in a straight channel with prescribed constant contact angle (30 degree) at the wall. What I expected is that the contact angle close to the wall at the front and the rear bubble is similar. But what I got from OF is not (you can see in the figure).

http://i189.photobucket.com/albums/z84/haduong83/30.jpg

I did check the code of constantAlphaContactAngle and found that there is no function to calculate new contact angle based on alpha_advance, alpha_reduce and u_theta as I saw in dynamicAlphaContactAngle. However, I found this line in interfaceProperties.C:

scalarField theta =
convertToRad*acap.theta(U_.boundaryField()[patchi], nHatp)

which calculate contact angle at the wall. Here, we can see the appearance of U_.boundaryField()[patchi]. So U must be taken into account in calculating contact angle in OF even we use constantAlphaContactAngle option. That might explain my observation.

Do you have any other ideas about that?

Regards,

Duong

mohsenkh599 October 12, 2010 08:29

compiling error
 
Hi every body,
I have a problem when I want to compile a model.
the error is:

************************************************** *********************
************************************************** *********************

viscoelasticLaws/mohsen/mohsen.C:159: error: no match for ‘operator==’ in ‘Foam::operator+(const Foam::tmp<Foam::fvMatrix<Type> >&, const Foam::tmp<Foam::fvMatrix<Type> >&) [with Type = Foam::SymmTensor<double>](((const Foam::tmp<Foam::fvMatrix<Foam::SymmTensor<double> > >&)((const Foam::tmp<Foam::fvMatrix<Foam::SymmTensor<double> > >*)(& Foam::fvm::div(const Foam::surfaceScalarField&, Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>&) [with Type = Foam::SymmTensor<double>](((Foam::GeometricField<Foam::SymmTensor<double>, Foam::fvPatchField, Foam::volMesh>&)(&((Foam::mohsen*)this)->Foam::mohsen::a_))))))) == Foam::operator+(const Foam::tmp<Foam::GeometricField<TypeR, PatchField, GeoMesh> >&, const Foam::tmp<Foam::GeometricField<Type1, PatchField, GeoMesh> >&) [with Type1 = Foam::Tensor<double>, Type2 = Foam::SymmTensor<double>, PatchField = Foam::fvPatchField, GeoMesh = Foam::volMesh](((const Foam::tmp<Foam::GeometricField<Foam::SymmTensor<do uble>, Foam::fvPatchField, Foam::volMesh> >&)((const Foam::tmp<Foam::GeometricField<Foam::SymmTensor<do uble>, Foam::fvPatchField, Foam::volMesh> >*)(& Foam::operator*(const Foam::dimensioned<double>&, const Foam::tmp<Foam::GeometricField<Type, PatchField, GeoMesh> >&) [with Type = Foam::SymmTensor<double>, PatchField = Foam::fvPatchField, GeoMesh = Foam::volMesh](((const Foam::tmp<Foam::GeometricField<Foam::SymmTensor<do uble>, Foam::fvPatchField, Foam::volMesh> >&)((const Foam::tmp<Foam::GeometricField<Foam::SymmTensor<do uble>, Foam::fvPatchField, Foam::volMesh> >*)(& Foam::operator&(const Foam::tmp<Foam::GeometricField<TypeR, PatchField, GeoMesh> >&, const Foam::tmp<Foam::GeometricField<Type1, PatchField, GeoMesh> >&) [with Type1 = Foam::SymmTensor<double>, Type2 = Foam::SymmTensor<double>, PatchField = Foam::fvPatchField, GeoMesh = Foam::volMesh](((const Foam::tmp<Foam::GeometricField<Foam::SymmTensor<do uble>, Foam::fvPatchField, Foam::volMesh> >&)((const Foam::tmp<Foam::GeometricField<Foam::SymmTensor<do uble>, Foam::fvPatchField, Foam::volMesh> >*)(& Foam::operator-(const Foam::dimensioned<Type>&, const Foam::tmp<Foam::GeometricField<Type1, PatchField, GeoMesh> >&) [with Form = Foam::SymmTensor<double>, Type = Foam::SymmTensor<double>, PatchField = Foam::fvPatchField, GeoMesh = Foam::volMesh](((const Foam::tmp<Foam::GeometricField<Foam::SymmTensor<do uble>, Foam::fvPatchField, Foam::volMesh> >&)((const Foam::tmp<Foam::GeometricField<Foam::SymmTensor<do uble>, Foam::fvPatchField, Foam::volMesh> >*)(& Foam::operator*(const Foam::scalar&, const Foam::GeometricField<Type, PatchField, GeoMesh>&) [with Type = Foam::SymmTensor<double>, PatchField = Foam::fvPatchField, GeoMesh = Foam::volMesh](((const Foam::GeometricField<Foam::SymmTensor<double>, Foam::fvPatchField, Foam::volMesh>&)((const Foam::GeometricField<Foam::SymmTensor<double>, Foam::fvPatchField, Foam::volMesh>*)(&((Foam::mohsen*)this)->Foam::mohsen::a_))))))))))))))))))))’

************************************************** ************
************************************************** ************

and the model is :

// Velocity gradient tensor
volTensorField L = fvc::grad( U() );

// Convected derivate term
// volTensorField C = tau_ & L;

// vorticity term
volTensorField Vor = L - L.T();

// Twice the rate of deformation tensor
volSymmTensorField twoD = twoSymm( L );


// fc
scalar fc = 0.5; //1 - 27 * det( a_ );

//double dot term

// volSymmTensorField mk = -( (2 / 35) * ( 1 - fc ) * twoD - (2 / 7) * ( a_ & twoD + twoD & a_ + (1 / 2) * tr( twoD & a_ ) * I_ ) - fc * tr( twoD & a_ ) * a_ );

// Stress transport equation
tmp<fvSymmTensorMatrix> aEqn
(
fvm::ddt(a_)
+ fvm::div(phi(), a_)
==
( 1 / 2 ) * ( Vor & a_ - a_ & Vor ) + ( 1 / 2 ) * keisi_ * ( twoD & a_ + a_ & twoD + 2 * ( ( (2 / 35) * ( 1 - fc ) * twoD - (2 / 7) * ( 1 - fc ) * ( a_ & twoD + twoD & a_ + (1 / 2) * tr( twoD & a_ ) * I_ ) - fc * tr( twoD & a_ ) * a_ ) ) + 2 * ci_ * twoD & ( I_ - 3 * a_ )
//+ twoSymm( C )
//- zeta_ / 2 * ( (tau_ & twoD) + (twoD & tau_) )
//- fvm::Sp( epsilon_ / etaP_ * tr(tau_) + 1/lambda_, tau_ )
);

aEqn().relax();
solve(aEqn);
// Viscoelastic stress
tau_ = (etaP_ / keisi_) * (a_ - I_);
}

************************************************** ****************
************************************************** ****************

where "a" , "tau" and "I" are symmetric tensors and "Vor" is asymmetric tensor. and the others are constant parameters.

best

kathrin_kissling October 12, 2010 08:39

Hi
likely the term
( 1 / 2 ) * ( Vor & a_ - a_ & Vor ) + ( 1 / 2 ) * keisi_ * ( twoD & a_ + a_ & twoD + 2 * ( ( (2 / 35) * ( 1 - fc ) * twoD - (2 / 7) * ( 1 - fc ) * ( a_ & twoD + twoD & a_ + (1 / 2) * tr( twoD & a_ ) * I_ ) - fc * tr( twoD & a_ ) * a_ ) ) + 2 * ci_ * twoD & ( I_ - 3 * a_ )

doesnt match the right side of your equation system. Obviously it should be of type of type symTensorMarix.

Can you check whether this is ok? Just put an Info staement including the term and look what field is printed.

I'm not really sure whether you are in the correct thread since this one focusses on contact angles while your problem seems to be viscoelastic...

Best
Kathrin

kathrin_kissling October 12, 2010 08:43

I just recognized now...

please do not crosspost!
People don't like their threads being spamed. Please keep to the forum netiquette otherwise it is likely you don't get answers at all!

Best Kathrin

herbert October 27, 2010 07:23

Quote:

Originally Posted by duongquaphim (Post 276956)
However, I found this line in interfaceProperties.C:

scalarField theta =
convertToRad*acap.theta(U_.boundaryField()[patchi], nHatp)

which calculate contact angle at the wall. Here, we can see the appearance of U_.boundaryField()[patchi]. So U must be taken into account in calculating contact angle in OF even we use constantAlphaContactAngle option. That might explain my observation.

Hi Duong,

I don't think so. U is not used when you use constantAlphaContactAngle. It has only to be given to the function theta() to match both cases: constant and dynamic contactAngle-BC.

There must be some other effects, that cause your bubble not being symmetric. Please keep in mind, that the contact angle is only corrected in the first cell layer. Of course inertia or some other effects influence the shape of the interface as well, and these effects are not symmetric.

Regards,
Stefan

idefix February 26, 2014 04:18

Hello,

I know the thread is older but I still didn´t get some points, perhaps someone can help me:

Quote:

Hi ala,
it is respected in UEqn.H (in case you enable the momentum corrector) and in the pEqn.
There is one term

fvc::interpolate(interface.sigmaK())*fvc::snGrad(g amma))

which is described by Brackbill 1992.

sigmaK is a member function of interfaceProperties class. And calls calculateContact angle.

Hope this helps

Kathrin

Does anyone know a good book/ paper where the constantAlphaContactAngle is explained?
I am not familiar with the code and I can´t understand what is going on by looking at the code.
I just want to know which equation is calculated when I use the boundary condition constantAlphaContactAngle .

Thanks a lot


All times are GMT -4. The time now is 11:25.