diff --git a/source/Applications/pChannel/configBombadilpChannel.cfg b/source/Applications/pChannel/configBombadilpChannel.cfg index c1ae0bf884dfc269a8abc343afb9648e5bb0ec75..e7ba220f45dc3b6321eefc37f88fed87f60fb23e 100644 --- a/source/Applications/pChannel/configBombadilpChannel.cfg +++ b/source/Applications/pChannel/configBombadilpChannel.cfg @@ -2,10 +2,10 @@ #Simulation parameters for porous channel # -pathname = d:/temp/pChannel6 +pathname = d:/temp/pChannel pathGeo = d:/Projects/SFB880/GeometrienPoroeseMedien/PA80-110 numOfThreads = 4 -availMem = 4e9 +availMem = 8e9 logToFile = false #porous media @@ -53,21 +53,11 @@ pmL = 1e-3 1e-3 1e-3 #grid refineLevel = 1 -#refineLevel = 1 -#refineLevel = 0 -#deltaXfine = 3.6496350365e-5 deltaXfine = 20e-6 -# x10 -#deltaXfine = 1.4598540146e-5# x4 -#3deltaXfine = 3.6496350365e-6 -#deltaXfine = 3.95e-6 -#deltaXfine = 8.0e-6 blocknx = 10 10 10 -#blocknx = 32 40 20 lengthFactor = 4 -thinWall = false -forcing = 7.83841e-09 -changeQs = false +thinWall = true +changeQs = true channelHigh = 0.002 @@ -82,15 +72,15 @@ channelHigh = 0.002 # for DLRF-15 Re = 102000/2 Re = 51000 #real velocity is 54.95 m/s -u_LB = 0.01 +u_LB = 0.1 -restartStep = 500 -restartStepStart = 500 +restartStep = 20000 +restartStepStart = 20000 -timeAvStart = 1 -timeAvStop = 5 +timeAvStart = 10000 +timeAvStop = 20000 -endTime = 1000 -outTime = 1000 +endTime = 20000 +outTime = 10000 diff --git a/source/Applications/pChannel/pChannel.cpp b/source/Applications/pChannel/pChannel.cpp index 6552f986663721e3495e9c6d7605763a0f4ad0ff..446f8f9f08ff3a779916e345199e6bcd44f912aa 100644 --- a/source/Applications/pChannel/pChannel.cpp +++ b/source/Applications/pChannel/pChannel.cpp @@ -41,7 +41,6 @@ void run(string configname) double Re = config.getDouble("Re"); double channelHigh = config.getDouble("channelHigh"); double lengthFactor = config.getDouble("lengthFactor"); - double forcing = config.getDouble("forcing"); bool changeQs = config.getBool("changeQs"); double timeAvStart = config.getDouble("timeAvStart"); double timeAvStop = config.getDouble("timeAvStop"); @@ -107,7 +106,7 @@ void run(string configname) double offsetMaxX1 = pmL[0]*lengthFactor; double offsetMaxX2 = pmL[1]*2.0; - double offsetMaxX3 = pmL[2] + channelHigh; //DLR-F15 //pmL[2]*2.0; + double offsetMaxX3 = channelHigh; // pmL[2] + channelHigh; //DLR-F15 //pmL[2]*2.0; //bounding box double g_minX1 = origin[0]; @@ -295,6 +294,7 @@ void run(string configname) intHelper.setBC(); ////porous media + if(false) { string samplePathname = pathGeo + "/" + sampleFilename; @@ -372,20 +372,59 @@ void run(string configname) grid->accept(bcVisitor); mu::Parser inflowProfile; - //inflowProfile.SetExpr("x3 < h ? 0.0 : uLB+1*x1-1*x2"); - //inflowProfile.SetExpr("uLB+1*x1-1*x2"); - inflowProfile.SetExpr("uLB"); - //inflowProfile.DefineConst("uLB", u_LB); - inflowProfile.DefineConst("uLB", 0.0116); + // inflowProfile.SetExpr("x3 < h ? 0.0 : uLB+1*x1-1*x2"); + ////inflowProfile.SetExpr("uLB+1*x1-1*x2"); + // //inflowProfile.SetExpr("uLB"); + // inflowProfile.DefineConst("uLB", u_LB); + // //inflowProfile.DefineConst("uLB", 0.0116); + // inflowProfile.DefineConst("h", pmL[2]); + + //inflowProfile = Utilities::getDuctParaboloidX(g_maxX2/2.0, g_maxX2, pmL[2]+channelHigh/2.0, channelHigh, 0.1); + + //poiseuille flow + //double Cy = g_maxX2 / 2.0; + //double Hy = g_maxX2; + //double Cz = pmL[2] + channelHigh / 2.0; + //double Hz = channelHigh; + //double V = (9.0/4.0)*u_LB; + //inflowProfile.SetExpr("x3 < h ? 0.0 : V*(((-(x2-Cy)^2.0+(Hy/2.0)^2.0)/(Hy/2.0)^2.0)*((-(x3-Cz)^2.0+(Hz/2.0)^2.0)/(Hz/2.0)^2.0))"); + //inflowProfile.DefineConst("Cy", Cy); + //inflowProfile.DefineConst("Hy", Hy); + //inflowProfile.DefineConst("Cz", Cz); + //inflowProfile.DefineConst("Hz", Hz); + //inflowProfile.DefineConst("V", V); //inflowProfile.DefineConst("h", pmL[2]); - //inflowProfile = Utilities::getDuctParaboloidX(g_maxX2/2.0, g_maxX2, pmL[2]+channelHigh/2.0, channelHigh, 0.116); + //log-law + double u_tau = 0.05;//sqrt(nu_LB * (u_LB/channel_high)); //0.0013252866627413104; + double z0 = 0.0001;//nu_LB / (9.0*u_tau); + double k = 0.4; + //inflowProfile.SetExpr("2.3*u_tau/k*log(x3/z0)"); + //inflowProfile.SetExpr("x3 > 0 && (zMax-x3) > 0 ? (x3 < h ? 2.3*u_tau/k*log(x3/z0) : 2.3*u_tau/k*log((zMax-x3)/z0)) : 0"); + //inflowProfile.SetExpr("x3 > 0 && (zMax-x3) > 0 ? (x3 < h ? (1.0/k)*log(9.8*x3/u_tau)*u_tau : (1.0/k)*log(9.8*(zMax-x3)/u_tau)*u_tau ) : 0"); + //inflowProfile.SetExpr("x3 < h && x3 > 0 ? 2.3*u_tau/k*log(x3/z0) : 0.0"); + + //inflowProfile.SetExpr("Uref/log((Href+z0)/z0)*log((x3-zg+z0)/z0)"); + //inflowProfile.SetExpr("x3 > 0 && (zMax-x3) > 0 ? (x3 < h ? Uref/log((Href+z0)/z0)*log((x3-zg+z0)/z0) : Uref/log((Href+z0)/z0)*log((zMax-x3-zg+z0)/z0)) : 0"); + inflowProfile.SetExpr("x3 < h ? Uref/log((Href+z0)/z0)*log((x3-zg+z0)/z0) : Uref/log((Href+z0)/z0)*log((zMax-x3-zg+z0)/z0)"); + + inflowProfile.DefineConst("Uref", u_LB); + inflowProfile.DefineConst("Href", channelHigh); + inflowProfile.DefineConst("zg", 0.0); + + inflowProfile.DefineConst("u_tau", u_tau); + inflowProfile.DefineConst("k", k); + inflowProfile.DefineConst("z0", z0); + inflowProfile.DefineConst("h", channelHigh / 2.0); + inflowProfile.DefineConst("zMax", channelHigh); + D3Q27ETInitDistributionsBlockVisitor initVisitor(nu_LB, rho_LB); initVisitor.setVx1(inflowProfile); //initVisitor.setVx1(u_LB); grid->accept(initVisitor); + ////set connectors D3Q27InterpolationProcessorPtr iProcessor(new D3Q27IncompressibleOffsetInterpolationProcessor()); D3Q27SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nu_LB, iProcessor); @@ -493,14 +532,14 @@ void run(string configname) UbSchedulerPtr AdjForcSch(new UbScheduler()); - AdjForcSch->addSchedule(10, 0, 10000000); + AdjForcSch->addSchedule(1, 0, 10000000); D3Q27IntegrateValuesHelperPtr intValHelp(new D3Q27IntegrateValuesHelper(grid, comm, coord[0], coord[1], coord[2], coord[3], coord[4], coord[5])); if (myid == 0) GbSystem3D::writeGeoObject(intValHelp->getBoundingBox().get(), pathname + "/geo/IntValHelp", WbWriterVtkXmlBinary::getInstance()); double vxTarget=u_LB; - AdjustForcingCoProcessor AdjForcPPPtr(grid, AdjForcSch, pathname, intValHelp, vxTarget, forcing, comm); + AdjustForcingCoProcessor AdjForcPPPtr(grid, AdjForcSch, pathname, intValHelp, vxTarget, comm); //mu::Parser decrViscFunc; //decrViscFunc.SetExpr("nue0+c0/(t+1)/(t+1)"); @@ -510,15 +549,15 @@ void run(string configname) //DecrViscSch->addSchedule(10, 0, 1000); //DecreaseViscosityCoProcessor decrViscPPPtr(grid, DecrViscSch, &decrViscFunc, comm); - //if (changeQs) - //{ - // double z1 = pmL[2]; - // D3Q27IntegrateValuesHelperPtr intValHelp2(new D3Q27IntegrateValuesHelper(grid, comm, - // coord[0], coord[1], z1 - deltaXfine, - // coord[3], coord[4], z1 + deltaXfine)); - // if (myid == 0) GbSystem3D::writeGeoObject(intValHelp2->getBoundingBox().get(), pathname + "/geo/intValHelp2", WbWriterVtkXmlBinary::getInstance()); - // Utilities::ChangeRandomQs(intValHelp2); - //} + if (changeQs) + { + double z1 = pmL[2]; + D3Q27IntegrateValuesHelperPtr intValHelp2(new D3Q27IntegrateValuesHelper(grid, comm, + coord[0], coord[1], z1 - deltaXfine, + coord[3], coord[4], z1 + deltaXfine)); + if (myid == 0) GbSystem3D::writeGeoObject(intValHelp2->getBoundingBox().get(), pathname + "/geo/intValHelp2", WbWriterVtkXmlBinary::getInstance()); + Utilities::ChangeRandomQs(intValHelp2); + } std::vector<double> levelCoords; std::vector<int> levels; @@ -537,7 +576,7 @@ void run(string configname) levelCoords.push_back(0.0024); levelCoords.push_back(0.003); UbSchedulerPtr tavSch(new UbScheduler(1, timeAvStart, timeAvStop)); - TimeAveragedValuesCoProcessorPtr tav(new TimeAveragedValuesCoProcessor(grid, pathname, WbWriterVtkXmlBinary::getInstance(), tavSch, + TimeAveragedValuesCoProcessorPtr tav(new TimeAveragedValuesCoProcessor(grid, pathname, WbWriterVtkXmlBinary::getInstance(), tavSch, comm, TimeAveragedValuesCoProcessor::Velocity | TimeAveragedValuesCoProcessor::Fluctuations | TimeAveragedValuesCoProcessor::Triplecorrelations, levels, levelCoords, bounds));