Accelerating a turbulent pipe flow on nekRS.

[isurug1997]

Hi all,
I hope you are all doing well. I have recently managed to run Nek5000 cases where I am accelerating and decelerating a turbulent pipe flow (Retau = 170 - 350 and vise-versa) and now I would like to move on to nekRS and simulate the same thing. I would really appreciate it if I could get a bit of assistance with a starting point so that I could improve it to suit my scenario.

The way I managed to accelerate the flow on Nek5000 was through implementing a bulk velocity change over time (Ub = param(55)) and if I could get some insight on how to do the same thing on nekRS (implementing a bulk velocity change over time), I could start from there.

As always, thank you for your time.

Yours sincerely,
Isuru.

[yslan]

Hi,

We do have the constant flow rate control for periodic flow. See the example

nekRS/examples/turbPipePeriodic/turbPipe.par

Line 16 in 4f87e0e

constFlowRate = meanVelocity=1.0 + direction=Z

where you can choose either meanVelocity=1.0 or meanvolumetricflow=1.0.

You should be able to modify the value at run time in the udf file

void UDF_ExecuteStep(nrs_t *nrs, dfloat time, int tstep) {
  const dfloat flowRate = sin(M_PI * time);
  platform->options.setArgs("FLOW RATE", to_string_f(flowRate));
}

You might also want to use the future time as well

timeNew = time + nrs->dt[0]

Search flowRate in logfile

[isurug1997]

Hi,
Thank you so much for your quick response. As you mentioned I have updated the .udf file as below. Unfortunately, once I inspected the logfile for "flowRate", it appears not to change over time. Any possible suggestions I could try to fix this issue? I really appreciate your kind assistance so far with this case.

.................................................................................. .udf file .............................................................................................................................................................
//

// nekRS User Defined File

//

#include "udf.hpp"
#include
static dfloat Re_tau;
static dfloat zLength;
static dfloat dt0;

#include "ci.inc"
#ifdef okl
#include "turbPipe.oudf"
#endif

void userf(nrs_t *nrs, dfloat time, occa::memory o_U, occa::memory o_FU)
{
mesh_t *mesh = nrs->meshV;
dfloat mue, rho;
platform->options.getArgs("VISCOSITY", mue);
platform->options.getArgs("DENSITY", rho);
const dfloat RE_B = rho / mue;
const dfloat DPDZ = 4 * (2 * Re_tau / RE_B) * (2 * Re_tau / RE_B);

occa::memory o_FUz = o_FU + 2 * nrs->fieldOffset * sizeof(dfloat);
platform->linAlg->fill(mesh->Nlocal, DPDZ, o_FUz);
}

void UDF_Setup0(MPI_Comm comm, setupAide &options)
{
platform->par->extract("casedata", "re_tau", Re_tau);
platform->par->extract("casedata", "zlength", zLength);
options.getArgs("CI-MODE", ciMode);
if (ciMode) {
ciSetup(comm, options);
}

if (platform->options.compareArgs("BUILD ONLY", "FALSE")) {
double *const nek_cb_scnrs = (double *)nek::ptr("cb_scnrs");
nek_cb_scnrs[0] = zLength;
}
platform->options.getArgs("DT", dt0);
}

void UDF_Setup(nrs_t *nrs)
{
//if (platform->options.compareArgs("CONSTANT FLOW RATE", "FALSE")) {
// assign function pointer to drive flow by constant mean pressure gradient
//udf.uEqnSource = &userf;
//}
}

**void UDF_ExecuteStep(nrs_t *nrs, dfloat time, int tstep)
{
dfloat t1 = 0.08;
dfloat t2 = 1.08;
dfloat t3 = 2.08;
dfloat delta_t = (t2-t1);
dfloat delta_v = 1;
dfloat delta_a = delta_v/delta_t;
dfloat Ub0 = 1;
dfloat Ub1 = 2;

const dfloat timeNew = time + nrs->dt[0];

if (timeNew >= t1) {
const dfloat flowRate = Ub0;
platform->options.setArgs("FLOW RATE", to_string_f(flowRate));
std::cout << "Flowrate = " << flowRate << std::endl;
}
else if (t1 < timeNew && timeNew <= t2) {
const dfloat flowRate = Ub0 + (timeNew - t1) * delta_a;
platform->options.setArgs("FLOW RATE", to_string_f(flowRate));
std::cout << "Flowrate = " << flowRate << std::endl;
}
else if (t2 < timeNew && timeNew <= t3) {
const dfloat flowRate = Ub1 - (time - t2) * delta_a;
platform->options.setArgs("FLOW RATE", to_string_f(flowRate));
std::cout << "Flowrate = " << flowRate << std::endl;
}
else if (timeNew < t3) {
const dfloat flowRate = Ub0;
platform->options.setArgs("FLOW RATE", to_string_f(flowRate));
std::cout << "Flowrate = " << flowRate << std::endl;
}**

if (tstep < 5) {
nrs->dt[0] = 0.25 * dt0;
} else {
nrs->dt[0] = dt0;
}
if (nrs->lastStep) {
if (ciMode) {
ciTestErrors(nrs, time, tstep);
}
}
}

Yours sincerely,
Isuru.

[yslan]

Did you turn the flowrate control on in par file?
Something like this?

 constFlowRate = meanVelocity=1.0 + direction=Z 

[isurug1997]

Hi,
Thank you for your quick response.

Yes, It's set as follows. Now I can see the flowRate and uBulk in the logfile as well. Initially what I thought was commenting it out in the .par file would allow me to apply a variable meanVelocity as intended rather than meanVelocity set to a constant. As I am not still observing uBulk/flowRate change over the designated for-loop time intervals, Instead of setting "meanVelocity=1.0" is it required to define "meanVelocity" so that it changes over time?

.................................................................................. .par file .............................................................................................................................................................
[GENERAL]
#verbose = true
polynomialOrder = 7
#startFrom = "restart.fld"
stopAt = endTime
endTime = 200
dt = 5.0e-3
timeStepper = tombo2
subCyclingSteps = 2

writeControl = simulationTime
writeInterval = 20

regularization = hpfrt + nModes=1 + scalingCoeff=10
constFlowRate = meanVelocity=1.0 + direction=Z

[PRESSURE]
residualTol = 1e-04

[VELOCITY]
boundaryTypeMap = wall
viscosity = -19000
density = 1.0
residualTol = 1e-06

[CASEDATA]
Re_tau = 550
zLength = 6

Again, I really appreciate your support on this.
Yours sincerely,
Isuru.

[yslan]

is it required to define "meanVelocity" so that it changes over time?

You choose either meanVelocity or meanvolumetricflow with initial value.
Then, overwrite the value in UDF_ExecuteStep

Even though we print the corrected flow rate, I still recommend to monitor the mean velocity to make sure it's computed the right value.

[isurug1997]

Hi,
Thank you for your valuble response. It is clear to me now.
The only issue that resides now is that the flowRate/uBulk is not appearing to change over time. I have been monitoring both the "flowRate" and "uBulk" and if possible, I would really appreciate your comments on the latest .par and the .udf settings of the current nekRS code.

.................................................................................. .par file .............................................................................................................................................................
[GENERAL]
#verbose = true
polynomialOrder = 7
#startFrom = "restart.fld"
stopAt = endTime
endTime = 200

dt = 5.0e-3
timeStepper = tombo2
subCyclingSteps = 2

writeControl = simulationTime
writeInterval = 20

regularization = hpfrt + nModes=1 + scalingCoeff=10
constFlowRate = meanVelocity=1.0 + direction=Z

[PRESSURE]
residualTol = 1e-04

[VELOCITY]
boundaryTypeMap = wall
viscosity = -19000
density = 1.0
residualTol = 1e-06

[CASEDATA]
Re_tau = 550
zLength = 6

.................................................................................. .udf file .............................................................................................................................................................
//
// nekRS User Defined File
//
#include "udf.hpp"
#include

static dfloat Re_tau;
static dfloat zLength;
static dfloat dt0;

#include "ci.inc"

#ifdef okl

#include "turbPipe.oudf"

#endif

void userf(nrs_t *nrs, dfloat time, occa::memory o_U, occa::memory o_FU)
{
mesh_t *mesh = nrs->meshV;
dfloat mue, rho;
platform->options.getArgs("VISCOSITY", mue);
platform->options.getArgs("DENSITY", rho);
const dfloat RE_B = rho / mue;
const dfloat DPDZ = 4 * (2 * Re_tau / RE_B) * (2 * Re_tau / RE_B);

occa::memory o_FUz = o_FU + 2 * nrs->fieldOffset * sizeof(dfloat);
platform->linAlg->fill(mesh->Nlocal, DPDZ, o_FUz);
}

void UDF_Setup0(MPI_Comm comm, setupAide &options)
{
platform->par->extract("casedata", "re_tau", Re_tau);
platform->par->extract("casedata", "zlength", zLength);

options.getArgs("CI-MODE", ciMode);
if (ciMode) {
ciSetup(comm, options);
}

if (platform->options.compareArgs("BUILD ONLY", "FALSE")) {
double *const nek_cb_scnrs = (double *)nek::ptr("cb_scnrs");
nek_cb_scnrs[0] = zLength;
}

platform->options.getArgs("DT", dt0);
}

void UDF_Setup(nrs_t *nrs)
{
//if (platform->options.compareArgs("CONSTANT FLOW RATE", "FALSE")) {
// assign function pointer to drive flow by constant mean pressure gradient
//udf.uEqnSource = &userf;
//}
}

void UDF_ExecuteStep(nrs_t *nrs, dfloat time, int tstep)
{
const dfloat t1 = 0.1;
const dfloat t2 = 1.1;
const dfloat t3 = 2.1;
dfloat delta_t = (t2-t1);
const dfloat delta_v = 1;
dfloat delta_a = delta_v/delta_t;
const dfloat Ub0 = 1;
const dfloat Ub1 = 2;

dfloat timeNew = time + nrs->dt[0];

if (timeNew >= t1) {
dfloat flowRate = Ub0;
platform->options.setArgs("FLOW RATE", to_string_f(flowRate));
//std::cout << "Flowrate = " << flowRate << std::endl;
}
else if (t1 < timeNew && timeNew <= t2) {
dfloat flowRate = Ub0 + (timeNew - t1) * delta_a;
platform->options.setArgs("FLOW RATE", to_string_f(flowRate));
//std::cout << "Flowrate = " << flowRate << std::endl;
}
else if (t2 < timeNew && timeNew <= t3) {
dfloat flowRate = Ub1 - (timeNew - t2) * delta_a;
platform->options.setArgs("FLOW RATE", to_string_f(flowRate));
//std::cout << "Flowrate = " << flowRate << std::endl;
}
else if (timeNew < t3) {
dfloat flowRate= Ub0;
platform->options.setArgs("FLOW RATE", to_string_f(flowRate));
//std::cout << "Flowrate = " << flowRate << std::endl;
}

if (tstep < 5) {
nrs->dt[0] = 0.25 * dt0;
} else {
nrs->dt[0] = dt0;
}

if (nrs->lastStep) {
if (ciMode) {
ciTestErrors(nrs, time, tstep);
}
}
}

Thank you.
Yours sincerely,
Isuru.