|
[Sponsors] |
August 21, 2014, 13:04 |
Adding porosity to UEqn.H in porousInterFoam
|
#1 |
Member
james wilson
Join Date: Aug 2014
Location: Orlando, Fl
Posts: 39
Rep Power: 12 |
Hi all,
I would like to modify UEqn.H in porousInterFoam to incorporate porosity in the solver. I defined porosity by linking it to a zone called "porous" in blockMeshDict and set the field using setFieldsDict and zoneToCell. I have defined porosity and porosityf in createFields.H as follows: /////////////////createFields.H//////////////////////////////// ... // Read porosity from /0 volScalarField porosity ( IOobject ( "porosity", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE ), mesh ); // Converting porosity to surfaceScalarField for operation compatability in UEqn. surfaceScalarField porosityf ( IOobject ( "porosityf", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), fvc::interpolate(porosity) ); ... /////////////////createFields.H//////////////////////////////// and now in UEqn.H I update the values every time step: /////////////////UEqn.H//////////////////////////////// // Calculate and cache mu for the porous media volScalarField mu(twoPhaseProperties.mu()); volScalarField rho_porosity ( "rho_porosity", rho/porosity ); surfaceScalarField rhoPhi_porosityf ( "rhoPhi_porosityf", rhoPhi/(porosityf*porosityf) ); surfaceScalarField muEff_porosityf ( "muEff_porosityf", muEff/porosityf ); fvVectorMatrix UEqn ( //pZones.ddt(rho, U) fvm::ddt(rho_porosity, U) + fvm::div(rhoPhi_porosityf, U) - fvm::laplacian(muEff_porosityf, U) - (fvc::grad(U) & fvc::grad(muEff)) //- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U)) ) & mesh.Sf())) == fvOptions(rho, U) ); ... /////////////////UEqn.H//////////////////////////////// This compiles fine and the solver runs but diverges shortly into the simulation as I am running a case identical (with the exception of the additional porosity terms) to: http://www.cfdandco.com/resources/JR...sInter_v11.pdf where the liquid is adjacent to the interface at t=0. Without the porosity terms (unmodified solver), this case runs fine. with the porosity terms (1 in free flow and 0.32 in porous zone) the solution diverges. My attempts to debug have led me to some good information that I do not know what to do with : ) When I make the porosity field uniform (both the volScalarField and the surfaceScalarField for porosity) with a value of 1 (free flow), the solution still diverges. This case is essentially dividing the modified terms by 1 and should have no effect on the solution, which is strange that it still diverges. something/1 = something, but the divergence here seems like something/1 = something_else. my example is crude and probably doesn't represent whats actually happening. Clearly there is something happening as a result of my modifications that I am too novice to realize. I have experienced success in eliminating diverging temperature fields in an unrelated case by streamlining the arguments of operators such as fvm::div(rhoCpPhi, T) by shifting to rhoCp*fvm::div(phi, T). This is also something I'd like to implement here with UEqn.H for example: fvVectorMatrix UEqn ( //pZones.ddt(rho, U) fvm::ddt(rho, U) / porosity + fvm::div(rhoPhi, U) /(porosityf*porosityf) - fvm::laplacian(muEff, U) / porosityf - (fvc::grad(U) & fvc::grad(muEff)) //- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U)) ) & mesh.Sf())) == fvOptions(rho, U) ); But i get a compile error when I try to compile porousInterFoam with the only modification being "fvm::ddt(rho, U) / porosity": UEqn.H:33:28: error: no match for ‘operator/’ in ‘Foam::fvm::ddt(const volScalarField&, const Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh>&) [with Type = Foam::Vector<double>, Foam::volScalarField = Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>]((*(const Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh>*)(& U))) / porosity’ I would appreciate any insight regarding operations within fvVectorMatrix and the effects of adding arguments to something like fvm::div() vs. moving them outside of the divergence operator. Also I am curious as to why I get the compile error when operating on this term fvm::div()*something/something_else and why I don't when I define the field first and then pass it to the fvm::div() operator. My over arching concern is naturally the divergence of the solution, but also, why dividing by a uniform field of 1 causes this divergence. Thanks in advance, James PS I'm new to the forums here and have benefited immensely from simply reading related problems! Thanks all |
|
August 21, 2014, 22:33 |
|
#2 |
Senior Member
Kyle Mooney
Join Date: Jul 2009
Location: San Francisco, CA USA
Posts: 323
Rep Power: 18 |
Hi James,
You might be better off replicating existing implementations of porosity models. For example, in foam 2.3 there is an extra term (fvOptions(U)) added onto the momentum equation in many of the solvers like this (from pimpleFoam): Code:
tmp<fvVectorMatrix> UEqn ( fvm::ddt(U) + fvm::div(phi, U) + turbulence->divDevReff(U) == fvOptions(U) ); Code:
porosity1 { type explicitPorositySource; active yes; selectionMode cellZone; cellZone stator; explicitPorositySourceCoeffs { type DarcyForchheimer; DarcyForchheimerCoeffs { d d [0 -2 0 0 0 0 0] (1e5 1000 1000); f f [0 -1 0 0 0 0 0] (0 0 0); coordinateSystem { type cartesian; origin (0 0 0); coordinateRotation { type axesRotation; e1 (1 0 0); e2 (0 1 0); } } } } } I hope that helps you some! Cheers, Kyle |
|
August 22, 2014, 22:17 |
|
#3 |
Member
james wilson
Join Date: Aug 2014
Location: Orlando, Fl
Posts: 39
Rep Power: 12 |
Thanks Kyle,
This is some great info! I looked into some models that incorporate porosity and they seem to modify the transient term only. Im set on modifying UEqn.H as the terms i need to modify are in this header. See the attached photo. Pi is already a feature in porousInterFoam but porosity is not. I would think, just as I have built custom functions and equations for something like the heat equation, I could just modify the terms in UEqn.H but the solver doesnt seem to like the modifications, even if it involves dividing by 1. Ill keep looking and continue to explore fvoptions as well. Thanks again! Any more suggestions? James |
|
August 23, 2014, 04:35 |
|
#4 |
Senior Member
Niels Gjoel Jacobsen
Join Date: Mar 2009
Location: Copenhagen, Denmark
Posts: 1,903
Rep Power: 37 |
Hallo James,
You might want to look into the waves2Foam package. It comes with a porous solver, which looks a lot like the set of equations sketched above. One of the tutorials is even identical with the one you are trying (though other dimensions), but it comes with experimental data and a post-processing script for the validation. Kind regards, Niels
__________________
Please note that I do not use the Friend-feature, so do not be offended, if I do not accept a request. |
|
August 23, 2014, 10:07 |
|
#5 |
Member
james wilson
Join Date: Aug 2014
Location: Orlando, Fl
Posts: 39
Rep Power: 12 |
Thanks Niels, I'll take a look at that asap.
James |
|
November 25, 2016, 23:19 |
|
#6 | |
Member
Anirudh Kulkarni
Join Date: May 2016
Posts: 62
Rep Power: 10 |
Quote:
|
||
Tags |
divergence, fvvectormatrix, porosity, porousinterfoam, solver compilation error |
Thread Tools | Search this Thread |
Display Modes | |
|
|
Similar Threads | ||||
Thread | Thread Starter | Forum | Replies | Last Post |
chtMultiRegionSimpleFoam | samiam1000 | OpenFOAM Running, Solving & CFD | 39 | March 31, 2016 09:43 |
Porosity and permeability in PorousInterFoam | jasouza1974 | OpenFOAM Pre-Processing | 14 | August 12, 2015 05:16 |
help adding the energy equation to porousinterfoam using the enthalpy formulation | nadine | OpenFOAM Programming & Development | 18 | June 17, 2014 09:39 |
Help with chtMultiRegionFoam | jbvw96 | OpenFOAM Running, Solving & CFD | 2 | December 26, 2010 18:16 |
Porosity field in Fluent | wojciech | FLUENT | 1 | September 20, 2010 12:19 |