[Lookid]

Hello everyone,

[Editing this post when some things are solved]

I'm totally new to OpenFOAM and i'm trying to modify laplacianFoam to calculate a parameter "alpha" () (cure degree for a chemical reaction). I'm trying to do simple things to begin but i'm lost .

So basically, I want to calculte this first (it calculates an induction time, because the reaction doesn't start immediatly)

(induction time)

And because the reaction is non-isotherm

And when this sum equals 1 (induction is over), the reaction begins :

with

or directly : (cure rate of the reaction)

Everything is constant except T (temperature) and t (time)

- t0 [s] : time constant
- T0 [K] : Temperature constant
- k0 [s-1] : a frequency factor
- Ea [J/mol] = [kg m2 s-2 mol-1] : activation energy
- n [-] : reaction degree

Thanks to alexeym, everything is solved, here are the lines of codes added to the solver :

        volScalarField ti(t0*Foam::exp(T0/T));             //calculation of isothermal induction time
        sumti += runTime.deltaT()/ti;                //calculation to know when we reach the induction time

        volScalarField xi(pos(sumti - 1));            //xi = 0 before induction reached for each cells, and 1 and it is reached
        forAll(xi, i)
        {
            // tsum[i] > 1.0 for the first time
            if (xi.v()[i] == 1.0 and tinduction.v()[i] == 0.0)    //if induction time is over AND tinduction has not been modified yet
            {
                tinduction[i] = runTime.value();        //assign the value of the induction in the volScalarField tinduction
            }
        }

        volScalarField k(k0*Foam::exp(-Ea/R/T));        //declaration and calculation of k, to help in the calcul of alpha


        Info<< "\nCalculating cure rate distribution\n" << endl;

        // For a cell, if induction time is not over, alpha = 0 (due to xi)
        //            if induction time is over, alpha is calculated
        alpha =
            xi*k*Foam::pow(runTime.value()-tinduction, n)/(k*Foam::pow(runTime.value()-tinduction, n)+1/k1);
    //the term k1 in the end is here to avoid dimensional problems, it doesn't change the result

[Lookid]

I'm trying to implement da/dt, but I have a problem with pow function.

Add this on laplacianFoam

solve
(
fvm::ddt(a)
- n*Foam:ow(k0*Foam::exp(-Ea/(R*T)),1/n)*Foam:ow(a,(n-1)/n)*Foam:ow((1-a),(n+1)/n)
);

It compiles well, but when I run it I have this error :

[a[0 0 -1 0 0 0 0] ] - [(((n*pow((k0*exp((-Ea|(R*T)))),(1|n)))*pow(a,((n-1)|n)))*pow((1-a),((n+1)|n)))[0 0 -0.465116 0 0 0 0] ]

From function void Foam::checkMethod(const Foam::fvMatrix<Type>&, const Foam:imensionedField<Type, Foam::volMesh>&, const char*) [with Type = double]
in file /home/lilian/OpenFOAM/OpenFOAM-4.1/src/finiteVolume/lnInclude/fvMatrix.C at line 1292.

FOAM aborting

#0 Foam::error:rintStack(Foam::Ostream&) at ??:?
#1 Foam::error::abort() at ??:?
#2 void Foam::checkMethod<double>(Foam::fvMatrix<double> const&, Foam:imensionedField<double, Foam::volMesh> const&, char const*) at ??:?
#3 ? at ??:?
#4 ? at ??:?
#5 __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#6 ? at ??:?
Abandon (core dumped)


This term pow((k0*exp((-Ea|(R*T)))),(1|n))) should be of dimension [0 0 -1 0 0 0 0].
k0 is like this, Ea/RT is dimensionless, and n is also dimensionless. But on the error : [0 0 -0.465116 0 0 0 0], 1/n = 0.465116, so it multiplies the value of n with the dimension of k0, I don't get it. If I put 1 on the pow functions, then it works.

[alexeym]

Hi,

Foam::pow(a,(n-1)/n)*Foam::pow((1-a),(n+1)/n)

1 in (1 - a) in fact should be dimensioned and has dimensions of a.

Why do you want to solve equation instead of direct calculation of a?

[Lookid]

Because in the initial formula, I have a term "t^n", and I don't know how to implement t.
Plus, I may add a source term in the energy equation later, which is proportional to da/dt.

Thanks for the tip, but a problem comes from before, when I try to solve this :

fvm::ddt(a) - n*pow((k0*exp((-Ea|(R*T)))),(1|n)))

I already have the dimension problem

[alexeym]

As I said, your dimension problem comes from the fact, that you have used just 1 instead of dimensioned 1. So 1 - a is in fact dimensionless instead of 1/s (since instead of complaining modern OpenFOAM will construct dimless dimensionedScalar from 1 and then during - operation it will again take dimensions of the first operand).

Use something like:

const dimensionedScalar ONE("1", a.dimensions(), 1.0);
... Foam::pow(a, (n - 1)/n)*Foam::pow((ONE - a), (n + 1)/n) ...

[Lookid]

Thanks again but it's still not working.
I tried stupid thing like declare "one" in transportproperties, and still the same.

Even when I try without the "1", like this :

            solve
            (

                fvm::ddt(a)
                - n*Foam::pow(k0*Foam::exp(-Ea/(R*T)),n)
        
            );

I have the same error :

--> FOAM FATAL ERROR: 

    [a[0 0 -1 0 0 0 0] ] - [(n*pow((k0*exp((-Ea|(R*T)))),n))[0 0 -2.15 0 0 0 0] ]

    From function void Foam::checkMethod(const Foam::fvMatrix<Type>&, const Foam::DimensionedField<Type, Foam::volMesh>&, const char*) [with Type = double]
    in file /home/lilian/OpenFOAM/OpenFOAM-4.1/src/finiteVolume/lnInclude/fvMatrix.C at line 1292.

It takes the dimension of k0, and multiply it by n

[alexeym]

Originally Posted by Lookid

It takes the dimension of k0, and multiply it by n

And what is wrong with it? You write pow(k0*exp(...), n); exp(...) is dimensionless, so dimensions of the expression are dimensions of k0 in n-th power. In terms of SI base unit powers it is multiplication by n.

[Lookid]

Yes of course, sorry for this .
I put another coefficient in front of it of dimension 1/n and value 1, but that's not how it should work.
With this both sides have the same dimensions but the result is nonsense, 'a' should vary from 0 to 1 (1 = reaction over), but it stays 0.
Maybe there's a way to say "solve this just with values and don't look at the dimensions".

I modified the first post to explain my whole problem clearly.

[Lookid]

Thanks to alexeym, i'm doing some steps in my problem ! I edited the first post again to show what has been done.

the calculation of alpha is made, and alpha = 0 when induction is not over, which is good. But the calculation of alpha starts at the current time and it has to start at t=0, which means implement t = current time - induction time for each nodes but I didn't success in it yet.
A volScalarField t_induction is created, and I'm trying to store the time value when injection is over for each nodes with something like

if (xi=1) {
          t_induction = runTime.deltaT()...

But I already have an error because xi (volScalarField) = 1 is not possible.

Do you have any clues ?

[Lookid]

Problem solved , the lines of code added are in the first post.