[wang219910611]

Hi Everyone,

Recently I am looking at the KHRT breakup model in OpenFoam, located at src/lagrangian/spray/submodels/BreakupModel.

While I am looking at it, I find out there is something different of the KHRT applied in OpenFoam and the one applied in KIVA. That is, there is no breakup up length concept in Openfoam, while in KIVA, within the breakup length Lb, only KH breakup can happen. This might lead to more RT breakup in OpenFoam than it should be.

Now I want to apply this limitation in Openfoam. But I am not should how to set the variable to indicate the length from droplets to the injector. Or is it simply x....?

Anyone gets any ideas here? Many thanks!

[wang219910611]

Hi everyone, I think I have got the answer. It seems this function

p.position().x() < Lb

would be appropriate for this.

However, I meet another problem before I could solve this one. In my opinion, the breakup model should iterate though all the parcels to finish its job. But I can't find the iteration process in the BreakupModel files. Nor can I find it in the SprayCloud class. So, where is it?

Any suggestion will be appreciated, many thanks.

[mturcios777]

Hi wang,

Look at the penetration funcion in Kinematic Cloud object, as .x() is on the the x coordinate of the current parcel position.

[wang219910611]

Dear Marco,
First of all, thank your for the suggestion! After looking at KinematicCloudl.H, I believe I should use this variable as breakup length indictor

scalar dis = mag(this->position()-this->position0());

And I tried to add this variable the same way as Urel is applied:
First, in SprayParcel.C

    
template<class ParcelType>
template<class TrackData>
void Foam::SprayParcel<ParcelType>::calcBreakup
(
   ......
   ......
    scalar rhoAv = this->pc()/(R*Tav);
    scalar muAv = this->muc();
    vector Urel = this->U() - this->Uc();
    scalar dis = mag(this->position() - this->position0());
    scalar Urmag = mag(Urel);
    scalar Re = this->Re(this->U(), this->d(), rhoAv, muAv);
    .......
    .......
td.cloud().breakup().update
        (
           ....
           ....
            this->U(),
            rhoAv,
            muAv,
            Urel,
            dis,
            Urmag,
            this->tMom(),
            dChild,
            massChild
        )

And then I add the "dis" in my BreakupModel.C,EKHRT.H and EKHRT.C file accordingly. However, although I am able to compile the EKHRT model without any problem, when I tried to run the aachenBomb case ( modified accordingly), it failed as soon as the second time step(when particle are injected). Could you please kindly provide some clues which steps I have done wrong? I have attached the main codes in the attachment.

[wang219910611]

Solving 3-D cloud sprayCloud

Program received signal SIGFPE, Arithmetic exception.
0x00000000004b81c0 in Foam::divideOp3<double, double, double>::operator() (this=0x7fffffff7770, x=@0x7fffffff7b30: -0.59621378633365985, 
    y=@0x7fffffff7778: 6.9533558061496476e-310) at /home/bo/OpenFOAM/OpenFOAM-2.3.0/src/OpenFOAM/lnInclude/ops.H:167
warning: Source file is more recent than executable.
167    Op(divide, x / y+0.00001)
(gdb) bt
#0  0x00000000004b81c0 in Foam::divideOp3<double, double, double>::operator() (this=0x7fffffff7770, x=@0x7fffffff7b30: -0.59621378633365985, 
    y=@0x7fffffff7778: 6.9533558061496476e-310) at /home/bo/OpenFOAM/OpenFOAM-2.3.0/src/OpenFOAM/lnInclude/ops.H:167
#1  0x000000000049efd1 in VectorSpaceOps<3, 0>::opVS<Foam::Vector<double>, Foam::VectorSpace<Foam::Vector<double>, double, 3>, double, Foam::divideOp3<double, double, double> > (vs=..., vs1=..., s=@0x7fffffff7778: 6.9533558061496476e-310, o=...)
    at /home/bo/OpenFOAM/OpenFOAM-2.3.0/src/OpenFOAM/lnInclude/VectorSpaceM.H:34
#2  0x000000000048cdbd in Foam::operator/<Foam::Vector<double>, double, 3> (vs=..., s=6.9533558061496476e-310)
    at /home/bo/OpenFOAM/OpenFOAM-2.3.0/src/OpenFOAM/lnInclude/VectorSpaceI.H:585
#3  0x00007fffe9a7a365 in Foam::EKHRT<Foam::SprayCloud<Foam::ReactingCloud<Foam::ThermoCloud<Foam::KinematicCloud<Foam::Cloud<Foam::SprayParcel<Foam::ReactingParcel<Foam::ThermoParcel<Foam::KinematicParcel<Foam::particle> > > > > > > > > >::update (this=0x466dbe0, 
    dt=4.313064053868374e-07, g=..., d=@0xb957af8: 3.7789377622424555e-05, tc=@0xb957c10: 0, ms=@0xb957c18: 0, 
    nParticle=@0xb957af0: 2.0142790546900065, KHindex=@0xb957bf8: 0, y=@0xb957c00: 0, yDot=@0xb957c08: 0, d0=3.7779345486647244e-05, 
    rho=662.15416587450784, mu=0.00031014123604303163, sigma=0.017555291094985808, U=..., rhoc=36.119402030502705, 
    muc=3.8978930586897314e-05, Urel=..., pos=..., injectionPos=..., Urmag=34.888565560017959, tMom=6.9533558061496476e-310, 
    dChild=@0x3e4460aa68000000: <error reading variable>, massChild=@0xb957c10: 0) at lnInclude/EKHRT.C:130
#4  0x000000000050db45 in Foam::SprayParcel<Foam::ReactingParcel<Foam::ThermoParcel<Foam::KinematicParcel<Foam::particle> > > >::calcBreakup<Foam::ReactingParcel<Foam::ThermoParcel<Foam::KinematicParcel<Foam::particle> > >::TrackingData<Foam::SprayCloud<Foam::ReactingCloud<Foam::ThermoCloud<Foam::KinematicCloud<Foam::Cloud<Foam::SprayParcel<Foam::ReactingParcel<Foam::ThermoParcel<Foam::KinematicParcel<Foam::particle> > > > > > > > > > > (this=0xb957a90, td=..., dt=4.313064053868374e-07, cellI=167259)
    at /home/bo/OpenFOAM/OpenFOAM-2.3.0/src/lagrangian/spray/lnInclude/SprayParcel.C:253
#5  0x00000000004fd682 in Foam::SprayParcel<Foam::ReactingParcel<Foam::ThermoParcel<Foam::KinematicParcel<Foam::particle> > > >::calc<Foam::ReactingParcel<Foam::ThermoParcel<Foam::KinematicParcel<Foam::particle> > >::TrackingData<Foam::SprayCloud<Foam::ReactingCloud<Foam::ThermoCloud<Foam::KinematicCloud<Foam::Cloud<Foam::SprayParcel<Foam::ReactingParcel<Foam::ThermoParcel<Foam::KinematicParcel<Foam::particle> > > > > > > > > > > (this=0xb957a90, td=..., dt=4.313064053868374e-07, cellI=167259)
    at /home/bo/OpenFOAM/OpenFOAM-2.3.0/src/lagrangian/spray/lnInclude/SprayParcel.C:135

I tried to use gdb to debug with the problem, then above error showed up. It seemed to me that this is a "divided by zero" problem. However, while I was trying to identify this according to the information provided, I am lost... could anyone please give any suggestions here? Thanks in advance.

[mturcios777]

If this is the first timestep, then it is possible that dis is still 0 for a parcel that hasn't moved yet (this depends on your timestep and lots of other parameters). It seems like the mag function is having issues, which doesn't make sense as it looks equipped to handle the zero vector.

Can you try find out if the call to mag is working? Otherwise it may be that when you divide by dis the error occurs. That is a bit simpler fix, as you can do x/(y + VSMALL) or similar. I would use the built in types like SMALL and VSMALL. Check your parentheses as well, as it looks like the operation is trying to do x/y + 0.00001 which is probably not what you want.

[wang219910611]

Hi Marco,

Firstly, thank you very much for your kindly suggestion. I tried firstly to leave mag function out of this problem by using a constant instead:

 scalar lambdaRT = constant::mathematical::twoPi*cRT_/(KRT + VSMALL);
 //   scalar diff = mag(pos - injectionPos); //here I left mag function out
    scalar diff=0.5;
    // if lambdaRT < diameter, then RT waves are growing 
    if (diff > Lbreakup_) // Lbreakup_ is set to 0.1 here
{
    if ((tc > 0) || (lambdaRT < d) )
    {
        tc += dt;
    }
    // characteristic RT breakup time
   
          .......

However, the error is still there (although the left code are exactly the same as original KHRT model...). :

Solving 3-D cloud sprayCloud

Program received signal SIGFPE, Arithmetic exception.
0x00000000004b81c0 in Foam::divideOp3<double, double, double>::operator() (
    this=0x7fffffff7770, x=@0x7fffffff7b30: -0.59621378633365985, 
    y=@0x7fffffff7778: 6.9533558061496476e-310)
    at /home/bo/OpenFOAM/OpenFOAM-2.3.0/src/OpenFOAM/lnInclude/ops.H:167
warning: Source file is more recent than executable.
167    Op(divide, x / y)
(gdb) bt
#0  0x00000000004b81c0 in Foam::divideOp3<double, double, double>::operator() (
    this=0x7fffffff7770, x=@0x7fffffff7b30: -0.59621378633365985, 
    y=@0x7fffffff7778: 6.9533558061496476e-310)
    at /home/bo/OpenFOAM/OpenFOAM-2.3.0/src/OpenFOAM/lnInclude/ops.H:167
#1  0x000000000049efd1 in VectorSpaceOps<3, 0>::opVS<Foam::Vector<double>, Foam::VectorSpace<Foam::Vector<double>, double, 3>, double, Foam::divideOp3<double, double, double> > (vs=..., vs1=..., s=@0x7fffffff7778: 6.9533558061496476e-310, 
    o=...)
    at /home/bo/OpenFOAM/OpenFOAM-2.3.0/src/OpenFOAM/lnInclude/VectorSpaceM.H:34
#2  0x000000000048cdbd in Foam::operator/<Foam::Vector<double>, double, 3> (
    vs=..., s=6.9533558061496476e-310)
    at /home/bo/OpenFOAM/OpenFOAM-2.3.0/src/OpenFOAM/lnInclude/VectorSpaceI.H:585
#3  0x00007fffe9a7a365 in Foam::EKHRT<Foam::SprayCloud<Foam::ReactingCloud<Foam::ThermoCloud<Foam::KinematicCloud<Foam::Cloud<Foam::SprayParcel<Foam::ReactingParcel<Foam::ThermoParcel<Foam::KinematicParcel<Foam::particle> > > > > > > > > >::update (this=0x466dbe0, dt=4.313064053868374e-07, g=..., 
    d=@0xb957af8: 3.7789377622424555e-05, tc=@0xb957c10: 0, ms=@0xb957c18: 0, 
    nParticle=@0xb957af0: 2.0142790546900065, KHindex=@0xb957bf8: 0, 
    y=@0xb957c00: 0, yDot=@0xb957c08: 0, d0=3.7779345486647244e-05, 
    rho=662.15416587450784, mu=0.00031014123604303163, 
    sigma=0.017555291094985808, U=..., rhoc=36.119402030502705, 
---Type <return> to continue, or q <return> to quit---

Then I tried to adapt the ops.H file as Op(divide, x / (y+VSMALL)), but still error. By the way, the calculation goes well for the first time step, but as so as it tries to calculate the 3D spray cloud in the second, things mess up.

[wang219910611]

After I have commented out all the changes I made in the EKHRT files(except for "pos" and "position0" in the update function, error still exits:

Solving 3-D cloud sprayCloud

Program received signal SIGFPE, Arithmetic exception.
0x00000000004b81c0 in Foam::divideOp3<double, double, double>::operator() (this=0x7fffffff7770, 
    x=@0x7fffffff7b30: -0.59621378633365985, y=@0x7fffffff7778: 6.9533558061496476e-310)
    at /home/bo/OpenFOAM/OpenFOAM-2.3.0/src/OpenFOAM/lnInclude/ops.H:167
warning: Source file is more recent than executable.
167    Op(divide, x / (y+VSMALL))
(gdb) bt
#0  0x00000000004b81c0 in Foam::divideOp3<double, double, double>::operator() (this=0x7fffffff7770, 
    x=@0x7fffffff7b30: -0.59621378633365985, y=@0x7fffffff7778: 6.9533558061496476e-310)
    at /home/bo/OpenFOAM/OpenFOAM-2.3.0/src/OpenFOAM/lnInclude/ops.H:167
#1  0x000000000049efd1 in VectorSpaceOps<3, 0>::opVS<Foam::Vector<double>, Foam::VectorSpace<Foam::Vector<double>, double, 3>, double, Foam::divideOp3<double, double, double> > (vs=..., vs1=..., s=@0x7fffffff7778: 6.9533558061496476e-310, o=...)
    at /home/bo/OpenFOAM/OpenFOAM-2.3.0/src/OpenFOAM/lnInclude/VectorSpaceM.H:34
#2  0x000000000048cdbd in Foam::operator/<Foam::Vector<double>, double, 3> (vs=..., s=6.9533558061496476e-310)
    at /home/bo/OpenFOAM/OpenFOAM-2.3.0/src/OpenFOAM/lnInclude/VectorSpaceI.H:585
#3  0x00007fffe9a7a375 in Foam::EKHRT<Foam::SprayCloud<Foam::ReactingCloud<Foam::ThermoCloud<Foam::KinematicCloud<Foam::Cloud<Foam::SprayParcel<Foam::ReactingParcel<Foam::ThermoParcel<Foam::KinematicParcel<Foam::particle> > > > > > > > > >::update
    (this=0x466dbe0, dt=4.313064053868374e-07, g=..., d=@0xb957af8: 3.7789377622424555e-05, tc=@0xb957c10: 0, 
    ms=@0xb957c18: 0, nParticle=@0xb957af0: 2.0142790546900065, KHindex=@0xb957bf8: 0, y=@0xb957c00: 0, 
    yDot=@0xb957c08: 0, d0=3.7779345486647244e-05, rho=662.15416587450784, mu=0.00031014123604303163, 
    sigma=0.017555291094985808, U=..., rhoc=36.119402030502705, muc=3.8978930586897314e-05, Urel=..., pos=..., 
    injectionPos=..., Urmag=34.888565560017959, tMom=6.9533558061496476e-310, 
    dChild=@0x3e4460aa68000000: <error reading variable>, massChild=@0xb957c10: 0) at lnInclude/EKHRT.C:130

So maybe there is something wrong with my modified SprayParcel class? The only place I have modified is as following, in which I added "pos" and "position0" so I can calculate " scalar diff = mag(pos-position0)" in the EKHRT.C, I don't know where can the Arithmetic exception error happen...:

    ......
    scalar muAv = this->muc();
    const vector& pos = this->position();
    const vector& injectionPos = this->position0();
    scalar Urmag = mag(Urel);
    scalar Re = this->Re(this->U(), this->d(), rhoAv, muAv);
    ......
    const vector g = td.cloud().g().value();

    scalar massChild = 0.0;
    scalar dChild = 0.0;
    if
    (
        td.cloud().breakup().update
        (
            dt,
            g,
            ......
            Urel,
            pos,
            injectionPos,
            Urmag,
            this->tMom(),
            dChild,
            massChild
        )
    )

[mturcios777]

Hmm, that is pretty unusual. I'll try have a look at the code when I get some time later this week. If you are able, can you try run your modified case with the standard breakup models and post the results. If your new model has commented out differences and the old one runs properly, we can eliminate some possible sources of problems.

[wang219910611]

Hi, Marco, in fact, I didn't do much change in the case file except for the breakup model selection and the coefficient selection, where I just added the "Lbreakup" :

EKHRTCoeffs
    {
        solveOscillationEq yes;
        B0              0.61;
        B1              40;
        Ctau            1;
        CRT             0.1;
        msLimit         0.2;
        WeberLimit      6;
        Lbreakup       0.1;
    }

After I changed the breakup model back to the normal KHRT, case rans well.
Now, after I have commented the changes out, all the difference between EKHRT and KHRT is in the update function during breakup calculation, in SprayParcel.C and EKHRT.C and EKHRT.H ( and BreakupModel file , where I added "pos" and "position0"). However, I still couldn't figure where I have done wrong. If you could have a look at the code, I would be very appreciated. Thank you in advance. : )

[mturcios777]

I had a look at the classes and couldn't find anything unusual, though I did notice a few extra spaces at the end of the EKRT model directory. I try keep a debug compiled version of libraries around to spot things like this, as they usually give the line numbers where errors occur, but I also see you are running on a Mac. Can you post the results from a debug compiled version of your library?

[wang219910611]

Hi, Marco

In fact, I always run the Openfoam 230 on my Ubuntu desk PC rather than the mac.I just use mac to upload the thread because the Internet connection in my PC is a bit strange... Anyway, here is the result of Debug model with gdb ( with the code I uploaded in the last thread):

Solving 3-D cloud sprayCloud

Program received signal SIGFPE, Arithmetic exception.
0x00000000004b81c0 in Foam::divideOp3<double, double, double>::operator() (this=0x7fffffff7770, x=@0x7fffffff7b30: -0.59621378633365985, 
    y=@0x7fffffff7778: 6.9533558061496476e-310) at /home/bo/OpenFOAM/OpenFOAM-2.3.0/src/OpenFOAM/lnInclude/ops.H:167
warning: Source file is more recent than executable.
167    Op(divide, x / (y+VSMALL))
(gdb) bt
#0  0x00000000004b81c0 in Foam::divideOp3<double, double, double>::operator() (this=0x7fffffff7770, x=@0x7fffffff7b30: -0.59621378633365985, 
    y=@0x7fffffff7778: 6.9533558061496476e-310) at /home/bo/OpenFOAM/OpenFOAM-2.3.0/src/OpenFOAM/lnInclude/ops.H:167
#1  0x000000000049efd1 in VectorSpaceOps<3, 0>::opVS<Foam::Vector<double>, Foam::VectorSpace<Foam::Vector<double>, double, 3>, double, Foam::divideOp3<double, double, double> > (vs=..., vs1=..., s=@0x7fffffff7778: 6.9533558061496476e-310, o=...)
    at /home/bo/OpenFOAM/OpenFOAM-2.3.0/src/OpenFOAM/lnInclude/VectorSpaceM.H:34
#2  0x000000000048cdbd in Foam::operator/<Foam::Vector<double>, double, 3> (vs=..., s=6.9533558061496476e-310)
    at /home/bo/OpenFOAM/OpenFOAM-2.3.0/src/OpenFOAM/lnInclude/VectorSpaceI.H:585
#3  0x00007fffe9a7a375 in Foam::EKHRT<Foam::SprayCloud<Foam::ReactingCloud<Foam::ThermoCloud<Foam::KinematicCloud<Foam::Cloud<Foam::SprayParcel<Foam::ReactingParcel<Foam::ThermoParcel<Foam::KinematicParcel<Foam::particle> > > > > > > > > >::update (this=0x466dbe0, dt=4.313064053868374e-07, g=..., 
    d=@0xb957af8: 3.7789377622424555e-05, tc=@0xb957c10: 0, ms=@0xb957c18: 0, nParticle=@0xb957af0: 2.0142790546900065, KHindex=@0xb957bf8: 0, 
    y=@0xb957c00: 0, yDot=@0xb957c08: 0, d0=3.7779345486647244e-05, rho=662.15416587450784, mu=0.00031014123604303163, sigma=0.017555291094985808, 
    U=..., rhoc=36.119402030502705, muc=3.8978930586897314e-05, Urel=..., pos=..., injectionPos=..., Urmag=34.888565560017959, 
    tMom=6.9533558061496476e-310, dChild=@0x3e4460aa68000000: <error reading variable>, massChild=@0xb957c10: 0) at lnInclude/EKHRT.C:130
#4  0x000000000050db45 in Foam::SprayParcel<Foam::ReactingParcel<Foam::ThermoParcel<Foam::KinematicParcel<Foam::particle> > > >::calcBreakup<Foam::ReactingParcel<Foam::ThermoParcel<Foam::KinematicParcel<Foam::particle> > >::TrackingData<Foam::SprayCloud<Foam::ReactingCloud<Foam::ThermoCloud<Foam::KinematicCloud<Foam::Cloud<Foam::SprayParcel<Foam::ReactingParcel<Foam::ThermoParcel<Foam::KinematicParcel<Foam::particle> > > > > > > > > > > (this=0xb957a90, 
    td=..., dt=4.313064053868374e-07, cellI=167259) at /home/bo/OpenFOAM/OpenFOAM-2.3.0/src/lagrangian/spray/lnInclude/SprayParcel.C:253
#5  0x00000000004fd682 in Foam::SprayParcel<Foam::ReactingParcel<Foam::ThermoParcel<Foam::KinematicParcel<Foam::particle> > > >::calc<Foam::ReactingParcel<Foam::ThermoParcel<Foam::KinematicParcel<Foam::particle> > >::TrackingData<Foam::SprayCloud<Foam::ReactingCloud<Foam::ThermoCloud<Foam::KinematicCloud<Foam::Cloud<Foam::SprayParcel<Foam::ReactingParcel<Foam::ThermoParcel<Foam::KinematicParcel<Foam::particle> > > > > > > > > > > (this=0xb957a90, td=..., 
    dt=4.313064053868374e-07, cellI=167259) at /home/bo/OpenFOAM/OpenFOAM-2.3.0/src/lagrangian/spray/lnInclude/SprayParcel.C:135
#6  0x00000000004ec886 in Foam::KinematicParcel<Foam::particle>::move<Foam::ReactingParcel<Foam::ThermoParcel<Foam::KinematicParcel<Foam::particle> > >::TrackingData<Foam::SprayCloud<Foam::ReactingCloud<Foam::ThermoCloud<Foam::KinematicCloud<Foam::Cloud<Foam::SprayParcel<Foam::ReactingParcel<Foam::ThermoParcel<Foam::KinematicParcel<Foam::particle> > > > > > > > > > > (this=0xb957a90, td=..., trackTime=1.2040637150382544e-06)
    at /home/bo/OpenFOAM/OpenFOAM-2.3.0/src/lagrangian/intermediate/lnInclude/KinematicParcel.C:325

Then I changed the ops.H file back to the original, the results are exactly the same as shown above ( except for Op(divide, x/y).
And then if I run the case with normal KHRT model, the solver works well.

[mturcios777]

I just realized you had already posted what I needed. Sorry about that. I don't think its the calculation that is wrong, but the arguments being fed to it. In the stack trace, you can see the arguments of the functions where the errors are occuring, as well as the line number where it is happeing. In this case, have a look at #3:

#3  0x00007fffe9a7a375 in Foam::EKHRT<Foam::SprayCloud<Foam::ReactingCloud<Foam::ThermoCloud<Foam::KinematicCloud<Foam::Cloud<Foam::SprayParcel<Foam::ReactingParcel<Foam::ThermoParcel<Foam::KinematicParcel<Foam::particle> > > > > > > > > >::update (this=0x466dbe0, dt=4.313064053868374e-07, g=..., 
    d=@0xb957af8: 3.7789377622424555e-05, tc=@0xb957c10: 0, ms=@0xb957c18: 0, nParticle=@0xb957af0: 2.0142790546900065, KHindex=@0xb957bf8: 0, 
    y=@0xb957c00: 0, yDot=@0xb957c08: 0, d0=3.7779345486647244e-05, rho=662.15416587450784, mu=0.00031014123604303163, sigma=0.017555291094985808, 
    U=..., rhoc=36.119402030502705, muc=3.8978930586897314e-05, Urel=..., pos=..., injectionPos=..., Urmag=34.888565560017959, 
    tMom=6.9533558061496476e-310, dChild=@0x3e4460aa68000000: <error reading variable>, massChild=@0xb957c10: 0) at lnInclude/EKHRT.C:130

The code in question is:

vector = acceleration = Urel/tMom

The value of tMom is machine zero (E-310!), so that is why we get the error. When you changed the code in ops.H, did you recompile all the dependant librararies? Because you modified the base OpenFOAM library, you would need to recompile libOpenFOAM.so and anything else that depended on ops.H. However, I don't think this is what you want, since if x/(y+SMALL) or something similar is what you use, then calls to vector/scalar where scalar = 0 may work even when they shouldn't (ie unphysical behaviour).

The most sensible solution is to find out why tMom = 0 when you are using your model. It feels like maybe you are calling the update function or calculating the breakup too soon, as it seems tMom hasn't been initialized or calculated. The line numbers in your stack trace are always super helpful.

[wang219910611]

Hi Marco,

Thank you for your advice! At last I find this error is linked to the sequence I set up variables in the update function!

If in previous thread, I set up update function like this, that's when all the error happens:

td.cloud().breakup().update
        (
            dt,
            g,
            ...
            pos,
            injectionPos,
            Urmag,
            this->tMom(),
            dChild,
            massChild
       )

However, when I correct the function like this, errors are gone magically! :

 (
            dt,
            g,
            ......
            dChild,
            massChild,
            pos,
            injectionPos
       )

[wang219910611]

(Following last thread)
Then I tried to add a unused scalar diff=mag(pos-injectionPos) inEKHRT.C, error happens again.
Then I tried to change variables positions in "update function" according to diff position in EKHRT.C like this:

td.cloud().breakup().update
        (
            dt,
            g,
            this->d(),
            this->tc(),
            this->ms(),
            this->nParticle(),
            this->KHindex(),
            this->y(),
            this->yDot(),
            this->d0(),
            this->rho(),
            mu_,
            sigma_,
            this->U(),
            rhoAv,
            muAv,
            Urel,
            Urmag,
            this->tMom(),
            pos,
            injectionPos,
            dChild,
            massChild)

Error still, though different from last time.

Now I realize all the error comes from the sequence of variables in update function in SprayParcel.C and EKHRT.C. However, I am not sure how this happens still as I think I have currently put everything in the right order but there is still some errors. Could you please have a look at this? I have attached current ETAB.C and SprayParcel.C in this thread.

[wang219910611]

Hi, after a deeper look into the SprayParcel file, I think to distinguish primary and second breakup in the KHRT model, the better way is to apply the KH part into the atomization model. The liquidcore should be able to work as the length limit. Is it alright?
Another concern is that why does OpenFoam does't make such a distinguish after all...