[Lee Kyun Bum]

Happy new year, openfoamers


Since the last November, I've been learning openFoam and at the same time, I have been developing a new viscosity code with openFoam7. It's about temperature-dependent viscosity code based on Casson viscosity model.

before doing that, I already implemented a temperature equation to simpleFoam solver and did some test for Newtonian viscosity case. That was successful.

I referred to the PDF file:
http://www.tfd.chalmers.se/~hani/kur...nFoam%20v2.pdf

Then, I tested the above link's temperature-dependent powerLaw model and that also worked well with my solver.



However, after implementing my own temperature-dependent viscosity model to my simpleFoam solver, my solver stopped without any message, as follows:
__________________________________________________ ___________
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create mesh for time = 0


SIMPLE: Convergence criteria found
p: tolerance 1e-07
U: tolerance 1e-07
T: tolerance 1e-07

Reading transportProperties

Reading field p

Reading field U

Reading field T

Reading/calculating face flux field phi

Selecting incompressible transport model CassonTemp
__________________________________________________ ___________

The above message log is only thing I saw when my solver stopped.

Anyway, I noticed that there might be some problem in my code, especially the part of calculation of viscosity nu in my viscosity source file(CassonTemp.C).

my viscosity model written in the source file:
__________________________________________________ ________________________________

// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //

Foam::tmp<Foam::volScalarField>
Foam::viscosityModels::CassonTemp::calcNu() const
{
const volScalarField& T_ = U_.mesh().lookupObject<volScalarField>("T"); // T_ is defined!

//For Casson yield stress tau0 as tau0(T)=a1/(T+b1)+c1
dimensionedScalar a1_("a1",dimTemperature, 3.924);
dimensionedScalar b1_("b1",dimTemperature, -306.4);
dimensionedScalar c1_("c1",dimless, 7.067);
//For Casson plastic viscosity m as m(T)=a2/(T+b2)+c2
dimensionedScalar a2_("a2",dimTemperature, 39.42);
dimensionedScalar b2_("b2",dimTemperature, -279.3);
dimensionedScalar c2_("c2",dimless, -0.3417);

//Calculation of tau0 field
volScalarField tau0_ =
a1_*pow(T_+b1_,-1.0)+c1_)*dimensionedScalar("one",dimViscosity/dimTime,1.0);

//Calculation of m field
volScalarField m_ = (a2_*pow(T_+b2_,-1.0)+c2_)*dimensionedScalar("one",dimViscosity,1.0 );

//Calculation of nu1_ field
volScalarField nu1_ = sqr(sqrt(tau0_/max(strainRate(),dimensionedScalar("VSMALL",dimles s/dimTime,VSMALL)))+sqrt(m_*pow(max(strainRate(),dim ensionedScalar("VSMALL",dimless/dimTime,VSMALL)),-1.0)));

return max(nuMin_,min(nuMax_,nu1_/rho_.value()));
}
__________________________________________________ ________________________________

My idea was, using some measured data of two parameters tau0 and m at different temperatures, I did some nonlinear regression and approximated the parameters, tau0 and m as functions of temperature, i.e tau0(T) and m(T). After calculating volScalarfield with the two functions of temperatures, I wanted to calculate the viscosity nu.

Originally, the apparent viscosity nu of Casson model is defined as:
nu = sqr[sqrt(tau0/strainRate)+sqrt(m)] where tau0: yield stress and m: Casson viscosity

-> For my case, using the mentioned functions of temperature, tau0(T) and m(T), the above formula would be rewritten as:
nu = sqr[sqrt( tau0(T) / strainRate ) + sqrt(m(T) )]


However, I was wrong.

Since I spent too much time for finding my incorrect parts and trying some modifications, I desperately ask you for advice and help.
I'm still confusing because some operation are possible between dimensioned scalar and scalar field.



Additionally, I have some questions,

1.Is there any way to define a new volScalarfield from an elementwise multiplication or division of a volScalarfield tau0 by different volScalarfield or strainRate()?

What I want to do is, i.e.
dimensionedScalar a, b, c, d, e, f, g (...);
volScalarField temperature = (...);
volScalarField viscosity = [(a*(temperature-b)+c)/strainRate() + (d*(temperature-f)+g)]^2;



2. If I define any volScalarfield, do they have the same 'celli'?
ex) volScalarField& T_ = (...);
volScalarField& shr_ = strainRate();

T_[celli] == shr_[celli]? or T_[celli] != shr_[celli]?

The reason why I'm asking about this is, as far as I know any volScalarfields have their 'innerfield' at the cell center and 'boundary field' defined on the boundary.



3.Can exponential function, exp(volScalarfield a) be a new volScalarfield? If so, they have the same dimension?

Thank you in advance!