CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM (http://www.cfd-online.com/Forums/openfoam/)
-   -   Conjugate Mass Foam Solver with reaction? (http://www.cfd-online.com/Forums/openfoam/86433-conjugate-mass-foam-solver-reaction.html)

chegdan March 22, 2011 15:51

Conjugate Mass Foam Solver with reaction?
 
Hello Foamers,

I managed to combine one of my solvers that solves for time dependent passive scalar transport in a turbulent flow field for at particular species with mass transport in a solid similar to conjugateHeatFoam. This works without a hitch, but when I try to attached a simple first order reaction based on a stationary species in the solid phase it gives some errors on runtime. At first, i defined the separate species as its own field, not part of the coupledFvScalarMatrix. This gave an error:


--> FOAM FATAL ERROR:
incompatible fields for operation
[C] + [Cmetal]

I am sure it is from the source term (fvm::Sp(0.01*Csolid,Cmetal)) where i tried to couple a non coupleFvScalarMatrix field into the mix. I then decided to try and make a third region (identical to the solid region) that is overlayed on the solid region into a different mesh.
Code:

        // Add solid-side equation
        CEqns.set
        (
            1,
            new fvScalarMatrix
            (
                fvm::ddt(Csolid) - fvm::laplacian(Dsolid, Csolid) + fvm::Sp(0.01*Csolid,Cmetal)
            )
        );
        //solve the metal oxidation reaction only in the solid

        CEqns.set
        (
            2,
            new fvScalarMatrix
            (
                fvm::ddt(Cmetal) + fvm::Sp(0.01*Csolid,Cmetal)
            )
        );

I then recieved an error:

--> FOAM FATAL ERROR:
incompatible fields for operation
[C] + [C]

I'm assuming that the error was because I was trying to solve a system in which one variable was defined on one mesh....and another defined on a separate mesh. how does one overcome this issue?

Dan

p.s. I am using 1.6-ext

benk March 22, 2011 16:14

I'm pretty sure all the fields need to be on the same mesh, if they're not then you'll have to find a way of mapping the terms from one mesh to another. I would say try implementing your source term using an explicit source term first to see if that works, or even just a scalar value to see if it runs.

chegdan March 22, 2011 18:40

Quote:

Originally Posted by benk (Post 300604)
I'm pretty sure all the fields need to be on the same mesh, if they're not then you'll have to find a way of mapping the terms from one mesh to another. I would say try implementing your source term using an explicit source term first to see if that works, or even just a scalar value to see if it runs.

Mapping was not required. it turned out to be simple, but I still have to validate my results before I speak too soon. They could not be on the same mesh so I had to make another region over top the solid region with exactly the same mesh (just copied the polymesh folder). No fancy Sp, Su, or SuSp source handling was required. Only explicit definition of the source term was required. here is a snippet that I used:

Code:

        // Add solid-side equation
        CEqns.set
        (
            1,
            new fvScalarMatrix
            (
                fvm::ddt(Csolid) - fvm::laplacian(Dsolid, Csolid) + k*Csolid*Cmetal
            )
        );

        //add the solid-side equation for another scalar in the solid
        CEqns.set
        (
            2,
            new fvScalarMatrix
            (
                fvm::ddt(Cmetal) + k*Csolid*Cmetal
            )
        );

        CEqns.solve();

Thanks for the suggestion of using just an explicit source term. That was the sticking point.

Dan

Cyp March 23, 2011 05:16

Hi!

If you want to use fvm::Sp() in the building of your matrix, you ought to inverse Csolid and Cmetal within the fonction Sp(). Indeed the second argument is the variable of the matrix:

Code:

fvm::Sp(0.01*Cmetal,Csolid)
The solution you mentioned (that is to say to build the matrix with k*Csolid*Cmetal as a source term which is not included in the matrix also work. However, I think the use of fvm::Sp() is advised when it is possible to use it.

Bests,
Cyp

chegdan March 23, 2011 10:09

Quote:

Originally Posted by Cyp (Post 300678)
Hi!

If you want to use fvm::Sp() in the building of your matrix, you ought to inverse Csolid and Cmetal within the fonction Sp(). Indeed the second argument is the variable of the matrix:

Code:

fvm::Sp(0.01*Cmetal,Csolid)
The solution you mentioned (that is to say to build the matrix with k*Csolid*Cmetal as a source term which is not included in the matrix also work. However, I think the use of fvm::Sp() is advised when it is possible to use it.

Bests,
Cyp

The programmer's guide states that stating a source term as (k*Cmetal*Csolid) is fully implicit (http://www.foamcfd.org/Nabla/guides/...x14-410002.4.9) and I'm not clear on the advantage of fvm::Sp() over k*Cmetal*Csolid (would be great to know is you have an idea). However, I received an error at runtime (but compiles correctly):

Code:

--> FOAM FATAL ERROR:
incompatible fields for operation
    [C] + [C]

    From function checkMethod(const fvMatrix<Type>&, const fvMatrix<Type>&)
    in file /home/dcombest/OpenFOAM/OpenFOAM-1.6-ext/src/finiteVolume/lnInclude/fvMatrix.C at line 1207.

FOAM aborting

when using this code snippet:

Code:

// Add solid-side equation
        CEqns.set
        (
            1,
            new fvScalarMatrix
            (
                fvm::ddt(Csolid) - fvm::laplacian(Dsolid, Csolid) + fvm::Sp(k*Cmetal,Csolid)//k*Csolid*Cmetal
            )
        );

        //add the solid-side equation for another scalar in the solid
        CEqns.set
        (
            2,
            new fvScalarMatrix
            (
                fvm::ddt(Cmetal) + fvm::Sp(k*Cmetal,Csolid)//k*Csolid*Cmetal
            )
        );

        CEqns.solve();

however, the last code snippet:

Code:

        // Add solid-side equation
        CEqns.set
        (
            1,
            new fvScalarMatrix
            (
                fvm::ddt(Csolid) - fvm::laplacian(Dsolid, Csolid) + k*Csolid*Cmetal
            )
        );

        //add the solid-side equation for another scalar in the solid
        CEqns.set
        (
            2,
            new fvScalarMatrix
            (
                fvm::ddt(Cmetal) + k*Csolid*Cmetal
            )
        );

        CEqns.solve();

compiles and executes at runtime without errors. If I substitute a simple 0.1 scalar instead of my dimensioned scalar k (read in at runtime from transport properties) then I receive an error about dimensions not matching...which is expected. I will do some more test today and make notes here if there is anything wrong. The help is much appreciated!

Dan

Cyp March 23, 2011 10:18

Hi!

Look at your second equation : if you want to construct a matrix for the variable Cmetal, the second argument of your source term fvm::Sp() must be Cmetal :

Code:

fvm::Sp(k*Csolid,Cmetal)

I think it should work now.

chegdan March 23, 2011 11:44

Excellent!
 
That did the trick! I'm still unclear as to why the fvm::Sp() is better than just an explicit k*C*C in the equation. Is this a situation where k*C*C is collected in the b of Ax=b and fvm::Sp() collects the term into A of Ax = b? Is there an advantage? Thanks again for the help.

Dan

chegdan March 23, 2011 11:51

I actually found the answer on the message boards here (http://www.cfd-online.com/Forums/ope...-fvm-sp-b.html) which states:

Quote:


fvm::Sp(A, B) is used for an implicit source term, so "B" is the variable you are solving for. Using A*B will leave the term on the right hand side of the equations. Depending on the sign of your source term implicit treatment can help with stability.


Cyp March 23, 2011 11:53

Quote:

Is this a situation where k*C*C is collected in the b of Ax=b and fvm::Sp() collects the term into A of Ax = b?
Yes indeed, it is exactly the point! Using fvm::Sp() makes your resolution fully implicit. And most of the time the resolution is better (= better stablility).


All times are GMT -4. The time now is 05:54.