[MikeWorth]

I'm trying to use a dynamicAlphaContactAngle boundary to 'pin' a contact point at the furthest that it reaches during advancement. As far as I understand dynamic contact angles this should work if the receding contact angle is zero, with a non-zero advancing contact angle.

I've tried setting this up to no avail, and having looked into the code it doesn't seem to do what I'd expect. The contact angle returned seems to be calcuated as: theta0_ + (thetaA_ - thetaR_)*tanh(uwall/uTheta_)

At low velocity (uwall<<uTheta) the limit is theta0, which makes perfect sense.
At high velocity in either direction (tanh(uwall/uTheta) = ±1)it doesn't limit to thetaA or thetaR as I'd expect. In fact thetaA and thetaR don't even seem to be independent parameters - only the difference between the two actually matters.

Does anyone know why this is calculated the way it is? As far as I can tell, to get independent control of the three different angles I'd need to either overwrite this definition, or create a new boundary condition that behaves as I'd expect...

[MikeWorth]

I've looked closer, and as far as I can tell the code for dynamicAlphaContactAngle doesn't behave at all as expected. At high velocity in either direction it doesn't tend towards the advancing or receding contact angle.

I don't think that there is a standard theory for exactly how to transition between advancing and receding angles (there doesn't seem to be anything in my textbook). Reusing the tanh approach, I think that this calculation gives results much more in line with what I'd expect:

    u_ratio=uWall/uTheta
    t=tanh(u_ratio)
    if u_ratio >= 0:
        return (1-t)*theta0 + t*thetaA
    else:
        return (1+t)*theta0 + -t*thetaR

As shown against the original in this plot:


Unless anyone has an explanation I'll register this as a bug. Just having a bit of trouble with verifying my account on gitlab...


And here is the implementation of my modified calculation:

Foam::tmp<Foam::scalarField>
Foam::dynamicAlphaContactAngleFvPatchScalarField::theta
(
    const fvPatchVectorField& Up,
    const fvsPatchVectorField& nHat
) const
{




    if (uTheta_ < SMALL)
    {
        return tmp<scalarField>::New(size(), theta0_);
    }

    const vectorField nf(patch().nf());

    // Calculated the component of the velocity parallel to the wall
    vectorField Uwall(Up.patchInternalField() - Up);
    Uwall -= (nf & Uwall)*nf;

    // Find the direction of the interface parallel to the wall
    vectorField nWall(nHat - (nf & nHat)*nf);

    // Normalise nWall
    nWall /= (mag(nWall) + SMALL);

    // Calculate Uwall resolved normal to the interface parallel to
    // the interface
    scalarField uwall(nWall & Uwall);

    Info << "USING MIKE WORTHS MODIFIED DYNAMIC CONTACT ANGLE CALCULATION\n";

    //Calculation will likely give nonsense without warning if values do not satisfy thetaA>=theta0>=thetaR
    return theta0_ * (1 - tanh(uwall/uTheta_) * sign(uwall/uTheta_) )
           - thetaA_ * neg(uwall) * tanh(uwall/uTheta_)
           + thetaR_ * pos(uwall) * tanh(uwall/uTheta_);
    //Original calculation:
    //return theta0_ + (thetaA_ - thetaR_)*tanh(uwall/uTheta_);
}

[MikeWorth]

I've also dug up this thread from last year, which was started with a bug report both that thetaA and thetaR are "the wrong way round" and that the limiting behaviour of the funciton is not to either thetaA or thetaR.

The discussion charged off into a dispute about which way round to define the angle, and missed the other (in my mind far less subjective) issue that the function ought to tend to one of them at each extreme...

Bug reported: https://develop.openfoam.com/Develop.../-/issues/1746

[MS.slh]

Hi Mike, Are you aware if this has been fixed/modified in the most recent version of OpenFOAM?

[MikeWorth]

The latest release version is 20.06; and if you look in the source code you can see that it is unchanged compared to my bug report.

Version 20.12 is due out in December, hopefully the change will be in there. If you really need it right now then you could swap in the appropriate lines from the bug report and recompile it yourself.

[alextau]

Do you know what is an acceptable value for uTheta could be? I choose 0.01 for my test case (sliding drops on the inclined surface) and it seems that both dynamic angles are the same as the static one.

Thanks!

[MS.slh]

uTheta is a scaling factor with same dimensions as velocity (LT^-1). It affects the rate at which the contact angle changes when the capillary number changes. I usually run the model without contact angle definition and see what is the range of velocity near to the wall at the interface and based on that will choose an appropriate scaling factor, uTheta. Maybe the name "uTheta" is misleading because it's a velocity not an angle/angular velocity.

[alextau]

Thank you for the guidance.

I am new to linux.. To my understanding I have to compile after I change .C file. Is there any relatively easy tutorial for how to do it properly?

Alex.

[alimea]

Originally Posted by MikeWorth

I've also dug up this thread from last year, which was started with a bug report both that thetaA and thetaR are "the wrong way round" and that the limiting behaviour of the funciton is not to either thetaA or thetaR.

The discussion charged off into a dispute about which way round to define the angle, and missed the other (in my mind far less subjective) issue that the function ought to tend to one of them at each extreme...

Bug reported: https://develop.openfoam.com/Develop.../-/issues/1746

Hi Mike

Is this problem solved for OpenFOAM-v7?
Or I should replace the above line on the code and recompile it?
thanks

[MikeWorth]

Looking at the same line in v7 suggests that it also has the same bug - https://github.com/OpenFOAM/OpenFOAM...arField.C#L152

Making the same alteration and recompiling will probably work (I've not tried it, but it doesn't interact with much other code). After making the change run some very simple simulations and verify that the results are as expected. Probably also worth adding to the openfoam foundation bug-tracker to try and get it resolved in future versions.

[alimea]

Originally Posted by MikeWorth

Looking at the same line in v7 suggests that it also has the same bug - https://github.com/OpenFOAM/OpenFOAM...arField.C#L152

Making the same alteration and recompiling will probably work (I've not tried it, but it doesn't interact with much other code). After making the change run some very simple simulations and verify that the results are as expected. Probably also worth adding to the openfoam foundation bug-tracker to try and get it resolved in future versions.

Dear Mike

Thanks for your answer.
I replace that line with the code and recompiled it, however, I did not face a large difference in my results. Maybe in my case, that bug was not very important.

My new question is related to uTheta.

Could you please give me a more explanation about it? What is it?

I read your statement "I usually run the model without contact angle definition and see what is the range of velocity near to the wall at the interface and based on that will choose an appropriate scaling factor, uTheta." but I didn't get it. without contact angle definition means you run the case with zeroGradient BC?

I solved my case one time with uTheta = 2 and another time with uTheta = 0.01, and found that my case is very sensitive to the uTheta value.

Thank you

[mhaghq]

Hi Mike, @MikeWorth

Can you please post the steps to reproduce this graph? How are you selecting a uWall from your patch, what is your criteria there?
My understanding is that from all contact angles at the patch you take the max and min of those at each time step, Am I right?

I am trying to validate the dynamic contact angle BC and came across your reported explanation.

Thank you ��

[mhaghq]

Originally Posted by MikeWorth

I've looked closer, and as far as I can tell the code for dynamicAlphaContactAngle doesn't behave at all as expected. At high velocity in either direction it doesn't tend towards the advancing or receding contact angle.

I don't think that there is a standard theory for exactly how to transition between advancing and receding angles (there doesn't seem to be anything in my textbook). Reusing the tanh approach, I think that this calculation gives results much more in line with what I'd expect:

    u_ratio=uWall/uTheta
    t=tanh(u_ratio)
    if u_ratio >= 0:
        return (1-t)*theta0 + t*thetaA
    else:
        return (1+t)*theta0 + -t*thetaR

As shown against the original in this plot:


Unless anyone has an explanation I'll register this as a bug. Just having a bit of trouble with verifying my account on gitlab...


And here is the implementation of my modified calculation:

Foam::tmp<Foam::scalarField>
Foam::dynamicAlphaContactAngleFvPatchScalarField::theta
(
    const fvPatchVectorField& Up,
    const fvsPatchVectorField& nHat
) const
{




    if (uTheta_ < SMALL)
    {
        return tmp<scalarField>::New(size(), theta0_);
    }

    const vectorField nf(patch().nf());

    // Calculated the component of the velocity parallel to the wall
    vectorField Uwall(Up.patchInternalField() - Up);
    Uwall -= (nf & Uwall)*nf;

    // Find the direction of the interface parallel to the wall
    vectorField nWall(nHat - (nf & nHat)*nf);

    // Normalise nWall
    nWall /= (mag(nWall) + SMALL);

    // Calculate Uwall resolved normal to the interface parallel to
    // the interface
    scalarField uwall(nWall & Uwall);

    Info << "USING MIKE WORTHS MODIFIED DYNAMIC CONTACT ANGLE CALCULATION\n";

    //Calculation will likely give nonsense without warning if values do not satisfy thetaA>=theta0>=thetaR
    return theta0_ * (1 - tanh(uwall/uTheta_) * sign(uwall/uTheta_) )
           - thetaA_ * neg(uwall) * tanh(uwall/uTheta_)
           + thetaR_ * pos(uwall) * tanh(uwall/uTheta_);
    //Original calculation:
    //return theta0_ + (thetaA_ - thetaR_)*tanh(uwall/uTheta_);
}

Hi Mike, @MikeWorth

Can you please post the steps to reproduce this graph? How are you selecting a uWall from your patch, what is your criteria there?
My understanding is that from all contact angles at the patch you take the max and min of those at each time step, Am I right?

I am trying to validate the dynamic contact angle BC and came across your reported explanation.

Thank you