TurbulentFlowSimulation-HPC/MixingLengthModel.cpp at master · ludcila/TurbulentFlowSimulation-HPC · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
#include "MixingLengthModel.h"
#include <iostream>


MixingLengthModel::MixingLengthModel(const Parameters & parameters) :
	_parameters(parameters)
{
}
MixingLengthModel::~MixingLengthModel() {}

/**************************************************************************************
	Lm = Kh
**************************************************************************************/

LmKh::LmKh(const Parameters & parameters) :
	MixingLengthModel(parameters)
{
}

FLOAT LmKh::at(TurbulentFlowField & flowField, int i, int j) {
	return 0.41 * flowField.getWallDistance(i, j);
}

FLOAT LmKh::at(TurbulentFlowField & flowField, int i, int j, int k) {
	return 0.41 * flowField.getWallDistance(i, j, k);
}


/**************************************************************************************
	Lm = min(Kh, 0.09 delta)
	Assuming laminar flat plate Blasius boundary layer
**************************************************************************************/

LmLaminarFlatPlate::LmLaminarFlatPlate(const Parameters & parameters) :
	MixingLengthModel(parameters)
{
}

FLOAT LmLaminarFlatPlate::at(TurbulentFlowField & flowField, int i, int j) {
	FLOAT kh = 0.41 * flowField.getWallDistance(i, j);
	FLOAT x = _parameters.meshsize->getPosX(i, j) + 0.5 * _parameters.meshsize->getDx(i, j);
	FLOAT centerLineVelocity = flowField.getCenterLineVelocity()[i];
	// FLOAT centerLineVelocity = flowField.getVelocity().getVector(i, _parameters.geometry.sizeY / 2)[0];
	FLOAT downstreamRe = centerLineVelocity * x * _parameters.flow.Re;
	FLOAT delta = 4.91 * x / pow(downstreamRe, 0.5);
	return std::min(kh, 0.09 * delta);
}

FLOAT LmLaminarFlatPlate::at(TurbulentFlowField & flowField, int i, int j, int k) {
	FLOAT kh = 0.41 * flowField.getWallDistance(i, j, k);
	FLOAT x = _parameters.meshsize->getPosX(i, j, k) + 0.5 * _parameters.meshsize->getDx(i, j, k);
	FLOAT centerLineVelocity = flowField.getCenterLineVelocity()[i];
	// FLOAT centerLineVelocity = flowField.getVelocity().getVector(i, _parameters.geometry.sizeY / 2, _parameters.geometry.sizeZ / 2)[0];
	FLOAT downstreamRe = centerLineVelocity * x * _parameters.flow.Re;
	FLOAT delta = 4.91 * x / pow(downstreamRe, 0.5);
	return std::min(kh, 0.09 * delta);
}


/**************************************************************************************
	Lm = min(Kh, 0.09 delta)
	Assuming turbulent flat plate Blasius boundary layer
**************************************************************************************/

LmTurbulentFlatPlate::LmTurbulentFlatPlate(const Parameters & parameters) :
	MixingLengthModel(parameters)
{
}

FLOAT LmTurbulentFlatPlate::at(TurbulentFlowField & flowField, int i, int j) {
	FLOAT kh = 0.41 * flowField.getWallDistance(i, j);
	FLOAT x = _parameters.meshsize->getPosX(i, j) + 0.5 * _parameters.meshsize->getDx(i, j);
	FLOAT centerLineVelocity = flowField.getCenterLineVelocity()[i];
	// FLOAT centerLineVelocity = flowField.getVelocity().getVector(i, _parameters.geometry.sizeY / 2)[0];
	FLOAT downstreamRe = centerLineVelocity * x * _parameters.flow.Re;
	FLOAT delta = 0.382 * x / pow(downstreamRe, 0.2);
	return std::min(kh, 0.09 * delta);
}

FLOAT LmTurbulentFlatPlate::at(TurbulentFlowField & flowField, int i, int j, int k) {
	FLOAT kh = 0.41 * flowField.getWallDistance(i, j, k);
	FLOAT x = _parameters.meshsize->getPosX(i, j, k) + 0.5 * _parameters.meshsize->getDx(i, j, k);
	FLOAT centerLineVelocity = flowField.getCenterLineVelocity()[i];
	// FLOAT centerLineVelocity = flowField.getVelocity().getVector(i, _parameters.geometry.sizeY / 2, _parameters.geometry.sizeZ / 2)[0];
	FLOAT downstreamRe = centerLineVelocity * x * _parameters.flow.Re;
	FLOAT delta = 0.382 * x / pow(downstreamRe, 0.2);
	return std::min(kh, 0.09 * delta);
}