[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);
}