Implementation of a turbulent inflow boundary condition for LES - - PowerPoint PPT Presentation

▶

May 31, 2023 272 likes •503 views

Background Test case Mean velocity Fluctuations Conlusion Implementation of a turbulent inflow boundary condition for LES based on a vortex method Nina Gall Jrgensen Technical University of Denmark Department of Civil Engineering

SLIDE 1

Background Test case Mean velocity Fluctuations Conlusion

Implementation of a turbulent inflow boundary condition for LES based on a vortex method

Nina Gall Jørgensen

Technical University of Denmark Department of Civil Engineering

22-10-2012

Nina Gall Jørgensen Implementation of a boundary condition in OpenFOAM 09-09-2012 1 / 21

SLIDE 2

Background Test case Mean velocity Fluctuations Conlusion

Motivation

Atmospheric boundary layer flows

Nina Gall Jørgensen Implementation of a boundary condition in OpenFOAM 09-09-2012 2 / 21

SLIDE 3

Background Test case Mean velocity Fluctuations Conlusion

Mean flow profile estimations

Obtain mean flow data from experiments or RANS simulations

Logarithmic law U(z) = u∗ κ ln(z − d0 z0 ) Power law U(z) = U10( z 10)α Relation between z0 and α α = 1 ln(zref/z0)) Turbulent kinetic energy k = 3 2(UIu)2

Mean velocity

Nina Gall Jørgensen Implementation of a boundary condition in OpenFOAM 09-09-2012 3 / 21

SLIDE 4

Background Test case Mean velocity Fluctuations Conlusion

Turbulent velocity inlet boundary conditions in OpenFoam

(University of Rostock, Germany - inflowGenerator)

tubulentInlet inflowGenerator mapped

Nina Gall Jørgensen Implementation of a boundary condition in OpenFOAM 09-09-2012 4 / 21

SLIDE 5

Background Test case Mean velocity Fluctuations Conlusion

Test Case - boxPowerLawTurbulentLES

(See Report for detailed information on the test case)

0 directory: U, p, k, nuSgs, nuTilda constant: LESProperties - LESModel oneEqEddy system: controlDict - application pisoFoam

Nina Gall Jørgensen Implementation of a boundary condition in OpenFOAM 09-09-2012 5 / 21

SLIDE 6

Background Test case Mean velocity Fluctuations Conlusion

Copy and Compile Parabolic Velocity BC

Compile guidelines given in course lecture notes: Copy and compile a boundary condition, and a deeper look at parabolicInlet

Copy the BC from the OpenFOAM-extend project at SourceForge:

OF21x mkdir -p $WM_PROJECT_USER_DIR/src/finiteVolume/fields/fvPatchFields/derived cd $WM_PROJECT_USER_DIR/src/finiteVolume/fields/fvPatchFields/derived svn checkout http://openfoam-extend.svn.sourceforge.net/svnroot/\

penfoam-extend/trunk/Core/OpenFOAM-1.5-dev/src/finiteVolume/\

fields/fvPatchFields/derived/parabolicVelocity/ cd parabolicVelocity

Create the Make/files file and include following:

parabolicVelocityFvPatchVectorField.C LIB = $(FOAM_USER_LIBBIN)/libmyBCs

Create the Make/options file and include following:

EXE_INC = \

I$(LIB_SRC)/finiteVolume/lnInclude

EXE_LIBS =

Compile the dynamic library:

wmake libso

Nina Gall Jørgensen Implementation of a boundary condition in OpenFOAM 09-09-2012 6 / 21

SLIDE 7

Background Test case Mean velocity Fluctuations Conlusion

Modifications to Parabolic Velocity BC

Copy parabolicVelocity to powerLawTurbulentVelocity:

cd $WM_PROJECT_USER_DIR/src/finiteVolume/fields/ fvPatchFields/derived cp -r parabolicVelocity powerLawTurbulentVelocity cd powerLawTurbulentVelocity rename .H and .C files: rename parabolic powerLawTurbulent * Modify the Make/files file: powerLawTurbulentVelocityFvPatchVectorField.C LIB = $(FOAM_USER_LIBBIN)/libmypowerLawTurbulent Change class name everywhere in .H and .C files:

sed -i s/parabolic/powerLawTurbulent/g powerLawTurbulentVelocityFvPatchVectorField.C sed -i s/parabolic/powerLawTurbulent/g powerLawTurbulentVelocityFvPatchVectorField.H

Compile the dynamic library with wmake libso

Nina Gall Jørgensen Implementation of a boundary condition in OpenFOAM 09-09-2012 7 / 21

SLIDE 8

Background Test case Mean velocity Fluctuations Conlusion

Modifications to Parabolic Velocity BC

Add private data for a and alpha in the .H and .C files in same way as scalar

maxValue. Erase everywhere maxValue

U(z) = azα

Change the operator to power law profile

scalarField coord = (c & y_); vectorField Uinf = (n_a_pow(coord,alpha_));

Compile the dynamic library wclean libso, wmake libso Run boxPowerLawTurbulentLES test case

(for this write ::operator instead of Uinf)

Nina Gall Jørgensen Implementation of a boundary condition in OpenFOAM 09-09-2012 8 / 21

SLIDE 9

Background Test case Mean velocity Fluctuations Conlusion

Vortex Method

Generation of 2D small scale flow structures

(Mathey et al., 2006 - Assessment of vortex method for LES inlet conditions)

Velocity at inlet of channel flow d¯ ω dt = ¯ ω · ∇¯ u + ν∇2¯ ω Vorticity equation solved with particle discretisation ¯ u(¯ x) = 1 2π

Γi (¯ xi − ¯ x) × ¯ z |¯ x − ¯ xi|2 (1 − e

|¯ x−¯ xi|2 2σ2

)e

|¯ x−¯ xi|2 2σ2

Nina Gall Jørgensen Implementation of a boundary condition in OpenFOAM 09-09-2012 9 / 21

SLIDE 10

Background Test case Mean velocity Fluctuations Conlusion

Vortex Method - Application

Inflow generation for LES, Embedded LES and Hybrid RANS-LES

ANSYS Fluent workshop: Embedded LES of flow over backward facing step

Nina Gall Jørgensen Implementation of a boundary condition in OpenFOAM 09-09-2012 10 / 21

SLIDE 11

Background Test case Mean velocity Fluctuations Conlusion

Calculation of tangential fluctuations in OpenFOAM

Input: turbulent kinetic energy, k and number of vortex particles, N

Calculate sigma and Ap Init number of random vortex particles xi at inlet patch Calculate circulation ghe in all positions xi from k Calculate tangential fluctuations ux at all inlet faces Add ux to mean profile Repeat calculation for each time step

Nina Gall Jørgensen Implementation of a boundary condition in OpenFOAM 09-09-2012 11 / 21

SLIDE 12

Background Test case Mean velocity Fluctuations Conlusion

Mean turbulent kinetic energy

Add private data zref and N to the .H and .C files as for a and alpha Calculate k

scalar z0 = zref_/exp(1/alpha_); scalarField Iu = 1/log(coord/z0); scalarField Uinfx = Uinf.component(0); scalarField k = 1.5(pow(Uinfx,2)pow(Iu,2));

Nina Gall Jørgensen Implementation of a boundary condition in OpenFOAM 09-09-2012 12 / 21

SLIDE 13

Background Test case Mean velocity Fluctuations Conlusion

Control of vortex particles with σ

Check size of N and calculate sigma and inlet patch area Ap

if (N_ >= c.size()) { Info << "Error: N should be less than or equal to " << c.size() << endl; } scalarField Af = patch().magSf(); scalar Ap = sum(Af); vector bbMax = bb.max(); scalar lx = bbMax[1]; scalar ly = bbMax[2]; scalar L = (2lxly)/(lx+ly); scalar sigma = (0.1L)/2; scalar sigma2 = 2sqr(sigma);

Nina Gall Jørgensen Implementation of a boundary condition in OpenFOAM 09-09-2012 13 / 21

SLIDE 14

Background Test case Mean velocity Fluctuations Conlusion

Store tangential fluctuations for all inlet patch faces

Start if statement to execute calculation for each time step

const scalar t = this->db().time().timeOutputValue(); if (t >= t) {

The tangential fluctuations ux_sum are calculated for all inlet patch faces and stored in a vectorField ux

vectorField ux = c*0; forAll(patch(),facei) { vector x1 = c[facei]; scalar x1x = x1[1]; scalar x1y = (x1 & y_);

Nina Gall Jørgensen Implementation of a boundary condition in OpenFOAM 09-09-2012 14 / 21

SLIDE 15

Background Test case Mean velocity Fluctuations Conlusion

Circulation Γi of all vortex particles

Generate random positions xi for all vortex positions, i=1..N

srand (time(NULL) ); vector ux_sum = vector(0, 0, 0); for (int iter = 1; iter <= N_; iter ++) { int RandIndex = rand() % c.size(); vector xi = c[RandIndex]; if (xi == x1) { scalar Lxi = sqrt(patch().magSf()[iter]); scalar addxi = 0.1*Lxi; xi[1] = xi[1]+addxi; xi[2] = xi[2]+addxi; } scalar xix = xi[1]; scalar xiy = (xi & y_);

Γt(¯ x) = 4

πSk(¯

xi) 3N(2ln(3) − 3ln(2))

Nina Gall Jørgensen Implementation of a boundary condition in OpenFOAM 09-09-2012 15 / 21

SLIDE 16

Background Test case Mean velocity Fluctuations Conlusion

Circulation Γi of all vortex particles

Calculate Γi in the random vortex positions from local k

scalar dist = sqr(x1x-xix)+sqr(x1y-xiy); vector xx = xi-x1; scalar ks = k[RandIndex]; scalar gamma = 4sqrt((Foam::constant::mathematical::piApks)/ (3N_*0.117783035656383));

Nina Gall Jørgensen Implementation of a boundary condition in OpenFOAM 09-09-2012 16 / 21

SLIDE 17

Background Test case Mean velocity Fluctuations Conlusion

Calculate and sum up the tangential fluctuations

For each inlet face calculate the contribution to the tangential fluctuatins from all vortex

particles. The constant sigma is bounded to the local face edge length.

¯ u(¯ x) = 1 2π

Γi (¯ xi − ¯ x) × ¯ z |¯ x − ¯ xi|2 (1 − e

|¯ x−¯ xi|2 2σ2

)e

|¯ x−¯ xi|2 2σ2

scalar delta = sqrt(patch().magSf()[RandIndex]); if (sigma <= delta) { sigma = delta; sigma2 = 2sqr(sigma); } vector z1 = vector(1,0,0); vector ux1 = (gamma((xx ^ z1)/dist)(1-exp(-dist/sigma2))exp(-dist/sigma2)); ux_sum += ux1; }

Nina Gall Jørgensen Implementation of a boundary condition in OpenFOAM 09-09-2012 17 / 21

SLIDE 18

Background Test case Mean velocity Fluctuations Conlusion

Store tangential fluctuations for all inlet patch faces

The tangential fluctuations ux_sum are calculated for all inlet patch faces and stored in a vectorField ux. Finally, the mean and fluctuating velocities are summed up in the operator.

ux[facei] = (1/(2Foam::constant::mathematical::pi))ux_sum; } vectorField::operator=(Uinf+ux); }

Compile with wmake libso

Nina Gall Jørgensen Implementation of a boundary condition in OpenFOAM 09-09-2012 18 / 21

SLIDE 19

Background Test case Mean velocity Fluctuations Conlusion

Velocity contours at inlet patch

boxPowerLawTurbulentVelocity Streamwise Spanwise Vertical

Nina Gall Jørgensen Implementation of a boundary condition in OpenFOAM 09-09-2012 19 / 21

SLIDE 20

Background Test case Mean velocity Fluctuations Conlusion

Velocity vectors at inlet patch

Nina Gall Jørgensen Implementation of a boundary condition in OpenFOAM 09-09-2012 20 / 21

SLIDE 21

Background Test case Mean velocity Fluctuations Conlusion

Conclusion

Further modifications Modify sigma to be calculated from k and ǫ Add fluctuations to streamwise velocity Vary randomly sign of the circulation, Γi Randomize generation of vortex particles in time Add fluctuations to any chosen mean velocity profile Move definiton of k to 0/k (Use also for RASmodel) Verification Space correlations Frequency spectra

Thank you for your attention

Nina Gall Jørgensen Implementation of a boundary condition in OpenFOAM 09-09-2012 21 / 21

Implementation of a turbulent inflow boundary condition for LES based on a vortex method

Nina Gall Jørgensen

Technical University of Denmark Department of Civil Engineering

22-10-2012

Motivation

Atmospheric boundary layer flows

Mean flow profile estimations

Obtain mean flow data from experiments or RANS simulations

Logarithmic law U(z) = u∗ κ ln(z − d0 z0 ) Power law U(z) = U10( z 10)α Relation between z0 and α α = 1 ln(zref/z0)) Turbulent kinetic energy k = 3 2(UIu)2

Mean velocity

Turbulent velocity inlet boundary conditions in OpenFoam

(University of Rostock, Germany - inflowGenerator)

tubulentInlet inflowGenerator mapped

Test Case - boxPowerLawTurbulentLES

(See Report for detailed information on the test case)

0 directory: U, p, k, nuSgs, nuTilda constant: LESProperties - LESModel oneEqEddy system: controlDict - application pisoFoam

Copy and Compile Parabolic Velocity BC

Compile guidelines given in course lecture notes: Copy and compile a boundary condition, and a deeper look at parabolicInlet

Copy the BC from the OpenFOAM-extend project at SourceForge:

OF21x mkdir -p $WM_PROJECT_USER_DIR/src/finiteVolume/fields/fvPatchFields/derived cd $WM_PROJECT_USER_DIR/src/finiteVolume/fields/fvPatchFields/derived svn checkout http://openfoam-extend.svn.sourceforge.net/svnroot/\

fields/fvPatchFields/derived/parabolicVelocity/ cd parabolicVelocity

Create the Make/files file and include following:

parabolicVelocityFvPatchVectorField.C LIB = $(FOAM_USER_LIBBIN)/libmyBCs

Create the Make/options file and include following:

EXE_INC = \

EXE_LIBS =

Compile the dynamic library:

wmake libso

Modifications to Parabolic Velocity BC

Copy parabolicVelocity to powerLawTurbulentVelocity:

Compile the dynamic library with wmake libso

Modifications to Parabolic Velocity BC

Add private data for a and alpha in the .H and .C files in same way as scalar

U(z) = azα

Change the operator to power law profile

scalarField coord = (c & y_); vectorField Uinf = (n_*a_*pow(coord,alpha_));

Compile the dynamic library wclean libso, wmake libso Run boxPowerLawTurbulentLES test case

(for this write ::operator instead of Uinf)

Vortex Method

Generation of 2D small scale flow structures

(Mathey et al., 2006 - Assessment of vortex method for LES inlet conditions)

Velocity at inlet of channel flow d¯ ω dt = ¯ ω · ∇¯ u + ν∇2¯ ω Vorticity equation solved with particle discretisation ¯ u(¯ x) = 1 2π

Γi (¯ xi − ¯ x) × ¯ z |¯ x − ¯ xi|2 (1 − e

)e

Vortex Method - Application

Inflow generation for LES, Embedded LES and Hybrid RANS-LES

ANSYS Fluent workshop: Embedded LES of flow over backward facing step

Calculation of tangential fluctuations in OpenFOAM

Input: turbulent kinetic energy, k and number of vortex particles, N

Calculate sigma and Ap Init number of random vortex particles xi at inlet patch Calculate circulation ghe in all positions xi from k Calculate tangential fluctuations ux at all inlet faces Add ux to mean profile Repeat calculation for each time step

Mean turbulent kinetic energy

Add private data zref and N to the .H and .C files as for a and alpha Calculate k

scalar z0 = zref_/exp(1/alpha_); scalarField Iu = 1/log(coord/z0); scalarField Uinfx = Uinf.component(0); scalarField k = 1.5*(pow(Uinfx,2)*pow(Iu,2));

Control of vortex particles with σ

Check size of N and calculate sigma and inlet patch area Ap

Store tangential fluctuations for all inlet patch faces

Start if statement to execute calculation for each time step

const scalar t = this->db().time().timeOutputValue(); if (t >= t) {

The tangential fluctuations ux_sum are calculated for all inlet patch faces and stored in a vectorField ux

vectorField ux = c*0; forAll(patch(),facei) { vector x1 = c[facei]; scalar x1x = x1[1]; scalar x1y = (x1 & y_);

Circulation Γi of all vortex particles

Generate random positions xi for all vortex positions, i=1..N

Γt(¯ x) = 4

xi) 3N(2ln(3) − 3ln(2))

Circulation Γi of all vortex particles

Calculate Γi in the random vortex positions from local k

scalar dist = sqr(x1x-xix)+sqr(x1y-xiy); vector xx = xi-x1; scalar ks = k[RandIndex]; scalar gamma = 4*sqrt((Foam::constant::mathematical::pi*Ap*ks)/ (3*N_*0.117783035656383));

Calculate and sum up the tangential fluctuations

For each inlet face calculate the contribution to the tangential fluctuatins from all vortex

¯ u(¯ x) = 1 2π

Γi (¯ xi − ¯ x) × ¯ z |¯ x − ¯ xi|2 (1 − e

)e

scalar delta = sqrt(patch().magSf()[RandIndex]); if (sigma <= delta) { sigma = delta; sigma2 = 2*sqr(sigma); } vector z1 = vector(1,0,0); vector ux1 = (gamma*((xx ^ z1)/dist)*(1-exp(-dist/sigma2))*exp(-dist/sigma2)); ux_sum += ux1; }

Store tangential fluctuations for all inlet patch faces

The tangential fluctuations ux_sum are calculated for all inlet patch faces and stored in a vectorField ux. Finally, the mean and fluctuating velocities are summed up in the operator.

ux[facei] = (1/(2*Foam::constant::mathematical::pi))*ux_sum; } vectorField::operator=(Uinf+ux); }

Compile with wmake libso

Velocity contours at inlet patch

boxPowerLawTurbulentVelocity Streamwise Spanwise Vertical

Velocity vectors at inlet patch

scalarField coord = (c & y_); vectorField Uinf = (n_a_pow(coord,alpha_));

scalar z0 = zref_/exp(1/alpha_); scalarField Iu = 1/log(coord/z0); scalarField Uinfx = Uinf.component(0); scalarField k = 1.5(pow(Uinfx,2)pow(Iu,2));

scalar dist = sqr(x1x-xix)+sqr(x1y-xiy); vector xx = xi-x1; scalar ks = k[RandIndex]; scalar gamma = 4sqrt((Foam::constant::mathematical::piApks)/ (3N_*0.117783035656383));

scalar delta = sqrt(patch().magSf()[RandIndex]); if (sigma <= delta) { sigma = delta; sigma2 = 2sqr(sigma); } vector z1 = vector(1,0,0); vector ux1 = (gamma((xx ^ z1)/dist)(1-exp(-dist/sigma2))exp(-dist/sigma2)); ux_sum += ux1; }

ux[facei] = (1/(2Foam::constant::mathematical::pi))ux_sum; } vectorField::operator=(Uinf+ux); }