Inviscid flow simulation problem

[sunny7bit]

Dear author, when I follow the tutorial (https://hystrath.github.io/guides/fleming/cfd/transport/#1-species-shear-viscosity-and-thermal-conductivity) for inviscid simulation, the simulation always gets stuck on the first iteration. I am wondering what else is needed besides setting the viscosity mu to 0. I am looking forward to your reply.

The lastest log information in the log.myhy2Foam：

Mean and max Courant Numbers = 0.00719930935881 0.5 deltaT = 3.98690608802e-09 Time = 3.986906088e-09

diagonal: Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0 diagonal: Solving for rhoUx, Initial residual = 0, Final residual = 0, No Iterations 0 diagonal: Solving for rhoUy, Initial residual = 0, Final residual = 0, No Iterations 0 diagonal: Solving for rhoUz, Initial residual = 0, Final residual = 0, No Iterations 0 diagonal: Solving for rhoE, Initial residual = 0, Final residual = 0, No Iterations 0

[sunny7bit]

This situation only happens when I edit the transport entry to constant in the thermophysicalProperties/thermoType dictionary. If I edit the transport entry to other models, the case will run successfully.

[fagabz]

Hi sunny7bit, I'm going through the same issue, were you able to solve it?

[LDK29]

Hi sunny7bit, I'm also going through the same issue, were you able to solve it?

[sunny7bit]

I’m sorry for not replying to you for such a long time. In fact, when we use the boundary condition such as nonEqSmoluchowskiJumpT and nonEqMaxwellSlipU, this problem will happen. This is due to a division by zero issue in the code for these boundary conditions. For example, in the folder applications\solvers\compressible\hy2Foam\BCs\U\nonEqMaxwellSlipUFvPatchVectorField.C, the following code exists:

void Foam::nonEqMaxwellSlipUFvPatchVectorField::updateCoeffs()
{
    if (updated())
    {
        return;
    }

const fvPatchScalarField& pmu =
        patch().lookupPatchField<volScalarField, scalar>(muName_);

    const fvPatchScalarField& prho =
        patch().lookupPatchField<volScalarField, scalar>(rhoName_);
    const fvPatchScalarField& pmfp =
        patch().lookupPatchField<volScalarField, scalar>(mfpName_); // NEW VINCENT 28/02/2016

    Field<scalar> pnu(pmu/prho);   

    // DELETION VINCENT 28/02/2016 ********************************************
    /*Field<scalar> C1
    (
        sqrt(ppsi*constant::mathematical::piByTwo)
      * (2.0 - accommodationCoeff_)/accommodationCoeff_
    );*/
    // END DELETION VINCENT 28/02/2016 ****************************************
    Field<scalar> C1
    (
        pmfp/pnu
      * (2.0 - accommodationCoeff_)/accommodationCoeff_
    ); // NEW VINCENT 28/02/2016: divided by nu to keep the same units as before

Notice these in the above code.

const fvPatchScalarField& pmu =
        patch().lookupPatchField<volScalarField, scalar>(muName_);

Field<scalar> pnu(pmu/prho);

pmfp/pnu
      * (2.0 - accommodationCoeff_)/accommodationCoeff_

/pnu is a divide-by-zero behavior if we use inviscid condition. I don’t know much about this kind of boundary, and perhaps this boundary itself is not suitable for inviscid situation. If that’s the case, an error output can actually be added to the code. It would report an error when this type of boundary condition is used along with inviscid condition. We can add a judgment and error output by mimicking the existing code.

Foam::nonEqMaxwellSlipUFvPatchVectorField::nonEqMaxwellSlipUFvPatchVectorField
(
    const fvPatch& p,
    const DimensionedField<vector, volMesh>& iF,
    const dictionary& dict
)
:
    mixedFixedValueSlipFvPatchVectorField(p, iF),
    TName_(dict.lookupOrDefault<word>("Tt", "Tt")),
    rhoName_(dict.lookupOrDefault<word>("rho", "rho")),
    muName_(dict.lookupOrDefault<word>("mu", "mu")),
    tauMCName_(dict.lookupOrDefault<word>("tauMC", "tauMC")),
    mfpName_(dict.lookupOrDefault<word>("mfp", "mfp")), // NEW VINCENT 28/02/2016
    **accommodationCoeff_(readScalar(dict.lookup("accommodationCoeff"))),**
    vectorUwall_(dict.lookupOrDefault<vector>("Uwall", vector::zero)), // NEW VINCENT 21/10/2016
    //Uwall_("Uwall_", dict, p.size()),
    Uwall_(p.size(), vectorUwall_), // NEW VINCENT 21/10/2016
    thermalCreep_(dict.lookupOrDefault("thermalCreep", true)),
    curvature_(dict.lookupOrDefault("curvature", true))
{
    if
    (
        mag(accommodationCoeff_) < SMALL
     || mag(accommodationCoeff_) > 1.0
    )
    {
        FatalIOErrorIn
        (
            "maxwellSlipUFvPatchScalarField::maxwellSlipUFvPatchScalarField"
            "("
                "const fvPatch&, "
                "const DimensionedField<vector, volMesh>&, "
                "const dictionary&"
            ")",
            dict
        )   << "unphysical accommodationCoeff_ specified"
            << "(0 < accommodationCoeff_ <= 1)" << endl
            << exit(FatalIOError);
    }

[LDK29]

I’m sorry for not replying to you for such a long time. In fact, when we use the boundary condition such as nonEqSmoluchowskiJumpT and nonEqMaxwellSlipU, this problem will happen. This is due to a division by zero issue in the code for these boundary conditions. For example, in the folder applications\solvers\compressible\hy2Foam\BCs\U\nonEqMaxwellSlipUFvPatchVectorField.C, the following code exists:

void Foam::nonEqMaxwellSlipUFvPatchVectorField::updateCoeffs() { if (updated()) { return; }

**const fvPatchScalarField& pmu =
    patch().lookupPatchField<volScalarField, scalar>(muName_);**
const fvPatchScalarField& prho =
    patch().lookupPatchField<volScalarField, scalar>(rhoName_);
const fvPatchScalarField& pmfp =
    patch().lookupPatchField<volScalarField, scalar>(mfpName_); // NEW VINCENT 28/02/2016

**Field<scalar> pnu(pmu/prho);**

// DELETION VINCENT 28/02/2016 ********************************************
/*Field<scalar> C1
(
    sqrt(ppsi*constant::mathematical::piByTwo)
  * (2.0 - accommodationCoeff_)/accommodationCoeff_
);*/
// END DELETION VINCENT 28/02/2016 ****************************************
Field<scalar> C1
(
    **pmfp/pnu**
  * (2.0 - accommodationCoeff_)/accommodationCoeff_
); // NEW VINCENT 28/02/2016: divided by nu to keep the same units as before

I don’t know much about this kind of boundary, and perhaps this boundary itself is not suitable for inviscid situation. If that’s the case, an error output can actually be added to the code. It would report an error when this type of boundary condition is used along with inviscid condition. We can add a judgment and error output by mimicking the existing code.

Foam::nonEqMaxwellSlipUFvPatchVectorField::nonEqMaxwellSlipUFvPatchVectorField ( const fvPatch& p, const DimensionedField<vector, volMesh>& iF, const dictionary& dict ) : mixedFixedValueSlipFvPatchVectorField(p, iF), TName_(dict.lookupOrDefault("Tt", "Tt")), rhoName_(dict.lookupOrDefault("rho", "rho")), muName_(dict.lookupOrDefault("mu", "mu")), tauMCName_(dict.lookupOrDefault("tauMC", "tauMC")), mfpName_(dict.lookupOrDefault("mfp", "mfp")), // NEW VINCENT 28/02/2016 accommodationCoeff_(readScalar(dict.lookup("accommodationCoeff"))), vectorUwall_(dict.lookupOrDefault("Uwall", vector::zero)), // NEW VINCENT 21/10/2016 //Uwall_("Uwall_", dict, p.size()), Uwall_(p.size(), vectorUwall_), // NEW VINCENT 21/10/2016 thermalCreep_(dict.lookupOrDefault("thermalCreep", true)), curvature_(dict.lookupOrDefault("curvature", true)) { if ( mag(accommodationCoeff_) < SMALL || mag(accommodationCoeff_) > 1.0 ) { FatalIOErrorIn ( "maxwellSlipUFvPatchScalarField::maxwellSlipUFvPatchScalarField" "(" "const fvPatch&, " "const DimensionedField<vector, volMesh>&, " "const dictionary&" ")", dict ) << "unphysical accommodationCoeff_ specified" << "(0 < accommodationCoeff_ <= 1)" << endl << exit(FatalIOError); }

Many thanks for your kindky reply! But I didin't used these two boundary in my simulation. I met this problem when a simple slip boundary was used. Assign a very small value for dynamic viscosity can partially solve this problem. And the results are reasonable for me.

[sunny7bit]

That feels a bit strange. I used the bluntedCone case that comes with hy2Foam, modified the temperature boundary to a fixed temperature wall, and changed the velocity to a slip wall. I found that it can be calculated normally.