CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Programming & Development

Modifying laplacianFoam

Register Blogs Members List Search Today's Posts Mark Forums Read

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   April 6, 2017, 06:18
Default Modifying laplacianFoam
  #1
Member
 
Lilian Chabannes
Join Date: Apr 2017
Posts: 58
Rep Power: 8
Lookid is on a distinguished road
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" (\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)

t_{i}=t_{0}exp(\frac{T{_{0}}^{}}{T}) (induction time)

And because the reaction is non-isotherm t_{induction}=\sum \frac{dt}{t_{i}}

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

\frac{\mathrm{d\alpha } }{\mathrm{d} t}= nk^{1/n}\alpha^\frac{n-1}{n}(1-\alpha)^\frac{n+1}{n} with k=k_{0}exp^{\frac{-E_{a}}{RT}}

or directly : \alpha =\frac{k(t-t_{induction})^{n}}{1+k(t-t_{induction})^{n}} (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 :

Code:
        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

Last edited by Lookid; April 11, 2017 at 10:48.
Lookid is offline   Reply With Quote

Old   April 7, 2017, 06:00
Default
  #2
Member
 
Lilian Chabannes
Join Date: Apr 2017
Posts: 58
Rep Power: 8
Lookid is on a distinguished road
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.
Lookid is offline   Reply With Quote

Old   April 7, 2017, 06:21
Default
  #3
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,930
Rep Power: 38
alexeym has a spectacular aura aboutalexeym has a spectacular aura about
Send a message via Skype™ to alexeym
Hi,

Code:
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?
alexeym is offline   Reply With Quote

Old   April 7, 2017, 06:26
Default
  #4
Member
 
Lilian Chabannes
Join Date: Apr 2017
Posts: 58
Rep Power: 8
Lookid is on a distinguished road
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
Lookid is offline   Reply With Quote

Old   April 7, 2017, 06:39
Default
  #5
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,930
Rep Power: 38
alexeym has a spectacular aura aboutalexeym has a spectacular aura about
Send a message via Skype™ to 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:

Code:
const dimensionedScalar ONE("1", a.dimensions(), 1.0);
... Foam::pow(a, (n - 1)/n)*Foam::pow((ONE - a), (n + 1)/n) ...
alexeym is offline   Reply With Quote

Old   April 7, 2017, 07:20
Default
  #6
Member
 
Lilian Chabannes
Join Date: Apr 2017
Posts: 58
Rep Power: 8
Lookid is on a distinguished road
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 :

Code:
            solve
            (

                fvm::ddt(a)
                - n*Foam::pow(k0*Foam::exp(-Ea/(R*T)),n)
        
            );
I have the same error :

Code:
--> 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
Lookid is offline   Reply With Quote

Old   April 7, 2017, 09:00
Default
  #7
Senior Member
 
Alexey Matveichev
Join Date: Aug 2011
Location: Nancy, France
Posts: 1,930
Rep Power: 38
alexeym has a spectacular aura aboutalexeym has a spectacular aura about
Send a message via Skype™ to alexeym
Quote:
Originally Posted by Lookid View Post
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.
alexeym is offline   Reply With Quote

Old   April 7, 2017, 10:38
Default
  #8
Member
 
Lilian Chabannes
Join Date: Apr 2017
Posts: 58
Rep Power: 8
Lookid is on a distinguished road
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.

Last edited by Lookid; April 7, 2017 at 12:03.
Lookid is offline   Reply With Quote

Old   April 10, 2017, 11:07
Default
  #9
Member
 
Lilian Chabannes
Join Date: Apr 2017
Posts: 58
Rep Power: 8
Lookid is on a distinguished road
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

Code:
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 is offline   Reply With Quote

Old   April 11, 2017, 10:49
Default
  #10
Member
 
Lilian Chabannes
Join Date: Apr 2017
Posts: 58
Rep Power: 8
Lookid is on a distinguished road
Problem solved , the lines of code added are in the first post.
Lookid is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Multi-region problem using laplacianFoam zfaraday OpenFOAM Pre-Processing 9 March 19, 2023 07:20
Heat source using fvOptions in laplacianFoam osha OpenFOAM Programming & Development 4 October 27, 2020 05:28
Add solar radiation to laplacianFoam aar_out OpenFOAM Programming & Development 2 April 30, 2017 16:12
laplacianFoam run-time error industrial0503 OpenFOAM Running, Solving & CFD 2 June 6, 2012 15:16
Serious error in calculating gradT (laplacianFoam) doubtsincfd OpenFOAM Bugs 1 December 30, 2011 03:56


All times are GMT -4. The time now is 21:40.