[KingKraut]

Hi all,

I am currently trying to implement a new transport model describing non-newtonian behaviour according to the formula (attached file BloodViscosity.PNG) with the functions (attached file BloodViscosity_Functions.PNG).

BloodViscosity.PNG

BloodViscosity_Functions.PNG

I created the file structure for the transportModel as described in http://www.tfd.chalmers.se/~hani/kur...lenceModel.pdf and took care of the division by the strainRate() similar to the implementation as it is done in the crossPowerLaw transportModel, which comes with the OpenFOAM compilation.

Foam::tmp<Foam::volScalarField>
Foam::viscosityModels::generalizedPowerLaw::calcNu () const
{
return
(nuInf_ + deltaNu_*exp(-(scalar(1) + strainRate()/a_)*exp(-b_/max(strainRate(), 1/dimensionedScalar("VGREAT", dimTime, VGREAT) ) ) ) )
* pow(max(strainRate()/one_, dimensionedScalar("VSMALL", dimless, VSMALL)),
(nInf_ - deltaN_*exp(-(scalar(1) + strainRate()/c_)*exp(-d_/max(strainRate(), 1/dimensionedScalar("VGREAT", dimTime, VGREAT) ) ) ) ) - 1.0);

}

The library is compiled fine with wmake libso and the inclusion in the controlDict works fine, too. If I then start a simulation run, the calculations work fine until after a couple of timesteps the core is dumped. Either the bounding k explodes or the time step continuity error becomes too large.
This happens with a custom solver (based on pisoFoam) as well as with pisoFoam itself.
I guess the implementation of the formula is somehow bogus, but I don't know how I could do this more elegantly.
Or possibly the boundary conditions for k and nuSgs are chosen wrong? However, this seems weird to me, because the same BCs work fine in the Newtonian case.
Can anybody give me a hint to a direction, where I should be looking for an improvement? Or maybe someone sees a clear mistake in the code above, which I am missing?

Thanks a lot to anyone looking into this!

Best regards
JOhannes