[Burak_1984]

Hi Foamers


I am trying to merge the TurbineSiting tutorial in Simple foam tutorials with the BosinessqSimpleFoam tutorials hot room tutorial.I got stuck somewhere eventually.The error start with reading/calculating flux phi

Basically I would like to have a turbinesiting case with the heat and rising thermal air. So I migrated most of my files from the turbinesiting.I am using the blockmesh dictionary and the snappyhexmesh dictionarias as it is.So my mesh is identical to this example

Then I also added the different files from the buoyantBoussinesqSimpleFoam example;These files include Alpat , p_rgh ,T .I have changed the names here to reflect the boundary condition names in the turbine siting tutorial.


In fvSchemes I only changed the Gauss Linear limited >Gauss linear limited corrected 0.33; as it is in the Turbine siting scheme.However I copied tis fvSchemes from the hotroom example which has only the addition of div(phi,T) bounded Gauss upwind; in divergence schemes.


For fvSolutions I took it directly from the hot room example without any changes.My control dictionary is taken from that example likewise.


In order to not to complicate thing I didn’t do any paralelization at first.So what I do while runing is pretty simple;


BlockMesh
snappyHexMesh -overwrite
buoyantBoussinesqSimpleFoam

I checked my mesh and it seems okay and understand the patches etc.I am getting this error at the end;


“Create time
Create mesh for time = 0

SIMPLE: convergence criteria
field p_rgh tolerance 0.01
field U tolerance 0.0001
field T tolerance 0.01
field "(k|epsilon|omega)" tolerance 0.001

Reading thermophysical properties

Reading field T
Reading field p_rgh
Reading field U
Reading/calculating face flux field phi


#0 Foam::error:rintStack(Foam::Ostream&) at ??:?
#1 Foam::sigFpe::sigHandler(int) at ??:?
#2 ? in "/lib/x86_64-linux-gnu/libc.so.6"
#3 Foam::surfaceInterpolation::makeWeights() const at ??:?
#4 Foam::surfaceInterpolation::weights() const at ??:?
#5 Foam::surfaceInterpolationScheme<Foam::Vector<doub le> >::dotInterpolate(Foam::GeometricField<Foam::Vecto r<double>, Foam::fvsPatchField, Foam::surfaceMesh> const&, Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> const&) const at ??:?
#6 Foam::fvc::flux(Foam::GeometricField<Foam::Vector< double>, Foam::fvPatchField, Foam::volMesh> const&) at ??:?
#7 ? in "/opt/openfoam5/platforms/linux64GccDPInt32Opt/bin/buoyantBoussinesqSimpleFoam"
#8 __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#9 ? in "/opt/openfoam5/platforms/linux64GccDPInt32Opt/bin/buoyantBoussinesqSimpleFoam"
Floating point exception (core dumped)
“

I first suspectedthere is an error in Alpat because when I open the Parafoam and click this value it crashes.Then I changed the inlet value from zero gradient to constant zero.I can look at it in the paraview but I still have the error.

Could you guys help and guide me in constructing a working example?

My second step in this effort is to divide the mesh into patches and give different patch T values for uneven heating.

I hope I was clear at explaning myself.

Thanks for the answers and effort.Greatly appriciate it.


Note:Here I also included my files in a zip file ready to run.(4.5MB).Also attached you can see my alphat and U files in txt format.


https://drive.google.com/open?id=1ZK...JSZ91o4Pah50Tp

Regard
Burak K


p { margin-bottom: 0.1in; line-height: 115%; }

[massive_turbulence]

Not sure how much this helps but using compressibleTwoPhaseEulerFoam might also work.

[Burak_1984]

Originally Posted by massive_turbulence

Not sure how much this helps but using compressibleTwoPhaseEulerFoam might also work.

Hi thanks for the rapid response.Evey bit helps I guess .

I checked the model but it seems to complicate my situation.I would like to maintain the structure of folders etc to some degree to still have a grip of it.

But thanks for the effort!

[massive_turbulence]

Originally Posted by Burak_1984

Hi thanks for the rapid response.Evey bit helps I guess .

I checked the model but it seems to complicate my situation.I would like to maintain the structure of folders etc to some degree to still have a grip of it.

But thanks for the effort!

To make changes like this you'd have to know the two solvers inside and out, meaning that you need to probably change the source and recompile it to work. Changing the object files will probably not work as something is obviously terribly wrong.

[mAlletto]

Hello. Your error comes from the file createPhi.H. This file creates the face fluxes. It interpolates the velocity U at the faces und make the do product with the surface normal vector.

See also https://www.openfoam.com/documentati...8H_source.html


So from a first guess I think something is wrong with U. (I did not have the time to look at the zip file you provided)

[mAlletto]

I just executed checkMesh on your mesh and this is the output:

  Checking basic cellZone addressing...
    No cellZones found.

Checking geometry...
    Overall domain bounding box (581321 4.78537e+06 938.71) (582290 4.78624e+06 1500)
    Mesh has 3 geometric (non-empty/wedge) directions (1 1 1)
    Mesh has 3 solution (non-empty) directions (1 1 1)
    Boundary openness (1.27842e-16 -3.45359e-17 0) OK.
 ***High aspect ratio cells found, Max aspect ratio: 1.81704e+152, number of cells 44006
  <<Writing 44006 cells with high aspect ratio to set highAspectRatioCells
 ***Zero or negative face area detected.  Minimum area: 0
  <<Writing 174456 zero area faces to set zeroAreaFaces
 ***Zero or negative cell volume detected.  Minimum negative volume: -4.83139e-08, Number of negative volume cells: 67083
  <<Writing 67083 zero volume cells to set zeroVolumeCells
    Mesh non-orthogonality Max: 180 average: 67.1048
   *Number of severely non-orthogonal (> 70 degrees) faces: 4094.
 ***Number of non-orthogonality errors: 210990.
  <<Writing 215084 non-orthogonal faces to set nonOrthoFaces
 ***Error in face pyramids: 22095 faces are incorrectly oriented.
  <<Writing 20667 faces with incorrect orientation to set wrongOrientedFaces
 ***Max skewness = 5.43226e+143, 38549 highly skew faces detected which may impair the quality of the results
  <<Writing 38549 skew faces to set skewFaces
    Coupled point location match (average 0) OK.

Failed 6 mesh checks.

I would suggest to make a better mesh first

[mAlletto]

I attached a working turbine siting tutorial. I removed the terrain.stl in order that i could upload the zipped file

[tobiasS]

Hey Burak,
The problem is in your initial value of p.. In such solvers (with rho calculation) it cannot be zero.. So it has to be something like 1e5..
Regards..


Hussam

[Burak_1984]

Originally Posted by mAlletto

I attached a working turbine siting tutorial. I removed the terrain.stl in order that i could upload the zipped file

Hi Michael thanks for the effort! that was nice of you.


I will try to do a better job at mesh building next time.I never thought it would be the problem.


Now I have to go back and see where I was short in my files.

[Burak_1984]

Originally Posted by tobiasS

Hey Burak,
The problem is in your initial value of p.. In such solvers (with rho calculation) it cannot be zero.. So it has to be something like 1e5..
Regards..


Hussam

Thank you Hussam I will check my example and try a small value of P instead of zero.Thanks!

[tobiasS]

Hey Burak,
You're welcome.. I would use just the real absolute pressure like 101300 Pa.. Because this value will be used to calculate the density..

[Burak_1984]

Originally Posted by tobiasS

Hey Burak,
You're welcome.. I would use just the real absolute pressure like 101300 Pa.. Because this value will be used to calculate the density..

Hi


I added the pressures and it actually runs;


p_rgh in OF 1.7


From here I understand the vector quantity of g is negative so any increase in height is increasing the p_rgh value above P. For my case with 935m height this value (P-rgh) corresponds to 112524 Pa.I entered these values and ran the case.


Now I ran into trouble again see the example below Any help would be appricated.

[Burak_1984]

Originally Posted by mAlletto

I attached a working turbine siting tutorial. I removed the terrain.stl in order that i could upload the zipped file

Hi Michael


Sorry I don't want to be the needy type but in your attached example when I change the terrain value from 300K to 320K it seems to run very fast and crash at the 8th time step.I also changed my pressure values P to 101300 and the P_rgh to 113000 ish.Otherwise no change to your example.



I think its a problem of P_rgh calculation since it was the next step when the iteration gets canceled and GAMG solver is for solving such.The floating point error is usually an indication of division by zero or a very big number probably.


time step continuity errors : sum local goes to very big numbers why is that? I guess at the end it's division by zero but in which boundary condition I couldn't figure out.



Can you get a quick and give some insight again? Thanks a lot in advance.


Time = 0.7

smoothSolver: Solving for Ux, Initial residual = 0.657214783401, Final residual = 0.0291041043696, No Iterations 2
smoothSolver: Solving for Uy, Initial residual = 0.905530396866, Final residual = 0.057515992305, No Iterations 2
smoothSolver: Solving for Uz, Initial residual = 0.677630489646, Final residual = 0.0540805547183, No Iterations 1
smoothSolver: Solving for T, Initial residual = 0.195037331875, Final residual = 0.00973081249398, No Iterations 2
GAMG: Solving for p_rgh, Initial residual = 0.329531552134, Final residual = 0.00375071802158, No Iterations 1
time step continuity errors : sum local = 1.21035331236e+42, global = -2.63344335341e+41, cumulative = -2.63344335341e+41
smoothSolver: Solving for epsilon, Initial residual = 1, Final residual = 0.0905756200369, No Iterations 1
smoothSolver: Solving for k, Initial residual = 1, Final residual = 5.63098879497e-06, No Iterations 1
ExecutionTime = 6.29 s ClockTime = 7 s

Time = 0.8

smoothSolver: Solving for Ux, Initial residual = 2.03765961123e-05, Final residual = 1.29769010614e-06, No Iterations 2
smoothSolver: Solving for Uy, Initial residual = 4.11436268038e-06, Final residual = 5.96627015157e-08, No Iterations 2
smoothSolver: Solving for Uz, Initial residual = 5.02213173347e-06, Final residual = 4.94774532627e-07, No Iterations 1
smoothSolver: Solving for T, Initial residual = 0.117876455954, Final residual = 0.00190403878936, No Iterations 2
#0 Foam::error:rintStack(Foam::Ostream&) at ??:?
#1 Foam::sigFpe::sigHandler(int) at ??:?
#2 ? in "/lib/x86_64-linux-gnu/libc.so.6"
#3 Foam::GAMGSolver::scale(Foam::Field<double>&, Foam::Field<double>&, Foam::lduMatrix const&, Foam::FieldField<Foam::Field, double> const&, Foam::UPtrList<Foam::lduInterfaceField const> const&, Foam::Field<double> const&, unsigned char) const at ??:?
#4 Foam::GAMGSolver::Vcycle(Foam::PtrList<Foam::lduMa trix::smoother> const&, Foam::Field<double>&, Foam::Field<double> const&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::Field<double>&, Foam::PtrList<Foam::Field<double> >&, Foam::PtrList<Foam::Field<double> >&, unsigned char) const at ??:?
#5 Foam::GAMGSolver::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at ??:?
#6 Foam::fvMatrix<double>::solveSegregated(Foam::dict ionary const&) at ??:?
#7 Foam::fvMatrix<double>::solve(Foam::dictionary const&) in "/opt/openfoam5/platforms/linux64GccDPInt32Opt/bin/buoyantBoussinesqSimpleFoam"
#8 Foam::fvMatrix<double>::solve() in "/opt/openfoam5/platforms/linux64GccDPInt32Opt/bin/buoyantBoussinesqSimpleFoam"
#9 ? in "/opt/openfoam5/platforms/linux64GccDPInt32Opt/bin/buoyantBoussinesqSimpleFoam"
#10 __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#11 ? in "/opt/openfoam5/platforms/linux64GccDPInt32Opt/bin/buoyantBoussinesqSimpleFoam"
Floating point exception (core dumped)

[Burak_1984]

After some hours spent I am convinced that my T boundary conditions don't hold up.The uneven situation causes a massive blow up because I think it's not mitigated somehow.I tried using the inletoutlet boundary condition but it didn't help.


Could somebody give some ideas;

#include "include/initialConditions"

dimensions [0 0 0 1 0 0 0];

internalField uniform 300;

boundaryField
{
inlet
{
type uniformFixedValue;
uniformValue constant 300;
}

outlet
{

type uniformFixedValue;
uniformValue constant 300;
}

terrain
{
//type fixedGradient; //fixed heat flux b.c.
//type uniformFixedValue;
//uniformValue constant 310;


type uniformFixedValue;
uniformValue constant 310;
}

ground
{
type zerogradient;
}

#include "include/sideAndTopPatches"
}



UPDATE On this one: I managed to keep the run stable by changing some relaxation factors in the fvsolutions.What I did was use the 1/3 relaxation factors for p_rgh and U and T.It didn't converge yet but it doesn'T blow up as well and the timestep continuity errors are in check as well.