[mecbe2002]

I am trying to add mass source using "scalarCodedSource" but the simulation diverges in first iteration itself.


Unit of MFR added here is kg/s
Solver i am using is "rhoPimpleFoam"


I tried with "scalarSemiImplicitSource" method using hard coded MFR value it works fine. But for my application i have to calculate MFR real time so I am using "codedSource".


It can be seen that pressure at the location where "mass Source" is being added goes negative.


Please let me know what is wrong with the code.

1. Unit of MFR is correct kg/s ?
2. source addition per cell volume? If so then MFR/cell volume has to be added as source?
3. filedName "rho" is correct for mass source?

massSourceEffusionHole
{
    type    scalarCodedSource;
    name    sourceTime0; 
    active  true;
    
    scalarCodedSourceCoeffs
    {            
        
        // Mass source variable
        fields          (rho);
        selectionMode   points;
           
        points
        (                                                                       
            // Sink location i.e. effusion hole inlet
            (0.035 0.0335 0.0)
            // Source location i.e. effusion hole exit
            (0.0385 0.0365 0.0)
        );
        
        codeInclude
        #{     
        #};
        
        codeCorrect
        #{                
        #};
        
        codeAddSup
        #{
            Pout<<"Mass Source Start"<<endl;
            const Time& time = mesh().time();
            const volScalarField& rhoJet = mesh_.lookupObject<volScalarField>("rho");
            const volVectorField& velJet = mesh_.lookupObject<volVectorField>("U");            
            const volScalarField& pJet = mesh_.lookupObject<volScalarField>("p");
            const scalarField& cellVolume = mesh_.V();
            
            // Adding mass source term
            scalarField& massSourceEffusionHole = eqn.source();
                                            
            // Finding the cell corresponding to sink i.e. effusion hole inlet
            label cellLabelSink = mesh_.findCell(points_[0]);
            // Finding the cell corresponding to source i.e. effusion hole exit
            label cellLabelSource = mesh_.findCell(points_[1]);      
            
            // Calculation of total pressure at Sink
            scalar totalPJet = pJet[cellLabelSink]+(0.5*rhoJet[cellLabelSink]*pow(mag(velJet[cellLabelSink]),2));
            Pout<<"Total Pressure at Sink = "<<totalPJet<<endl;
            Pout<<"Pressure Source = "<<pJet[cellLabelSource]<<endl;
            
            // Calculation of MFR at Sink
            scalar mfrSink = totalPJet*pow((pJet[cellLabelSource]/totalPJet),0.8571)*1.3212*1e-6*sqrt(((2.0*1.4)/(0.4*287.05*300))*(pow((totalPJet/pJet[cellLabelSource]),0.2857)-1));
                        
            
            Pout<<"MFR Sink = "<<mfrSink<<endl;       
            massSourceEffusionHole[cellLabelSource] += mfrSink;
            Pout<<"MFR Sink = "<<massSourceEffusionHole[cellLabelSource]<<endl;
                       
            Pout<<"Mass Source End"<<endl;
        #};
        
        codeSetValue
        #{        
        #}
        
        // Dummy entry. Make dependent on above to trigger recompilation code
        #{
            $codeInclude
            $codeCorrect
            $codeAddSup
            $codeSetValue
        #};
    }
    
    sourceTime0Coeffs
    {
        $scalarCodedSourceCoeffs;
    }
}

Time = 1

Using dynamicCode for fvOption:: sourceTime0 at line 72 in "/home/system/Projects/Siragugal/GasTubineCombustor/EffusionHoleModel/singleHoleFlatPlate_EffusionBC/constant/fvOptions.massSourceEffusionHole.scalarCodedSourceCoeffs"
Selecting finite volume options model type sourceTime0
    Source: sourceTime0
    - selecting cells using points
    - selected 2 cell(s) with volume 2.3806788e-10
Mass Source Start
Total Pressure at Sink = 220796.55
Pressure Source = 200000
MFR Sink = 0.00040910048
MFR Sink = 0.00040910048
Mass Source End
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
smoothSolver:  Solving for Ux, Initial residual = 0.0032980687, Final residual = 8.0245382e-07, No Iterations 16
smoothSolver:  Solving for Uy, Initial residual = 0.99999999, Final residual = 8.8311766e-07, No Iterations 23
smoothSolver:  Solving for Uz, Initial residual = 0.99999999, Final residual = 6.778778e-07, No Iterations 15
smoothSolver:  Solving for e, Initial residual = 4.2580813e-05, Final residual = 3.1907154e-07, No Iterations 1
Mass Source Start
Total Pressure at Sink = 220796.55
Pressure Source = 200000
MFR Sink = 0.00040910048
MFR Sink = 0.00040910048
Mass Source End
GAMG:  Solving for p, Initial residual = 0.00022212445, Final residual = 1.8957863e-05, No Iterations 17
GAMG:  Solving for p, Initial residual = 0.00029171794, Final residual = 2.5683858e-05, No Iterations 2
Mass Source Start
Total Pressure at Sink = 220796.54
Pressure Source = 61287.67
MFR Sink = 0.0005830456
MFR Sink = 0.0005830456
Mass Source End
diagonal:  Solving for rho, Initial residual = 0, Final residual = 0, No Iterations 0
time step continuity errors : sum local = 1.0513006e-05, global = -7.3666193e-06, cumulative = -7.3666193e-06
Mass Source Start
Total Pressure at Sink = 220796.56
Pressure Source = 61287.67
MFR Sink = 0.00058304562
MFR Sink = 0.00058304562
Mass Source End
GAMG:  Solving for p, Initial residual = 0.00027760611, Final residual = 2.3000663e-05, No Iterations 11
GAMG:  Solving for p, Initial residual = 0.00017927573, Final residual = 9.6137966e-07, No Iterations 31
Mass Source Start
Total Pressure at Sink = 220796.57
Pressure Source = -50175.002


#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigFpe::sigHandler(int) at ??:?
#2  ? in "/lib/x86_64-linux-gnu/libc.so.6"
#3  ? in "/lib/x86_64-linux-gnu/libm.so.6"
#4  pow in "/lib/x86_64-linux-gnu/libm.so.6"
#5  Foam::fv::sourceTime0FvOptionscalarSource::addSup(Foam::fvMatrix<double>&, int) at constant/fvOptions.massSourceEffusionHole.scalarCodedSourceCoeffs:114
#6  ? in "/home/system/OpenFOAM/OpenFOAM-6/platforms/linux64GccDPInt32Opt/bin/rhoPimpleFoam"
#7  ? in "/home/system/OpenFOAM/OpenFOAM-6/platforms/linux64GccDPInt32Opt/bin/rhoPimpleFoam"
#8  ? in "/home/system/OpenFOAM/OpenFOAM-6/platforms/linux64GccDPInt32Opt/bin/rhoPimpleFoam"
#9  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#10  ? in "/home/system/OpenFOAM/OpenFOAM-6/platforms/linux64GccDPInt32Opt/bin/rhoPimpleFoam"
Floating point exception (core dumped)

[mecbe2002]

Still I am having the same problem.


Is there any specific code to be added for "rhoPimpleFoam" solver?


Simulation diverging in first iteration itself.. bcz of pressure going negative at the cell where mass source is being added..