//----------------------------------------------------------------------------- // NastNavierStokes2dMAC.cpp //----------------------------------------------------------------------------- // // Copyright (C) 1998 Technische Universitaet Muenchen, Germany // written by Bernhard Brueck // // This file is part of Nast++ // //----------------------------------------------------------------------------- // CNastNavierStokes2dMAC implementiert den Loeser fuer die // Navier-Stokes-Gleichungen mit dem Marker-Cell-Verfahren. // (nicht verwechseln mit MAC fuer die Geometrie) // Dazu wird ein vorlaeufiges Geschwindigkeitsfeld berechnet und danach // Druck angepasst. //----------------------------------------------------------------------------- // Changes: // #include "NastNavierStokes2dMAC.h" #include "NastDebug.h" #include "NastGeometry2d.h" #include "NastSolverPoisson2d.h" #include "NastDiffQuot2d.h" #include "NastOutputContextIDL.h" #include // Workaround fuer einen Bug in HP-Includefiles #include // DBL_MAX #ifdef _MSC_VER # include #endif //----------------------------------------------------------------------------- // Konstruktor + Destruktor //----------------------------------------------------------------------------- // geometry Geometriebeschreibung // solver Loeser fue die Druckgleichungen // diffQuot die Differenzenquotienten // nSizeI Anzahl der Zellen in X-Richtung // nSizeJ Anzahl der Zellen in Y-Richtung // CNastNavierStokes2dMAC::CNastNavierStokes2dMAC(CNastSolverPoisson2d &solver, CNastDiffQuot2d &diffQuot, const CNastGeometry2d &geom, const CNastNavierStokesParameter para, int nSizeI, int nSizeJ) :m_grid( geom, geom.boundingBox(), nSizeI, nSizeJ ), m_gridF( m_grid.gridU().box(), m_grid.gridU().sizeI(), m_grid.gridU().sizeJ()), m_gridG( m_grid.gridV().box(), m_grid.gridV().sizeI(), m_grid.gridV().sizeJ()), m_gridRhs( m_grid.gridP().box(), m_grid.gridP().sizeI(), m_grid.gridP().sizeJ()), m_para ( para ), m_solver ( solver ), m_geom ( geom ), m_diffQuot( diffQuot ) { // den Differenzenquotienten noch das aktuelle Gitter mitteilen m_diffQuot.useGrid( m_grid ); // for( int j = m_grid.sizeJ() / 2+1; j < m_grid.sizeJ()+1; j++) // for( int i = 2; i < m_grid.sizeI()+1; i++) // { // NAST_DUMP_VARIABLE(i); // NAST_DUMP_VARIABLE(j); // if( m_grid.cell(i,j).isFluid( CNastCell2d::THIS | CNastCell2d::WEST )) // m_grid.setU( i, j, 1); // } } // Destruktor CNastNavierStokes2dMAC::~CNastNavierStokes2dMAC() { NAST_ASSERT_VALID( this ); } //------------------------------------------------------------------------- // Memberfunktionen //------------------------------------------------------------------------- static double nastMin( double x, double y) { if( x < y ) return x; return y; } // Berechnung der maximal moeglichen Schrittweite double CNastNavierStokes2dMAC::maxTimeStep() const { NAST_ASSERT_VALID( this ); double dx = m_grid.dx(); double dy = m_grid.dy(); double maxU = m_grid.gridU().absMax(); double maxV = m_grid.gridV().absMax(); double deltU = DBL_MAX; double deltV = DBL_MAX; if( maxU ) deltU = dx / maxU; if( maxV ) deltV = dy / maxV; double deltRe = 0.5 * m_para.reynolds() / (1 / (dx * dx) + 1/(dy * dy)); // der Sicherheitsfaktor tau ist in CNastSimulator2d::calcTimeStep return nastMin(nastMin( deltU, deltV), deltRe); } void CNastNavierStokes2dMAC::solve( double t ) { NAST_ASSERT_VALID( this ); // bei Zeitschrittweitensteuerung darf diese nicht // ueberschritten werden double dt = t - m_grid.time(); NAST_ASSERT_WARNING( dt < maxTimeStep(),"zu grosse Schrittweite" ); NAST_ASSERT( dt > 0 ); // vorlaeufiges Geschwindigkeitsfeld berechnen compFG( dt ); // die rechte Seite fuer die Druckgleichungen berechnen compRhs( dt ); // die Druckgleichungen loesen m_solver.solve( m_grid, m_gridRhs ); // die Geschwindigkeitswerte anpassen adapUV( dt ); } const CNastStaggeredGrid2d& CNastNavierStokes2dMAC::grid() const { NAST_ASSERT_VALID( this ); return m_grid; } //------------------------------------------------------------------------- // Hilfsfunktionen //------------------------------------------------------------------------- // Berechnung des vorlaeufigen Geschwindigkeitsfeldes fuer // den Zeitschritt dt void CNastNavierStokes2dMAC::compFG( double dt ) { double rre = 1 / m_para.reynolds(); for( int i = 1; i <= m_grid.sizeI()+1; i++) for( int j = 1; j <= m_grid.sizeJ()+1; j++) if( m_grid.cell(i,j).isFluid() ) { m_gridF(i , j) = m_grid.U(i, j); m_gridF(i+1, j) = m_grid.U(i+1, j); m_gridG(i, j ) = m_grid.V(i, j ); m_gridG(i, j+1) = m_grid.V(i, j+1); if( m_grid.cell(i,j).isFluid( CNastCell2d::THIS | CNastCell2d::WEST )) // Kante zwischen Fluidzellen ? { const double du2_dx = m_diffQuot.diff_du2_dx( i, j ); const double duv_dy = m_diffQuot.diff_duv_dy( i, j ); const double lapl_u = m_diffQuot.diff_lapl_u( i, j ); const double u = m_grid.U(i,j); m_gridF(i,j) = u + dt * (lapl_u * rre - du2_dx - duv_dy + m_para.gx()); } if( m_grid.cell(i,j).isFluid( CNastCell2d::THIS | CNastCell2d::SOUTH )) // Kante zwischen Fluidzellen ? { const double dv2_dy = m_diffQuot.diff_dv2_dy( i, j ); const double duv_dx = m_diffQuot.diff_duv_dx( i, j ); const double lapl_v = m_diffQuot.diff_lapl_v( i, j ); const double v = m_grid.V(i,j); m_gridG(i,j) = v + dt * (lapl_v * rre - dv2_dy - duv_dx + m_para.gy()); } } // zum Debuggen alles anzeigen if(0) { static int nCount = 0; if( nCount++ % 1 == 0) { int nSizeI = m_grid.gridP().sizeI(); int nSizeJ = m_grid.gridP().sizeJ(); CNastBox2d box = m_grid.gridP().box(); // zum debuggen alle diffquot anzeigen CNastGrid2d grid_du2_dx( box, nSizeI, nSizeJ ); CNastGrid2d grid_duv_dy( box, nSizeI, nSizeJ ); CNastGrid2d grid_lapl_u( box, nSizeI, nSizeJ ); CNastGrid2d grid_dv2_dy( box, nSizeI, nSizeJ ); CNastGrid2d grid_duv_dx( box, nSizeI, nSizeJ ); CNastGrid2d grid_lapl_v( box, nSizeI, nSizeJ ); for( int i = 1; i <= m_grid.sizeI()+1; i++) for( int j = 1; j <= m_grid.sizeJ()+1; j++) if( m_grid.cell(i,j).isFluid() ) { if( m_grid.cell(i,j).isFluid( CNastCell2d::THIS | CNastCell2d::WEST )) // Kante zwischen Fluidzellen ? { grid_du2_dx(i,j) = m_diffQuot.diff_du2_dx( i, j ); grid_duv_dy(i,j) = m_diffQuot.diff_duv_dy( i, j ); grid_lapl_u(i,j) = m_diffQuot.diff_lapl_u( i, j ); } else m_gridF(i,j) = m_grid.U(i,j); if( m_grid.cell(i,j).isFluid( CNastCell2d::THIS | CNastCell2d::SOUTH )) // Kante zwischen Fluidzellen ? { grid_dv2_dy(i,j) = m_diffQuot.diff_dv2_dy( i, j ); grid_duv_dx(i,j) = m_diffQuot.diff_duv_dx( i, j ); grid_lapl_v(i,j) = m_diffQuot.diff_lapl_v( i, j ); } else m_gridG(i,j) = m_grid.V(i,j); } // NAST_DUMP_VARIABLE(m_grid); NAST_DUMP_VISUAL(m_grid); NAST_DUMP_VISUAL(m_grid.gridU()); NAST_DUMP_VISUAL(m_grid.gridV()); NAST_DUMP_VISUAL(m_grid.gridP()); // NAST_DUMP_VISUAL(m_gridF); // NAST_DUMP_VISUAL(m_gridG); // NAST_DUMP_VISUAL( grid_lapl_u ); // NAST_DUMP_VISUAL( grid_lapl_v ); // NAST_DUMP_VISUAL( grid_du2_dx ); // NAST_DUMP_VISUAL( grid_dv2_dy ); // NAST_DUMP_VISUAL( grid_duv_dy ); // NAST_DUMP_VISUAL( grid_duv_dx ); NAST_BREAK(); } } } // rechte Seite fuer die Druckkorrektur berechnen void CNastNavierStokes2dMAC::compRhs( double dt) { const double rdxdt = 1 / (m_grid.dx() * dt); const double rdydt = 1 / (m_grid.dy() * dt); for ( int i=1; i <= m_grid.sizeI()+1; i++) for ( int j=1; j <= m_grid.sizeJ()+1; j++) { if( m_grid.cell(i,j).isFluid() ) // Zelle ist Fluidzelle ? { m_gridRhs(i, j) = (m_gridF( i+1, j ) - m_gridF( i, j )) * rdxdt + (m_gridG( i , j+1 ) - m_gridG( i, j )) * rdxdt; } else { m_gridRhs(i, j) = 0; } } // NAST_DUMP_VISUAL( m_gridRhs ); } // neue Geschwindigkeit berechnen void CNastNavierStokes2dMAC::adapUV( double dt ) { int i,j; const double dt_dx = dt / m_grid.dx(); const double dt_dy = dt / m_grid.dy(); for( i = 1; i <= m_grid.sizeI()+1; i++ ) for( j = 1; j <= m_grid.sizeJ()+1; j++ ) { CNastCell2d cell = m_grid.cell(i,j); // nur bei Kanten zwischen zwei Fluidzellen if( cell.isFluid( CNastCell2d::THIS | CNastCell2d::WEST )) { const double diff = m_grid.gridP()(i-1, j) - m_grid.gridP()(i, j); const double u = m_gridF(i,j) + diff * dt_dx; // und den Wert setzen m_grid.setU( i, j, u ); } // nur bei Kanten zwischen zwei Fluidzellen if( cell.isFluid( CNastCell2d::THIS | CNastCell2d::SOUTH )) { const double diff = m_grid.gridP()(i, j-1) - m_grid.gridP()(i, j); const double v = m_gridG(i,j) + diff * dt_dy; // und den Wert setzen m_grid.setV( i, j, v ); } } // und die Zeit im Gitter weitersetzen m_grid.time( dt + m_grid.time()); } // ---------------------------------------------------------------------------- // Debug // ---------------------------------------------------------------------------- void CNastNavierStokes2dMAC::debugDump( CNastDumpContext &dumpContext ) const { CNastNavierStokes2d::debugDump( dumpContext ); dumpContext << "CNastNavierStokes2dMAC" << "\n"; dumpContext << "\tm_solver = " << m_solver << "\n"; dumpContext << "\tm_diffQuot = " << m_diffQuot << "\n"; dumpContext << "\tm_grid = " << m_grid << "\n"; dumpContext << "\tm_gridF = " << m_gridF << "\n"; dumpContext << "\tm_gridG = " << m_gridG << "\n"; } void CNastNavierStokes2dMAC::assertValid() const { // zuerst einmal AssertValid der Basisklasse aufrufen CNastNavierStokes2d::assertValid(); }