[mykkujinu2201]

Dear Forum,


I have a question of adding force term.

It is related to the previous thread (Weird result from 'phase field solver'. Need help.).

I am trying to use phase-field method as in the previous thread.

I think the previous post focus on the phase-field method, while this is about a more general one. Thus I post another thread.

I modified pimpleFoam by adding the part to solve Cahn-Hilliard equation and chemical potential.

As a result, I could get capillary force, .


I added this capillary force in to momentum equation as

fvVectorMatrix UEqn
(
    fvm::ddt(rho,U) + fvm::div(rhoPhi,U)
  + MRF.DDt(rho,U)
  - fvm::laplacian(muf,U)-(fvc::grad(U)&fvc::grad(mu))///fvc::div((mu*dev(T(fvc::grad(U)))))//(fvc::grad(U)&fvc::grad(mu))//(fvc::grad(U)&fvc::grad(muf))//fvc::div((mu*dev(T(fvc::grad(U)))))//(fvc::grad(U)&fvc::grad(mu))
== 
  fvc::Sp(pot,gradC)+ // f_st=pot*fvc::grad(C) or f_st +
fvOptions(rho, U)
);

The test case is simulating a static bubble.

As time goes, all the residuals seem converged, however, it suddenly diverges.

I tested smaller grids and smaller time steps.
However, they ends at the same or similar moment, .

I also know there is another method to apply force term by adding force flux in pEqn.H.

However, it does not help at all and simulation ends as above.

Changing time scheme (Euler, CN, backward) and other discretization schemes did not solve the issue at all.


Though I am not sure, I think the equation is unstable, and eventually, it diverges.

Also, one advisor recommends adding a stabilization term in the momentum equation.

However, I do not know how to derive a stabilization term.


My questions are

1. Is it due to the stability of the momentum equation as I thought?

2. If so, how can I solve this issue? (How can I increase stability?)

3. How can I derive the stabilization term for the source term in the momentum equation?

4. I know fvm::Sp is used to solve source terms in an implicit manner. However, I also found that fvc::Sp works.
As far as I know, fvc stands for explicit expression.
What is the difference between fvm::Sp and fvc::Sp. I also wonder what is the difference between fvc::Sp and source term without using fvc::Sp.


Thank you for reading this thread.

Best regards,
Junkyu Kim

[msaravia]

Just a hint... Did you tried under-relaxing the source term? Also you could try to calculated the gradient using reconstruction, it has worked for me in the past for problems with discontinuities. The gradient should be something like

fvc::reconstruct(fvc::snGrad(C) * mesh.magSf())

[mykkujinu2201]

Dear Martin,

Thank you for the advice.

I will try to calculate the gradient with reconstruction.

Is it kind of smoothing by interpolation from cell to face and face to cell?

About under-relaxing, is there a way to relax force term?

It would help me a lot if you provide me any material to relax the force term in OpenFOAM.

I only know relaxing equations and fields such as velocity, pressure, and scalar field in scalar transport equation but not force term.

Best,
Jun

[mykkujinu2201]

Dear Martin,

I tried the run with the reconstruction method but it diverges at the same moment.

Now, I will try the run with relaxing force term as below.

I added below just before fvVectorMatrix UEqn.

f_st.relax()

And I put relaxation value at fvSolution.

Is it correct?

Best,
Jun

[mAlletto]

Just a question regarding the static bubble: Should in this case the chemical source term be in equilibrium with the pressure gradient?

If this is the case the gradient of the pressure and of the chemical potential should be discretized in exaclty the same manner in order to avoid spurious currents

[mykkujinu2201]

Dear Alletto,

First of all, thank you for the reply.

When I just directly added capillary force term, they are not discretized in the same way.

Chemical potential is discretized following fvc::grad in fvScheme, while pressure is discretized following fvc::snGrad in fvScheme.

After I read Martin's comment above, I modified the formulation.
Both gradient terms are discretized in the same way now as far as I know.

By the way, is there a way to discretize each term other than using fvScheme in OpenFOAM?

Best,
Jun

[mAlletto]

The choice of the discretion schemes goes easiest over fvschemes

[msaravia]

Originally Posted by mykkujinu2201

Dear Martin,

I tried the run with the reconstruction method but it diverges at the same moment.

Now, I will try the run with relaxing force term as below.

I added below just before fvVectorMatrix UEqn.

f_st.relax()

And I put relaxation value at fvSolution.

Is it correct?

Best,
Jun

Hi Ju, I think it is correct. Please note I have no expertise in your problem, I just gave you an idea of things that worked for me in the past with magnetostatics solvers.