Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F84145712
NozzleBC.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sat, Sep 21, 00:04
Size
11 KB
Mime Type
text/x-c
Expires
Mon, Sep 23, 00:04 (2 d)
Engine
blob
Format
Raw Data
Handle
20946671
Attached To
R3721 ConSol
NozzleBC.cpp
View Options
//
// Created by Fabian Moenkeberg on 2019-12-19.
//
#include "NozzleBC.hpp"
void NozzleBC::initialize(MeshBase *mesh, Configuration config) {
gamma = 1.4;
U_in.resize(4);
U_out.resize(4);
double r_in = 5.8;
double u1_in = 77*r_in;
double u2_in =0;
p_in = 101325*13.5;
double E_in = p_in/(gamma - 1) + 0.5*(u1_in*u1_in+u2_in*u2_in)*r_in;
double r_out = 1.3;
double u1_out = 197;
double u2_out =0;
p_out = 101325;
double E_out = p_out/(gamma - 1) + 0.5*(u1_out*u1_out+u2_out*u2_out)*r_out;
U_in = {r_in, r_in*u1_in, r_in*u2_in, E_in};
U_out = {r_out, r_out*u1_out, r_out*u2_out, E_out};
}
void NozzleBC::updateBoundary(std::vector<std::vector<double> > &U, MeshBase *mesh){
std::vector<double> n;
int nBC = mesh->getNbc();
for (int i = 0; i < nBC; i++){
U[mesh->getBCMap(3*i+0)] = U[mesh->getBCMap(3*i+1)];
int iEdge = mesh->getBCMap(3*i+2);
if (iEdge == -1){ // internal hybrid update
U[mesh->getBCMap(3 * i + 0)] = U[mesh->getBCMap(3 * i + 1)];
}else {
std::vector<int> edgesAtBC = mesh->getEdgesAtBC(iEdge);
switch (edgesAtBC[1]) {
case 2: { // slipwall BC
// std::vector<int> edgesAtBC = mesh->getEdgesAtBC(mesh->getBCMap(3*i+2));
n = mesh->getNormal(edgesAtBC[0]);
if (edgesAtBC[2] == 1) { //means left side of edge is out & normal is INWARDS pointing
n[0] = -n[0];
n[1] = -n[1];
}
flux_slipWall(U[mesh->getBCMap(3 * i + 1)], n, U[mesh->getBCMap(3 * i + 0)]);
break;
}
case 0: { // means inflow BC
// std::vector<int> edgesAtBC = mesh->getEdgesAtBC(mesh->getBCMap(3*i+2));
n = mesh->getNormal(edgesAtBC[0]);
if (edgesAtBC[2] == 1) { //means left side of edge is out & normal is INWARDS pointing
n[0] = -n[0];
n[1] = -n[1];
}
U[mesh->getBCMap(3 * i + 0)] = flux_farfield(U[mesh->getBCMap(3 * i + 1)], n, U_in);
break;
}
case 1: { // means outflow BC
// std::vector<int> edgesAtBC = mesh->getEdgesAtBC(mesh->getBCMap(3*i+2));
n = mesh->getNormal(edgesAtBC[0]);
if (edgesAtBC[2] == 1) { //means left side of edge is out & normal is INWARDS pointing
n[0] = -n[0];
n[1] = -n[1];
}
// std::vector<double> test = flux_farfield(std::vector<double>{1,.59160797*5,0,2.675}, n,std::vector<double>{1,.59160797*5,0,2.675});
// U[mesh->getBCMap(3 * i + 0)] = flux_farfield(U[mesh->getBCMap(3 * i + 1)], n, U[mesh->getBCMap(3 * i + 1)]);
U[mesh->getBCMap(3 * i + 0)] = flux_farfield(U[mesh->getBCMap(3 * i + 1)], n, U_out);
// U[mesh->getBCMap(3 * i + 0)] = flux_outflow(U[mesh->getBCMap(3 * i + 1)], n,p_in);
// U[mesh->getBCMap(3 * i + 0)] = U[mesh->getBCMap(3 * i + 1)];
break;
}
}
}
}
int nBC2 = mesh->getNbc2()/2;
int dimSys = mesh->getNelem();
std::vector<double> tmp;
std::vector<double> tmp2;
for (int i = 0; i < nBC2; i++){
tmp = U[mesh->getBCMap2(3*i + 1)];
tmp2 = U[mesh->getBCMap2(3*(i + nBC2) + 1)];
for (int j = 0; j < dimSys; j++){
U[mesh->getBCMap2(3*i + 0)][j] = 0.5*(tmp[j] + tmp2[j]);
}
}
}
void NozzleBC::updateEdges(std::vector<std::vector<std::vector<std::vector<double>>>> &UReconstr, MeshBase *mesh, double t){
std::vector<std::vector<int>> edgesAtBC = mesh->getEdgesAtBC();
int nelem = mesh->getNelem();
int nQuad = UReconstr[0][0].size();
int nEdgesBC = edgesAtBC.size();
std::vector<double> n;
for(int i = 0; i < nEdgesBC; i++) {
// if (edgesAtBC[i][0] == 2644){
// std::cout << i ;
// }
nQuad = UReconstr[0][edgesAtBC[i][0]].size();
for (int k = 0; k < nQuad; k++){
// if ((edgesAtBC[i][1] == 1)||(edgesAtBC[i][1] == 3)){// top and bottom wall with slipwall BC
if ((edgesAtBC[i][1] == 2)){// slipwall BC
if (edgesAtBC[i][2] == 0) { //means right side of edge is out & normal is OUTWARDS pointing
n = mesh->getNormal(edgesAtBC[i][0]);
flux_slipWall(UReconstr[1][edgesAtBC[i][0]][k], n, UReconstr[0][edgesAtBC[i][0]][k]);
}else{//means left side of edge is out & normal is INWARDS pointing
n = mesh->getNormal(edgesAtBC[i][0]);
n[0] = -n[0];
n[1] = -n[1];
flux_slipWall(UReconstr[0][edgesAtBC[i][0]][k], n, UReconstr[1][edgesAtBC[i][0]][k]);
}
}else if ((edgesAtBC[i][1] == 0)){ // means inflow BC
if (edgesAtBC[i][2] == 0) { //means right side of edge is out & normal is OUTWARDS pointing
n = mesh->getNormal(edgesAtBC[i][0]);
UReconstr[0][edgesAtBC[i][0]][k] = flux_farfield(UReconstr[1][edgesAtBC[i][0]][k], n, U_in);
}else{//means left side of edge is out & normal is INWARDS pointing
n = mesh->getNormal(edgesAtBC[i][0]);
n[0] = -n[0];
n[1] = -n[1];
UReconstr[1][edgesAtBC[i][0]][k] = flux_farfield(UReconstr[0][edgesAtBC[i][0]][k], n, U_in);
}
// }else{ //right wall, means outflow BC
}else if ((edgesAtBC[i][1] == 1)){ // means outflow BC
if (edgesAtBC[i][2] == 0) { //means right side of edge is out & normal is OUTWARDS pointing
n = mesh->getNormal(edgesAtBC[i][0]);
// UReconstr[0][edgesAtBC[i][0]][k] = UReconstr[1][edgesAtBC[i][0]][k];
UReconstr[0][edgesAtBC[i][0]][k] = flux_farfield(UReconstr[1][edgesAtBC[i][0]][k], n, U_out);
// UReconstr[0][edgesAtBC[i][0]][k] = flux_farfield(UReconstr[1][edgesAtBC[i][0]][k], n, UReconstr[1][edgesAtBC[i][0]][k]);
// UReconstr[0][edgesAtBC[i][0]][k] = flux_outflow(UReconstr[1][edgesAtBC[i][0]][k], n, p_in);
}else{//means left side of edge is out & normal is INWARDS pointing
n = mesh->getNormal(edgesAtBC[i][0]);
n[0] = -n[0];
n[1] = -n[1];
// UReconstr[1][edgesAtBC[i][0]][k] = UReconstr[0][edgesAtBC[i][0]][k];
UReconstr[1][edgesAtBC[i][0]][k] = flux_farfield(UReconstr[0][edgesAtBC[i][0]][k], n, U_out);
// UReconstr[1][edgesAtBC[i][0]][k] = flux_farfield(UReconstr[0][edgesAtBC[i][0]][k], n, UReconstr[0][edgesAtBC[i][0]][k]);
// UReconstr[1][edgesAtBC[i][0]][k] = flux_outflow(UReconstr[0][edgesAtBC[i][0]][k], n, p_in);
}
}
}
if ((edgesAtBC[i][1] == 11)){ // means hybrid Connection BC
if (edgesAtBC[i][2] == 0) { //means right side of edge is out & normal is OUTWARDS pointing
for (int k = 0; k < nQuad; k++){
UReconstr[0][edgesAtBC[i][0]][k] = UReconstr[0][edgesAtBC[i][3]][0];
}
}else{//means left side of edge is out & normal is INWARDS pointing
for ( int l = 0; l < nelem; l++){
UReconstr[1][edgesAtBC[i][0]][0][l] = 0;
for (int k = 0; k < nQuad; k++){
UReconstr[1][edgesAtBC[i][0]][0][l] += UReconstr[1][edgesAtBC[i][3]][k][l]/nQuad;
}
}
}
}else if ((edgesAtBC[i][1] == 14)){ // means hybrid QUADQUAD 2 QUAD or QUAD 2 QUADQUAD Connection BC
if (edgesAtBC[i][2] == 0) { //means right side of edge is out & normal is OUTWARDS pointing
UReconstr[0][edgesAtBC[i][0]]= UReconstr[1][edgesAtBC[i][3]];
}else{//means left side of edge is out & normal is INWARDS pointing
UReconstr[1][edgesAtBC[i][0]]= UReconstr[0][edgesAtBC[i][3]];
}
}
}
}
void NozzleBC::updateEdges(std::vector<std::vector<std::vector<std::vector<double>>>> &UReconstr, MeshBase *mesh, double t, int i){
std::vector<std::vector<int>> edgesAtBC = mesh->getEdgesAtBC();
int nQuad = UReconstr[0][0].size();
int nEdgesBC = edgesAtBC.size();
std::vector<double> n;
for (int k = 0; k < nQuad; k++){
if ((edgesAtBC[i][1] == 2)){// slipwall BC
if (edgesAtBC[i][2] == 0) { //means right side of edge is out & normal is OUTWARDS pointing
n = mesh->getNormal(edgesAtBC[i][0]);
flux_slipWall(UReconstr[1][edgesAtBC[i][0]][k], n, UReconstr[0][edgesAtBC[i][0]][k]);
}else{//means left side of edge is out & normal is INWARDS pointing
n = mesh->getNormal(edgesAtBC[i][0]);
n[0] = -n[0];
n[1] = -n[1];
flux_slipWall(UReconstr[0][edgesAtBC[i][0]][k], n, UReconstr[1][edgesAtBC[i][0]][k]);
}
}else if ((edgesAtBC[i][1] == 0)){ // means inflow BC
if (edgesAtBC[i][2] == 0) { //means right side of edge is out & normal is OUTWARDS pointing
n = mesh->getNormal(edgesAtBC[i][0]);
UReconstr[0][edgesAtBC[i][0]][k] = flux_farfield(UReconstr[1][edgesAtBC[i][0]][k], n, U_in);
}else{//means left side of edge is out & normal is INWARDS pointing
n = mesh->getNormal(edgesAtBC[i][0]);
n[0] = -n[0];
n[1] = -n[1];
UReconstr[1][edgesAtBC[i][0]][k] = flux_farfield(UReconstr[0][edgesAtBC[i][0]][k], n, U_in);
}
}else if ((edgesAtBC[i][1] == 1)){ // means outflow BC
if (edgesAtBC[i][2] == 0) { //means right side of edge is out & normal is OUTWARDS pointing
n = mesh->getNormal(edgesAtBC[i][0]);
// UReconstr[0][edgesAtBC[i][0]][k] = UReconstr[1][edgesAtBC[i][0]][k];
UReconstr[0][edgesAtBC[i][0]][k] = flux_farfield(UReconstr[1][edgesAtBC[i][0]][k], n, U_out);
// UReconstr[0][edgesAtBC[i][0]][k] = flux_farfield(UReconstr[1][edgesAtBC[i][0]][k], n, UReconstr[1][edgesAtBC[i][0]][k]);
// UReconstr[0][edgesAtBC[i][0]][k] = flux_outflow(UReconstr[1][edgesAtBC[i][0]][k], n, p_in);
}else{//means left side of edge is out & normal is INWARDS pointing
n = mesh->getNormal(edgesAtBC[i][0]);
n[0] = -n[0];
n[1] = -n[1];
// UReconstr[1][edgesAtBC[i][0]][k] = UReconstr[0][edgesAtBC[i][0]][k];
UReconstr[1][edgesAtBC[i][0]][k] = flux_farfield(UReconstr[0][edgesAtBC[i][0]][k], n, U_out);
// UReconstr[1][edgesAtBC[i][0]][k] = flux_farfield(UReconstr[0][edgesAtBC[i][0]][k], n, UReconstr[0][edgesAtBC[i][0]][k]);
// UReconstr[1][edgesAtBC[i][0]][k] = flux_outflow(UReconstr[0][edgesAtBC[i][0]][k], n, p_in);
}
}
}
}
void NozzleBC::flux_slipWall(std::vector<double> U, std::vector<double> n, std::vector<double> &ret){
ret = U;
std::vector<double> V_in{U[1]/U[0],U[2]/U[0]};
double vn = V_in[0]*n[0] + V_in[1]*n[1];
ret[1] = U[0]*(V_in[0]-2*vn*n[0]);
ret[2] = U[0]*(V_in[1]-2*vn*n[1]);
}
Event Timeline
Log In to Comment