|
[Sponsors] |
February 6, 2019, 05:15 |
Modified laplacian solver
|
#1 |
Member
Join Date: Oct 2018
Location: France
Posts: 34
Rep Power: 7 |
Dear all,
I tried to modify laplacianFoam to - be able to use temperature dependent material properties and - obtain the generated heat flux at the boundaries. I'm sure there already exist many solutions, function objects etc. to get this done, but for exercise I wanted to try that by my own. My code compiles and the solver in principle works. But during coding some questiosn came up I wanted to ask. But to start with the codes of the three modified files: laplacianFoam.C - modified Code:
#include "fvCFD.H" #include "fvOptions.H" #include "simpleControl.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // int main(int argc, char *argv[]) { #include "setRootCaseLists.H" #include "createTime.H" #include "createMesh.H" simpleControl simple(mesh); bool varProps = false; #include "createFields.H" //First try: code snippet example to read information from controlDict - not working if (runTime.controlDict().lookup("variable") == "yes" || runTime.controlDict().lookup("variable") == "true") { varProps = true; } // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Info<< "\nCalculating temperature distribution\n" << endl; while (simple.loop(runTime)) { Info<< "Time = " << runTime.timeName() << nl << endl; while (simple.correctNonOrthogonal()) { if(varProps == true) //variable material properties { volScalarField k = c0_k * pow(T,0) + c1_k * pow(T,1) + c2_k * pow(T,2) + c3_k * pow(T,3) + c4_k * pow(T,4) + c5_k * pow(T,5) + c6_k * pow(T,6) + c7_k * pow(T,7); volScalarField c = c0_c * pow(T,0) + c1_c * pow(T,1) + c2_c * pow(T,2) + c3_c * pow(T,3) + c4_c * pow(T,4) + c5_c * pow(T,5) + c6_c * pow(T,6) + c7_c * pow(T,7); volScalarField rho = c0_rho * pow(T,0) + c1_rho * pow(T,1) + c2_rho * pow(T,2) + c3_rho * pow(T,3) + c4_rho * pow(T,4) + c5_rho * pow(T,5) + c6_rho * pow(T,6) + c7_rho * pow(T,7); volScalarField thermDiff = k / (c * rho); fvScalarMatrix TEqn ( fvm::ddt(T) - fvm::laplacian(thermDiff, T) == fvOptions(T) ); fvOptions.constrain(TEqn); TEqn.solve(); fvOptions.correct(T); } else //constant material properties { fvScalarMatrix TEqn ( fvm::ddt(T) - fvm::laplacian(DT, T) == fvOptions(T) ); fvOptions.constrain(TEqn); TEqn.solve(); fvOptions.correct(T); } } #include "write.H" Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" << " ClockTime = " << runTime.elapsedClockTime() << " s" << nl << endl; } Info<< "End\n" << endl; return 0; } Code:
Info<< "Reading field T\n" << endl; volScalarField T ( IOobject ( "T", runTime.timeName(), mesh, IOobject::MUST_READ, IOobject::AUTO_WRITE ), mesh ); Info<< "Reading transportProperties\n" << endl; IOdictionary transportProperties ( IOobject ( "transportProperties", runTime.constant(), mesh, IOobject::MUST_READ_IF_MODIFIED, IOobject::NO_WRITE ) ); //Second try: code snippet example to read information from transportDict - not working if(transportProperties.lookup("variable")=="yes" || transportProperties.lookup("variable")=="true") { varProps = true; } Info<< "Determining material properties" << nl << endl; Info<< "Material properties variable (phi = f(T))" << " Reading polynomial coefficients for thermal conductivity" << endl; dimensionedScalar c0_k (transportProperties.lookup("c0_k")); dimensionedScalar c1_k (transportProperties.lookup("c1_k")); dimensionedScalar c2_k (transportProperties.lookup("c2_k")); dimensionedScalar c3_k (transportProperties.lookup("c3_k")); dimensionedScalar c4_k (transportProperties.lookup("c4_k")); dimensionedScalar c5_k (transportProperties.lookup("c5_k")); dimensionedScalar c6_k (transportProperties.lookup("c6_k")); dimensionedScalar c7_k (transportProperties.lookup("c7_k")); Info<< " Reading polynomial coefficients for specific heat capacity" << endl; dimensionedScalar c0_c (transportProperties.lookup("c0_c")); dimensionedScalar c1_c (transportProperties.lookup("c1_c")); dimensionedScalar c2_c (transportProperties.lookup("c2_c")); dimensionedScalar c3_c (transportProperties.lookup("c3_c")); dimensionedScalar c4_c (transportProperties.lookup("c4_c")); dimensionedScalar c5_c (transportProperties.lookup("c5_c")); dimensionedScalar c6_c (transportProperties.lookup("c6_c")); dimensionedScalar c7_c (transportProperties.lookup("c7_c")); Info<< " Reading polynomial coefficients for density" << nl << endl; dimensionedScalar c0_rho (transportProperties.lookup("c0_rho")); dimensionedScalar c1_rho (transportProperties.lookup("c1_rho")); dimensionedScalar c2_rho (transportProperties.lookup("c2_rho")); dimensionedScalar c3_rho (transportProperties.lookup("c3_rho")); dimensionedScalar c4_rho (transportProperties.lookup("c4_rho")); dimensionedScalar c5_rho (transportProperties.lookup("c5_rho")); dimensionedScalar c6_rho (transportProperties.lookup("c6_rho")); dimensionedScalar c7_rho (transportProperties.lookup("c7_rho")); Info<< "Thermal properties constant" << nl << " Reading thermal diffusivity DT\n" << endl; dimensionedScalar DT (transportProperties.lookup("DT")); #include "createFvOptions.H" Code:
if (runTime.writeTime()) { volVectorField gradT(fvc::grad(T)); volScalarField gradTx ( IOobject ( "gradTx", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), gradT.component(vector::X) ); volScalarField gradTy ( IOobject ( "gradTy", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), gradT.component(vector::Y) ); volScalarField gradTz ( IOobject ( "gradTz", runTime.timeName(), mesh, IOobject::NO_READ, IOobject::AUTO_WRITE ), gradT.component(vector::Z) ); Info << nl << endl; //CALCULATE AND SAVE BOUNDARY HEAT FLUXES: forAll(mesh.boundary(), iPatch) { //CELL CENTERS: //list containing the indices of the boundary neighbour cells (patch cells) labelList neighbourCells = mesh.boundary()[iPatch].faceCells(); //scalar field containing the inverse values of the cell center-face distances of the patch cells scalarField dCentreCell = mesh.boundary()[iPatch].deltaCoeffs(); //scalar field containing the temperatures of the cell centers scalarField cellTemperatures; cellTemperatures.resize(T.boundaryField()[iPatch].size()); //scalar field containing the thermal conductivities of the patch cells scalarField thermConductivity; thermConductivity.resize(T.boundaryField()[iPatch].size()); //loop over patch cells to obtain cell centre temperatures forAll(mesh.boundary()[iPatch], iFace) { label iCell = neighbourCells[iFace]; //index of patch cell cellTemperatures[iFace] = T.internalField()[iCell]; //temperature of patch cell scalar TC = cellTemperatures[iFace]; //temperature of patch cell //scalar TC = T.internalField()[iCell]; //should be identical with the two code lines above - does not work thermConductivity[iFace] = c0_k.value() * pow(TC ,0) + c1_k.value() * pow(TC ,1) + c2_k.value() * pow(TC ,2) + c3_k.value() * pow(TC ,3) + c4_k.value() * pow(TC ,4) + c5_k.value() * pow(TC ,5) + c6_k.value() * pow(TC ,6) + c7_k.value() * pow(TC ,7); //th. conductivity of patch cell } //CELL FACES: //scalar field containing the areas of the faces on the patch scalarField faceAreas = mesh.boundary()[iPatch].magSf(); //scalar field containing the temperatures of the cell faces scalarField faceTemperatures = T.boundaryField()[iPatch]; //BOUNDARY HEAT FLUX: //scalar field containing the face normals of the face cells scalarField patchHeatFlux = thermConductivity * (faceTemperatures - cellTemperatures) * dCentreCell * faceAreas; //overall heat flux over patch scalar Q = gSum(patchHeatFlux); Info<< "Total heat flux passing patch " << mesh.boundary()[iPatch].name() << " = " << Q << endl; //GENERATE IO-OBJECT WITH BOUNDARY HEAT FLUX: scalarField boundaryHeatFlux ( IOobject ( mesh.boundary()[iPatch].name(), runTime.timeName(), T.boundaryField()[iPatch].mesh, //problem argument - what to put here to make that work? IOobject::NO_READ, IOobject::AUTO_WRITE ), patchHeatFlux ); } Info << nl << endl; runTime.write(); } Now to ask some specific questions: 1) How do I implement and read a new boolean statement in the controlDict and/or the transportProperties? I tried both and my trials are coloured in red in the code snippets (laplacian.C & createFields.H). Unfortunately none of them work. 2) I would like to create an output file with a list of the heat fluxes of all boundary cells for each patch. Unfortunately this partt of the code doesn't compile, the problem is the registryObject and I don't get what I have to put there to make it work (coloured in blue in write.H). Could somebody briefly explain what I have to put in there? 3) For me the two upper lines to obtain the cell center temperatures coloured in green (write.H) are identical with the single (commented) line below. However, if I use the single line command, the results differ (and are obviously wrong). Does anyone know why that is? Thank you very much in advance, stockzahn Last edited by stockzahn; February 6, 2019 at 08:14. |
|
February 9, 2019, 10:37 |
|
#2 | |
Senior Member
|
Hi stockzahn,
Quote:
If you use a bool type, - you can use only “true” or “false”, not “yes|no|on|off”, - there is no need to use an equal-to operator (==) because the bool variable itself will be the answer (true or false). The function of a dictionary class returns ITstream, not bool. https://github.com/OpenFOAM/OpenFOAM...tionary.H#L320 It is impossible - to compare different class and type, and - to assign a type to another class. If you prefer to use a bool type, lookupOrDefault function with a bool argument returns bool. (what you give will be returned as a template class) The following line will work as what you want (I did NOT tested it , sorry…) Code:
bool varProps = runTime.controlDict().lookupOrDefault(”variable", false) https://github.com/OpenFOAM/OpenFOAM...witch/Switch.H “runTimeModifiable” entry in controlDict is declared as Swith. “adjustTimeStep” entry is dealt as Swith by explicitly specifying it. Code:
lookupOrDefault<Switch>("adjustTimeStep", false) The code of Time class will be a good example of dictionary handling. |
||
February 13, 2019, 10:59 |
|
#3 | |||
Member
Join Date: Oct 2018
Location: France
Posts: 34
Rep Power: 7 |
Thanks for your response snak,
unfortunately I don't quite understand your explanation. I tried to get the code in the links you send me, but I wasn't able to correct my mistake. Quote:
Thanks for these clarifications, but to which equal-to operator (==) you are referring to? In my code Code:
if (runTime.controlDict().lookup("variable") == "yes" || runTime.controlDict().lookup("variable") == "true") { varProps = true; } Quote:
So the comparison of the return value of "lookup"-function with a word/string ("yes", "true") seems to be a problem? What kind of types are possible to return with the "lookup"-function? Quote:
Actually I don't need to return a bool, I just want to make a new entry into the dictionary, like I did for the polynomial coefficients, and dependent on the value/word next to it (preferably a word like "yes" or "true" or "on" instead of numbers 1 or 0), the solver should decide, if variable thermal properties are used or not. If a certain keyword is found, then the bool "varProp" is switched to "true". Since you pointed out that in the controlDict the option "runTimeModifiable" is defined as Switch: Do all the possible option in the controlDict have to be defined (like in Time.H)? In the transportDict that obviously is not necessary, therefore there should be a possibilty to define keywords. I suppose. Thanks, stockzahn |
||||
February 13, 2019, 11:03 |
|
#4 |
Member
Join Date: Oct 2018
Location: France
Posts: 34
Rep Power: 7 |
I've just got a hint from Mowgli, for all of you who calculate polynoms like I do it in the first code snippet of my post #1: Have a look to the Horner scheme! Decreased the calculation time distinctly.
|
|
February 13, 2019, 23:42 |
|
#5 | |||
Senior Member
|
Hi, stockzahn
Quote:
Quote:
If there is no entry, the value is false. lookupOrDefault function can read the entry or can set the default value even if there is no entry. lookup function will cause an error with message when the entry does not exist. https://github.com/OpenFOAM/OpenFOAM...emplates.C#L32 --- The other example (not tested) Code:
Switch varProps = runTime.controlDict().lookupOrDefault<Switch>(”variable", false) if (varProps) { //do something } Quote:
https://github.com/OpenFOAM/OpenFOAM...tionary.H#L320 --- lookup-func can be used (not tested) Code:
bool varProps = false; varProps << runTime.controlDict().lookup("variable") if (varProps) { // do something } |
||||
February 15, 2019, 08:03 |
|
#6 |
Member
Join Date: Oct 2018
Location: France
Posts: 34
Rep Power: 7 |
Thanks again snak,
that one: Code:
bool varProps = runTime.controlDict().lookupOrDefault(”variable", false) worked perefctly - as I wanted it. I put it in the transportProperties though: Code:
bool varProps = transportProperties.lookupOrDefault(”variable", false); Thanks a lot again! stockzahn |
|
February 16, 2019, 07:08 |
|
#7 |
Senior Member
|
Congrats on your success, stockzahn.
bool can only handle yes or no. Using Switch class will improve the usability with on/off and yes/no. If you have any chance, could you test the Switch version? Code:
Switch varProps = transportProperties.lookupOrDefault<Switch>("variable", false); |
|
February 18, 2019, 03:41 |
|
#8 |
Member
Join Date: Oct 2018
Location: France
Posts: 34
Rep Power: 7 |
Dear snak,
Both versions (bool and Switch) work with all keyword pairs (on/off, true/false, yes/no). Switch: Code:
Switch varProps = transportProperties.lookupOrDefault<Switch>("variable", false); bool: Code:
bool varProps = transportProperties.lookupOrDefault("variable", false); Last edited by stockzahn; February 18, 2019 at 08:23. |
|
Thread Tools | Search this Thread |
Display Modes | |
|
|
Similar Threads | ||||
Thread | Thread Starter | Forum | Replies | Last Post |
New Solver for a modified pUCoupledFoam | noiseMaker | OpenFOAM | 0 | August 2, 2016 06:14 |
Error with modified interFoam solver | robertmaier9 | OpenFOAM | 3 | September 30, 2011 06:26 |
Working directory via command line | Luiz | CFX | 4 | March 6, 2011 20:02 |
why the solver reject it? Anyone with experience? | bearcat | CFX | 6 | April 28, 2008 14:08 |
Laplacian viscous stress term in compressible solver | jelmer | OpenFOAM Running, Solving & CFD | 3 | June 23, 2006 07:31 |