[adib_mohammad]

I have created a modified solver based on simpleFoam and called it twosimpleFoam. I added one scalar tranport equation to track number density of particles () in a laminar flow.
Below is the transport equation I am trying to solve.



I have created "NEqn.H" file similar and included it in twosimpleFoam.C. I have two versions of "NEqn.H". The only difference how I solve the source source term.
The first one ends up with "core dump", but second one converges. I don't understand which one is the correct way to solve this equation and how these two implementations are different.

First Implementation (DIVERGES!):

{
    fvScalarMatrix NEqn
    (
        fvm::ddt(N_agg)
        + fvm::div(phi, N_agg)
        - fvm::laplacian(diff_prefactor*diff_p, N_agg)
    ==
        -fvm::Sp(0.5 * beta * N_agg * Av * rho, N_agg))

    );

    NEqn.relax();
    fvOptions.constrain(NEqn);
    NEqn.solve();
    fvOptions.correct(N_agg);
}

Second Implementation (CONVERGES!):

{
    fvScalarMatrix NEqn
    (
        fvm::ddt(N_agg)
        + fvm::div(phi, N_agg)
        - fvm::laplacian(diff_prefactor*diff_p, N_agg)
    ==
        fvOptions(N_agg)

    );

    NEqn.relax();
    fvOptions.constrain(NEqn);
    solve(NEqn == -fvm::Sp(0.5 * beta * N_agg * Av * rho, N_agg));
    fvOptions.correct(N_agg);
}

[adib_mohammad]

Any suggestion?

[dlahaye]

I would like to understand this as well.

In the case where the simulations do converge, does the computed solution make physical sense? And if it does, do you obtain convergence by flipping the sign, i.e., solve == + fvm::Sp(...).

Same question in case where simulations diverge. Should you flip the sign? Or should you replace the function fvm::Sp() by fvm::SuSp()? See e.g. this post fvMatrix, fvOptions and SuSp: automatic implicit/explicit source-term treatment.

How good is converge in case you do converge?

More advanced, and possibly not worth the effort: can you run in a debugger, print the matrix to file and check diagonal dominance of the matrix?

Let us know.

[MamboJambo]

Hello,

Regarding the printing of the matrix to a file.



You can just use.diag(), lower() and uppper() and .source(). However, this does not include the boundary contributions, right? Is there a way to print the matrix with boundary conditions contributions?

[adib_mohammad]

The converged solution is not physical. I am looking into it to figure out how to get physical results from this. I tested fvm::SuSp, and I am fairly familiar with its difference from fvm::Sp.

The convergence is acheived only by placing the source term after NEqn.relax();

I have also tested the signs, but it did not affect the results.

Originally Posted by dlahaye

I would like to understand this as well.

In the case where the simulations do converge, does the computed solution make physical sense? And if it does, do you obtain convergence by flipping the sign, i.e., solve == + fvm::Sp(...).

Same question in case where simulations diverge. Should you flip the sign? Or should you replace the function fvm::Sp() by fvm::SuSp()? See e.g. this post fvMatrix, fvOptions and SuSp: automatic implicit/explicit source-term treatment.

How good is converge in case you do converge?

More advanced, and possibly not worth the effort: can you run in a debugger, print the matrix to file and check diagonal dominance of the matrix?

Let us know.

[Anup Singh]

I think you need to linearize your equation first or look you into what you have written, as for RHS you are putting Np^2 both as value to be calculated and value to be used.
As for why its working for second implementation, I think that might be because the matrix is formed with the declaration done with fvScalarMatrix which is only modified later on with solve function but as the implementation is incorrect (As it appears to me) you are getting nonphysical results.

[adib_mohammad]

I think the currect implementation of the source term (in both cases) can be considered as a standard linearization.

-fvm::Sp(0.5 * beta * N_agg * Av * rho, N_agg))

I got the idea of the second implementation from U equation (UEqn.h) in simpleFoam. Please see how the gradient of pressure is treated in U equation.

    if (simple.momentumPredictor())
    {
        solve(UEqn == -fvc::grad(p));

        fvOptions.correct(U);
    }

Originally Posted by Anup Singh

I think you need to linearize your equation first or look you into what you have written, as for RHS you are putting Np^2 both as value to be calculated and value to be used.
As for why its working for second implementation, I think that might be because the matrix is formed with the declaration done with fvScalarMatrix which is only modified later on with solve function but as the implementation is incorrect (As it appears to me) you are getting nonphysical results.

[dlahaye]

@adib_mohammad:
pls. observe that your implementation employs fvm (implicit, method) while the simpeFoam example you refer to employs fvc (explit, calculus).
question: do you have a reference solution in literature that one could employ as a reference? Thx.

@mambojambo:
for the boundary contribution. I found this post ldUMatrix building indexing in parallel to explain the matter well.

@anup:
could you pls. further elaborate your point of view on this matter. Thx.

[adib_mohammad]

I see your point about fvc and fvm.
I have a refernece analytical solution for a case without the source term, and I have validated the numerical results with that. It is basically a monodisperse population balance model. The problem occurs when I add coagulation source term.

Originally Posted by dlahaye

@adib_mohammad:
pls. observe that your implementation employs fvm (implicit, method) while the simpeFoam example you refer to employs fvc (explit, calculus).
question: do you have a reference solution in literature that one could employ as a reference? Thx.

@mambojambo:
for the boundary contribution. I found this post ldUMatrix building indexing in parallel to explain the matter well.

@anup:
could you pls. further elaborate your point of view on this matter. Thx.

[Anup Singh]

adib .. I get your point but the thing is that you have to be clear on how to use your source terms in its implementation..

you can either employ your source term entirely explicitly ie by using fvc...
or employ your source term implicitly ie. by using fvm..

Though you have to keep into account that you need to linearize your source term before its implementation in case you are implementing it implicitly...

For the implicit implementation of your source term, I would suggest you refer to the source term linearization by SV Patanker..

or there are other examples as well...

for example, you can look into the implementation of NS equation in OpenFOAM...
you can notice that the div term should be div(u, u) but it is implemented as div(phi, u) which actually linearizes the non-linear term by taking phi as know value...

you also have to implement the same approach in your problem....

As for how to visualize it ...

I assume you know how that the final solution matrix takes the form
AX = B

Then you can roughly assume that fvm contributes to the formation of the coefficient matrix A whereas fvc contributes the formation of matrix B..

Also if you still have some doubts you can simply discretize your equation and see what you get out of it..

[MamboJambo]

Originally Posted by Anup Singh

you can either employ your source term entirely explicitly ie by using fvc...
or employ your source term implicitly ie. by using fvm..

Just a complement here. If you want to implement the source term explicitly you just need to equate it to a DimensionedField<Type, volMesh>. You do not need to call any particular function from fvc.

E.g.

            fvScalarMatrix TEqn
            (
                fvm::ddt(T) - fvm::laplacian(DT, T)
             ==
                someVolScalarField
            );

[dlahaye]

@adib_mohammad:
can you please share note on the analytical reference solution? We can try to extend the model with a linear source term first and a non-linear source term subsequently.

Cheers.

[Anup Singh]

Yup, that's true in case you have already calculated your source term but if you have to calculate it at the instance then you have to use fvc utility to calculate on the spot which also saves space for you.

[adib_mohammad]

@dlahaye

The number density equation (mentioned in the post) without the source term is used to describe particle deposition on pipe walls, which is a widely used engineering. You can find analytical solution in "Smoke, Dust, and Haze: Fundamentals of Aerosol Behavior" or "Aerosol Technology: Properties, Behavior, and Measurement of Airborne Particles". The so-called peneration curves are plotted used the analytical solution of this equation (penetration tells you the pecentage of initial number of particles that passes through a channel without being deposited on the walls).

[dlahaye]

Thanks. How do you expect the plot to change in case that the source term is added? I.e., what is the influence of the additional source term.