View Issue Details
ID  Project  Category  View Status  Date Submitted  Last Update 

0003935  OpenFOAM  Bug  public  20221127 12:12  20221128 14:58 
Reporter  nucerl  Assigned To  henry  
Priority  normal  Severity  minor  Reproducibility  always 
Status  resolved  Resolution  fixed  
Platform  GNU/Linux  OS  Other  OS Version  (please specify) 
Product Version  10  
Fixed in Version  dev  
Summary  0003935: Discrepancy on the implementation of the IATE Turbulent Breakup source, Rti  
Description  1 Foam::tmp<Foam::fvScalarMatrix> 2 Foam::diameterModels::IATEsources::turbulentBreakUp::R 3 ( 4 const volScalarField& alphai, 5 volScalarField& kappai 6 ) const 7 { 8 volScalarField::Internal R 9 ( 10 IOobject 11 ( 12 "turbulentBreakUp:R", 13 iate_.phase().time().timeName(), 14 iate_.phase().mesh() 15 ), 16 iate_.phase().mesh(), 17 dimensionedScalar(kappai.dimensions()/dimTime, 0) 18 ); 19 20 const scalar Cti = Cti_.value(); 21 const scalar WeCr = WeCr_.value(); 22 const volScalarField Ut(this>Ut()); 23 const volScalarField We(this>We()); 24 25 forAll(R, celli) 26 { 27 if (We[celli] > WeCr) 28 { 29 R[celli] = 30 2*Cti*Ut[celli]*sqrt(1  WeCr/We[celli])*exp(WeCr/We[celli]); 31 } 32 } 33 34 return fvm::Su(R, kappai); 35 } The implementation of the only source term, Rti, is given in the snippet above. It can be seen that the defined 1/ms unit of Rti in line 17 of the snippet is not consistent with the implemented unit in line 30 which is 1/ms}. It can also be seen in line 34 of the snippet that the Su() function is utilized which treats the source term explicitly by modifying the source coefficients. The returned unit of the function is m/s*1/m = 1/s which is $not consistent with the other sink terms and the rest of the Interfacial Curvature Transport Equation (ICTE) implemented in OpenFOAM. The rate comparison is shown in the Snap#1 of the provided snaps after canceling out the long common terms in the reference and implemented equations. The effect of the abovementioned unit inconsistencies can be seen in the rate comparison which revealed that the implemented Rti is not consistent with the reference model [Ishii et. al, 2005]. Just considering the reference model, the required model can be derived as shown in Snap#2. It looks like the discrepancy comes from a missing length of a unit in the denominator (i.e. d: Sauter Mean Diameter) and a different constant. Changing the current model to the required one ensures that the local unit of Rti is consistent with the unit of other source terms which is 1/s and the unit of the returned value consequently becomes 1/ms which is also consistent with rest of the terms in the implemented ICTE. The discrepancy among the sources indicates an error unless it is fixed somewhere else in the code for Rti that I am not aware. It should also be noted the required form given in Snap#2 was also used in twoPhaseEulerFoam solver in OpenFOAM V7 (https://github.com/OpenFOAM/OpenFOAM7/blob/master/applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/diameterModels/IATE/IATEsources/turbulentBreakUp/turbulentBreakUp.C). However, the Su() function was not used to calculate Rti in OpenFOAM V7. I set up a test case to see the difference between the two models shown in Snap#2. The required model in Snap#2 is implemented to OpenFOAM V10 as shown in lines 24 and 32 of the snippet given in the Steps To Reproduce section (ignore the newIATE). However, the effect was negligible. Regardless of the implementation results, further investigation is needed to clarify TI modeling.  
Steps To Reproduce  1 Foam::tmp<Foam::fvScalarMatrix> 2 Foam::diameterModels::newIATEsources::turbulentBreakUp::R 3 ( 4 const volScalarField& alphai, 5 volScalarField& kappai 6 ) const 7 { 8 volScalarField::Internal R 9 ( 10 IOobject 11 ( 12 "turbulentBreakUp:R", 13 newiate_.phase().time().timeName(), 14 newiate_.phase().mesh() 15 ), 16 newiate_.phase().mesh(), 17 dimensionedScalar(kappai.dimensions()/dimTime, 0) 18 ); 19 20 const scalar Cti = Cti_.value(); 21 const scalar WeCr = WeCr_.value(); 22 const volScalarField Ut(this>Ut()); 23 const volScalarField We(this>We()); 24 const volScalarField& d(newiate_.d()()); //ebicer 25 26 forAll(R, celli) 27 { 28 if (We[celli] > WeCr) 29 { 30 R[celli] = 31 //2*Cti*Ut[celli]*sqrt(1  WeCr/We[celli])*exp(WeCr/We[celli]); 32 (1.0/3.0)*Cti/d[celli]*Ut[celli]*sqrt(1  WeCr/We[celli])*exp(WeCr/We[celli]); //ebicer 33 } 34 } 35 36 return fvm::Su(R, kappai); 37 }  
Tags  No tags attached.  



I agree that the dimensions are incorrect but by a factor of sqr(kappai) rather than 1/d (I am not sure what "d" is in the code snippet). Also according to eq. 11163 in Mamoru Ishii, Takashi Hibiki ThermoFluid Dynamics of TwoPhase Flow Second Edition the expression should be divided by 18 rather than 3 so I think the correct expression is: R[celli] = (Cti/18)*Ut[celli]*sqr(kappai[celli]) *sqrt(1  WeCr/We[celli])*exp(WeCr/We[celli]); 

Dear Henry, I am not sure how the unit sqr(kappai) missing since the local R (dimensionedScalar(kappai.dimensions()/dimTime, 0)) is 1/ms and the returned R (fvm::Su(R, kappai)) is in m/s*1/m = 1/s. kappai has the unit of 1/m. What am I missing here? All the literature of IATE is combined in Ishii's book, however the details of the derivations are usually missing. Although it is the same in onegroup IATE, the equation you are referring to is for twogroup IATE. I'd recommend checking Wu et. al, 1998 (Onegroup interfacial area transport in vertical bubbly flow) or Dr. Kim's thesis in 1999 (Interfacial Area Transport Equation and Measurement of Local Interfacial Characteristics) for more details. Back to the expression you provided, the 1/18 comes from how the sphericity (shape factor) is defined for assuming the bubbles to be in spherical shape, i.e. 1/36*PI(). However, please remember that the expression is for IATE not ICTE which is implemented in OpenFOAM. The "d" in the snippet refers to the Sauter mean diameter (shown below) which is updated after the kappaiEqn is solved in IATE.C. Foam::tmp<Foam::volScalarField> Foam::diameterModels::IATE::d() const { return d_; } If you want to calculate R[celli] in terms of kappai but not d then the you can replace d with 6/k which ensures both the units and the constant. Please see the equations 4 and 19 through 22 in the following link (https://snip.mathpix.com/nucerl/notes/11e456fc17f4427cab48b927c863844a) for detailed explanation. Thank you, Erol. 

fvm::Su(R, kappai) returns R as an explicit source which should and does have the units of kappai/time according to R[celli] = (Cti/18)*Ut[celli]*sqr(kappai[celli]) *sqrt(1  WeCr/We[celli])*exp(WeCr/We[celli]); From the references I have the term would be divided by 18. eq. 11163 in Mamoru Ishii, Takashi Hibiki ThermoFluid Dynamics of TwoPhase Flow Second Edition I cannot work out the correspondence between your notes and the references I have. 

Resolved in OpenFOAMdev by commit 3aea199ebebd440e0032ad40e61e6c968b248ba5 

It's still not correct. Regardless of the explicit/implicit treatment, the IATE source must return 1/ms. The partial differential equation (dk/dt) has a unit of 1/ms. As I already mentioned, the term having a factor of 1/18 does not mean we can directly use it. OpenFOAM does not solve for interfacial are concentration (ai), it solves for interfacial curvature (kappai). You can see this in the implementation of the sink terms as well; Random Collision (Rrc) and Wake Entrainment (Rwe). The sources in the uploaded image are the same as your reference but the implementation is different. For the Turbulent Breakup, there are 2 ways: 1) Using kappai: R[celli] = (Cti/18)*Ut[celli]*kappai[celli] *sqrt(1  WeCr/We[celli])*exp(WeCr/We[celli]); 2) Using d_: const volScalarField& d(iate_.d()()); R[celli] = (Cti/3)*Ut[celli]/d[celli] *sqrt(1  WeCr/We[celli])*exp(WeCr/We[celli]); 

The dimensions of (Cti/18)*Ut[celli]*kappai[celli] *sqrt(1  WeCr/We[celli])*exp(WeCr/We[celli]); are 1/s, but R has dimensions of 1/ms so your proposed expression is incorrect. 

I don't agree with your assessment of the 1/18. 

From my understanding the internal R should have the unit of 1/s and once it’s multiplied by kappai through Su() then source R will have the unit of 1/ms. This is how it is done for the Rrc source as well. Though it uses implicit source treatment, Sp(). So, unless Su() somehow multiplies R with kappai in a way that I don’t understand, then I don’t know why the defined unit of R[celli] should be 1/ms instead of 1/s. 

Where in Su() is R multiplied by kappai? 

I guess I have a misunderstanding of what return fvm::Su(R, kappai); does. I think that’s where my confusion comes from. 

Resolved by commit 3aea199ebebd440e0032ad40e61e6c968b248ba5 
Date Modified  Username  Field  Change 

20221127 12:12  nucerl  New Issue  
20221127 12:12  nucerl  File Added: Snap#1.png  
20221127 12:12  nucerl  File Added: Snap#2.png  
20221127 20:29  henry  Note Added: 0012886  
20221127 20:38  henry  Note Edited: 0012886  
20221127 23:36  nucerl  Note Added: 0012887  
20221128 07:28  henry  Note Added: 0012889  
20221128 08:08  henry  Note Edited: 0012889  
20221128 08:09  henry  Note Edited: 0012889  
20221128 08:12  henry  Note Edited: 0012889  
20221128 08:44  henry  Note Edited: 0012889  
20221128 08:48  henry  Note Added: 0012890  
20221128 12:40  nucerl  Note Added: 0012893  
20221128 12:40  nucerl  File Added: IATE General Form.png  
20221128 13:54  henry  Note Added: 0012894  
20221128 13:57  henry  Note Added: 0012896  
20221128 14:11  nucerl  Note Added: 0012897  
20221128 14:18  henry  Note Added: 0012899  
20221128 14:31  henry  Note Edited: 0012899  
20221128 14:34  nucerl  Note Added: 0012900  
20221128 14:58  henry  Assigned To  => henry 
20221128 14:58  henry  Status  new => resolved 
20221128 14:58  henry  Resolution  open => fixed 
20221128 14:58  henry  Fixed in Version  => dev 
20221128 14:58  henry  Note Added: 0012901 