//----------------------------------------------------------------------------- // NastStaggeredGrid2d.cpp //------------------------- ---------------------------------------------------- // // Copyright (C) 1998 Technische Universitaet Muenchen, Germany // written by Bernhard Brueck // // This file is part of Nast++ // //----------------------------------------------------------------------------- // CNastStaggered Grid enthaelt die Werte fuer Druck und Geschwindigkeit // und ermoeglicht den Zugriff auf die geometrischen Daten. // //----------------------------------------------------------------------------- // Changes: // #include "NastStaggeredGrid2d.h" #include "NastDebug.h" #include "NastBoundaryPoint2d.h" #include "NastGeometry2d.h" //------------------------------------------------------------------------- // Diskretisierung der Werte //------------------------------------------------------------------------- // | | | | // --+---------------------+---------------------+---------------------+-- // | | | | // | | | | // | | | | // | | | | // U(i-1,j+1) U(i-1,j+1) | | // | | | | // | | | | // | | | | // | | | | // --+---------------------+-------V(i,j+1)------+---------------------+-- // | | | | // | | | | // | | | | // | | | | // U(i-1,j) U(i,j) i,j U(i+1,j) | // | | | | // | | | | // | | | | // | | | | // --+-------V(i-1,j)------+--------V(i,j)-------+-------V(i+1,j)------+-- // | | | | // | | | | // | | | | // | | | | // | U(i-1,j) U(i,j-1) | // | | | | // | | | | // | | | | // | | | | // --+---------------------+-------V(i,j-1)------+------V(i+1,j-1)-----+-- // | | | | // ^ // | // Y CNastCell2d::NORTH // | // ----+--X--> CNastCell2d::WEST CNastCell2d::THIS CNastCell2d::EAST // | // | CNastCell2d::SOUTH //----------------------------------------------------------------------------- // Konstruktor + Destruktor //----------------------------------------------------------------------------- // Konstruktion eines versetzten Gitters unter Angabe der Ausdehnung, der // Geometrie und der Anzahl der Zellen // Achtung die Angaben box, nSizeI, nSizeJ beziehen sich auf die Ausdehnung // der Fluidzellen im inneren. Es werden auch noch die Randzellen // benoetigt und damit hat das Gitter die nSizeI+2, nSizeJ+2 Zellen // und auch die Ausdehnung wird um eine zwei Zellenbreiten groesser. CNastStaggeredGrid2d::CNastStaggeredGrid2d( const CNastGeometry2d &geom, const CNastBox2d &box, int nSizeI, int nSizeJ ) :m_dx( box.sizeX() / (nSizeI)), m_dy( box.sizeY() / (nSizeJ)), m_sizeI( nSizeI ), m_sizeJ( nSizeJ ), m_gridU( box.pntMin() - CNastVector2d( m_dx, m_dy * 0.5 ), box.pntMax() + CNastVector2d( m_dx, m_dy * 0.5 ), nSizeI+3, nSizeJ+2 ), m_gridV( box.pntMin() - CNastVector2d( m_dx * 0.5, m_dy), box.pntMax() + CNastVector2d( m_dx * 0.5, m_dy), nSizeI+2, nSizeJ+3), m_gridP( box.pntMin() - CNastVector2d( m_dx * 0.5, m_dy * 0.5), box.pntMax() + CNastVector2d( m_dx * 0.5, m_dy * 0.5), nSizeI+2, nSizeJ+2 ), m_cells( CNastBox2d(box).expand( m_dx, m_dy ), nSizeI+2, nSizeJ+2 ), m_geom( geom ), m_box( box ), m_time(0) { // und gleich mal die Zellen richtig setzen updateAll(); NAST_DUMP_VARIABLE(m_cells); } // Destruktor CNastStaggeredGrid2d::~CNastStaggeredGrid2d() { NAST_ASSERT_VALID( this ); } //------------------------------------------------------------------------- // Hilfsfunktionen //------------------------------------------------------------------------- // Wenn die Werte ausserhalb des Fluids liegen, dann muessen sie in Ab- // haengikeit der jeweiligen Randbedingung bestimmt werden. Dazu wird // die Geometrie nach der gueltigen Randbedingung gefragt und der // Wert berechnet. //------------------------------------------------------------------------- //------------------------------------------------------------------------- // Den V Wert westlich von V(i,j) berechnen, wenn er // ausserhalb des Fluids liegt // // | | | | // --+---------------------+---------------------+---------------------+-- // | | | | // | | | | // | | i,j | | // | Zelle mit | | | // | Randbedingung | Fluidzelle | | // | | | | // | | | | // | | | | // | +----------+ | | | // --+-----| V(i-1,j) |----+--------V(i,j)-------+---------------------+-- // | +----------+ | | | // | | | | // | | | | // | Zelle mit | | | // | Randbedingung | Fluidzelle | | // | | | | // | | | | // | | | | // | | | | // --+---------------------+---------------------+---------------------+-- // | | | | // double CNastStaggeredGrid2d::calcWestV( int i, int j ) const { NAST_ASSERT( cell(i,j).isFluid( CNastCell2d::THIS | CNastCell2d::SOUTH )); // den Randpunkt von der Geometrie erfragen const CNastBoundaryPoint2d pnt( m_geom.boundaryPoint( m_gridV.pos(i,j), CNastVector2d(-1,0))); double v = 0; switch( pnt.boundaryCondition() ) { case CNastBoundaryPoint2d::SLIP: // Rutschbedingung case CNastBoundaryPoint2d::OUTFLOW: // Ausstroembedingung v = m_gridV(i,j); break; case CNastBoundaryPoint2d::NOSLIP: // Haftbedingung v = -m_gridV(i,j); break; case CNastBoundaryPoint2d::INFLOW: // Einstroembedingung v = 2 * pnt.velocity().y() - m_gridV(i,j); break; default: NAST_ASSERT( false ); // unbekannte Randbedingung } return v; } //------------------------------------------------------------------------- // Den V Wert oestlich von V(i,j) berechnen, wenn er // ausserhalb des Fluids liegt // // | | | | // --+---------------------+---------------------+---------------------+-- // | | | | // | | | | // | | i,j | | // | | | Zelle mit | // | | Fluidzelle | Randbedingung | // | | | | // | | | | // | | | | // | | | +----------+ | // --+---------------------+--------V(i,j)-------+-----| V(i+1,j) |----+-- // | | | +----------+ | // | | | | // | | | | // | | Fluidzelle | Zelle mit | // | | | Randbedingung | // | | | | // | | | | // | | | | // | | | | // --+---------------------+---------------------+---------------------+-- // | | | | // double CNastStaggeredGrid2d::calcEastV( int i, int j ) const { NAST_ASSERT( cell(i,j).isFluid( CNastCell2d::THIS | CNastCell2d::SOUTH )); // den Randpunkt von der Geometrie erfragen const CNastBoundaryPoint2d pnt( m_geom.boundaryPoint( m_gridV.pos(i, j), CNastVector2d(+1, 0))); double V = 0; switch( pnt.boundaryCondition() ) { case CNastBoundaryPoint2d::SLIP: // Rutschbedingung case CNastBoundaryPoint2d::OUTFLOW: // Ausstroembedingung V = m_gridV(i,j); break; case CNastBoundaryPoint2d::NOSLIP: // Haftbedingung V = -m_gridV(i,j); break; case CNastBoundaryPoint2d::INFLOW: // Einstroembedingung V = 2 * pnt.velocity().y() - m_gridV(i,j); break; default: NAST_ASSERT( false ); // unbekannte Randbedingung } return V; } //------------------------------------------------------------------------- // Den U Wert noerdlich von U(i,j) berechnen, wenn er // ausserhalb des Fluids liegt // | | | | // --+---------------------+--------------------------------------------- // | | | | // | | Zelle mit | Zelle mit | // | | Randbedingung | Randbedingung | // | | +----------+ | // | | | U(i,j+1) | | // | | +----------+ | // | | | | // | | | | // | | | | // --+---------------------+--------------------------------------------- // | | | | // | | | | // | | Fluid | Fluid | // | | | | // | | U(i,j) i,j | // | | | | // | | | | // | | | | // | | | | // --+---------------------+--------------------------------------------- // | | | | // double CNastStaggeredGrid2d::calcNorthU( int i, int j ) const { NAST_ASSERT( cell(i,j).isFluid( CNastCell2d::THIS | CNastCell2d::WEST )); // den Randpunkt von der Geometrie erfragen const CNastBoundaryPoint2d pnt( m_geom.boundaryPoint( m_gridU.pos(i,j), CNastVector2d(0,+1))); double u = 0; switch( pnt.boundaryCondition() ) { case CNastBoundaryPoint2d::SLIP: // Rutschbedingung case CNastBoundaryPoint2d::OUTFLOW: // Ausstroembedingung u = m_gridU(i,j); break; case CNastBoundaryPoint2d::NOSLIP: // Haftbedingung u = -m_gridU(i,j); break; case CNastBoundaryPoint2d::INFLOW: // Einstroembedingung u = 2 * pnt.velocity().x() - m_gridU(i,j); break; } return u; } // Den U Wert suedlich von U(i,j) berechnen, wenn er // ausserhalb des Fluids liegt // | | | | // --+---------------------+--------------------------------------------- // | | | | // | | | | // | | Fluid | Fluid | // | | | | // | | U(i,j) i,j | // | | | | // | | | | // | | | | // | | | | // --+---------------------+--------------------------------------------- // | | | | // | | | | // | | | | // | | +----------+ | // | | | U(i,j-1) | | // | | Zelle mit +----------+ Zelle mit | // | | Randbedingung | Randbedingung | // | | | | // | | | | // --+---------------------+--------------------------------------------- // | | | | // double CNastStaggeredGrid2d::calcSouthU( int i, int j ) const { NAST_ASSERT( cell(i,j).isFluid( CNastCell2d::THIS | CNastCell2d::WEST )); // den Randpunkt von der Geometrie erfragen const CNastBoundaryPoint2d pnt( m_geom.boundaryPoint( m_gridU.pos(i,j), CNastVector2d(0,-1))); double u = 0; switch( pnt.boundaryCondition() ) { case CNastBoundaryPoint2d::SLIP: // Rutschbedingung case CNastBoundaryPoint2d::OUTFLOW: // Ausstroembedingung u = m_gridU(i,j); break; case CNastBoundaryPoint2d::NOSLIP: // Haftbedingung u = -m_gridU(i,j); break; case CNastBoundaryPoint2d::INFLOW: // Einstroembedingung u = 2* pnt.velocity().x() - m_gridU(i,j); break; } return u; } // interpolierte Geschwindigkeit an einem bel. Punkt innerhalb // des Gitters erfragen CNastVector2d CNastStaggeredGrid2d::velocity( const CNastPoint2d &pnt ) const { NAST_ASSERT_VALID( this ); // funktioniert vorerst nur, wenn der Punkte innhalb der Gitter // fur U,V ist. fuer ausserhalb muesste erst eine Extrapolation // implementiert werden. const double u = m_gridU.interpolate(pnt); const double v = m_gridV.interpolate(pnt); return CNastVector2d( u, v ); } // interpolierten Druck an einem bel. Punkt innerhalb // des Gitters erfragen double CNastStaggeredGrid2d::pressure( const CNastPoint2d &pnt ) const { NAST_ASSERT_VALID( this ); // funktioniert vorerst nur, wenn der Punkte innhalb des Gitters // fur P ist. Fuer ausserhalb muesste erst eine Extrapolation // implementiert werden. const double p = m_gridP.interpolate(pnt); return p; } // ---------------------------------------------------------------------------- // setzen der Werte // ---------------------------------------------------------------------------- // | | | | // | | | | // | | | | // --+---------------------+-------V(i,j+1)------+---------------------+-- // | | | | // | | | | // | | | | // | | | | // | | i,j | | // | | | | // | | isFluidNorth | | // | | | | // | | +--------+ | | // --+---------------------+------| V(i,j) |-----+---------------------+-- // | | +--------+ | | // | | | | // | | iSFluidSouth | | // | | | | // | | i,j-1 | | // | | | | // | | | | // | | | | // | | | | // --+---------------------+-------V(i,j-1)------+---------------------+-- // | | | | // | | | | // | | | | // V(i,j) auf dem Rand neu berechnen void CNastStaggeredGrid2d::updateV(int i, int j) { bool isFluidNorth = isFluidNorth = cell(i,j).isFluid(); bool isFluidSouth = isFluidSouth = cell(i,j).isFluid(CNastCell2d::SOUTH ); // haeufigsten Fall zuerst // ist die Kante zwischen zwei Fluidzellen, oder zwei Hinderniszellen ? // Dann ist nichts zu tun if( isFluidNorth == isFluidSouth ) return; // dann nichts wie raus // Genau eine der beiden Zellen ist Hinderniszelle, die andere // ist eine Fluidzelle if( isFluidNorth ) // ist das Hindernis suedlich V(i,j) ? { // die Suche nach dem Randpunkt wird vom Mittelpunkt der Fluid- // zelle aus gestartet const CNastPoint2d pnt = m_gridP.pos(i,j); const CNastVector2d vecSouth = CNastVector2d(0,-1); // Randpunkt bestimmen const CNastBoundaryPoint2d pntBorder = m_geom.boundaryPoint( pnt, vecSouth ); switch( pntBorder.boundaryCondition() ) { case CNastBoundaryPoint2d::SLIP: // Rutschbedingung und Haftbedingung case CNastBoundaryPoint2d::NOSLIP: m_gridV(i,j) = 0; // bei beiden stroemt kein Fluid ueber den Rand hinaus break; case CNastBoundaryPoint2d::OUTFLOW: // Ausstroembedingung m_gridV(i,j) = m_gridV(i,j+1); // der Wert muss auf den Rand uebertragen werden break; case CNastBoundaryPoint2d::INFLOW: // Einstroembedingung m_gridV(i,j) = pntBorder.velocity().y(); // Geschwindigkeit fest break; } } if( isFluidSouth ) // ist das Hindernis noerdlich V(i,j) ? { const CNastPoint2d pnt = m_gridP.pos(i,j-1); const CNastVector2d vecNorth = CNastVector2d(0,+1); // Randpunkt bestimmen const CNastBoundaryPoint2d pntBorder = m_geom.boundaryPoint( pnt, vecNorth ); switch( pntBorder.boundaryCondition() ) { case CNastBoundaryPoint2d::OUTFLOW: // Ausstroembedingung m_gridV(i,j) = m_gridV(i,j-1); break; case CNastBoundaryPoint2d::SLIP: // Rutschbedingung case CNastBoundaryPoint2d::NOSLIP: // Haftbedingung m_gridV(i,j) = 0; break; case CNastBoundaryPoint2d::INFLOW: // Einstroembedingung m_gridV(i,j) = pntBorder.velocity().y(); break; } } } // | | | | // | | | | // | | | | // --+---------------------+---------------------+---------------------+-- // | | | | // | | | | // | isFluidWest | isFluidEast | | // | +--------+ | | // U(i-1,j) | U(i,j) | U(i+1,j) | // | +--------+ | | // | i-1,j | i,j | | // | | | | // | | | | // --+---------------------+---------------------+---------------------+-- // | | | | // | | | | // | | | | // U(i,j) auf dem Rand neu berechnen void CNastStaggeredGrid2d::updateU(int i, int j) { bool isFluidWest = cell(i,j).isFluid( CNastCell2d::WEST ); bool isFluidEast = cell(i,j).isFluid(); // ist die Kante zwischen zwei Fluidzellen, oder zwei Hinderniszellen ? if( isFluidWest == isFluidEast ) return; // dann nichts wie raus if( isFluidWest ) // ist das Hindernis oestlich von U(i,j) ? { // die Suche nach dem Randpunkt wird vom Mittelpunkt der Fluid- // zelle aus gestartet, da der auf jedem Fall innerhalb des // Gebiets liegt const CNastPoint2d pnt = m_gridP.pos(i-1,j); const CNastVector2d vecEast = CNastVector2d( +1, 0); // Randpunkt bestimmen const CNastBoundaryPoint2d pntBorder = m_geom.boundaryPoint( pnt, vecEast ); switch( pntBorder.boundaryCondition() ) { case CNastBoundaryPoint2d::SLIP: // Rutschbedingung case CNastBoundaryPoint2d::NOSLIP: // Haftbedingung m_gridU(i,j) = 0; break; case CNastBoundaryPoint2d::OUTFLOW: // Ausstroembedingung m_gridU(i,j) = m_gridU(i-1,j); break; case CNastBoundaryPoint2d::INFLOW: // Einstroembedingung m_gridU(i,j) = pntBorder.velocity().x(); break; } } if( isFluidEast ) // ist das Hindernis westlich von U(i,j) ? { const CNastPoint2d pnt = m_gridP.pos(i,j); const CNastVector2d vecWest = CNastVector2d( -1,0); // Randpunkt bestimmen const CNastBoundaryPoint2d pntBorder = m_geom.boundaryPoint( pnt, vecWest ); switch( pntBorder.boundaryCondition() ) { case CNastBoundaryPoint2d::OUTFLOW: // Ausstroembedingung m_gridU(i,j) = m_gridU(i+1,j); break; case CNastBoundaryPoint2d::SLIP: // Haftbedingung case CNastBoundaryPoint2d::NOSLIP: // Rutschbedingung m_gridU(i,j) = 0; break; case CNastBoundaryPoint2d::INFLOW: // Einstroembedingung m_gridU(i,j) = pntBorder.velocity().x(); break; } } } // | | | | // | | | | // | | | | // | | | | // --+---------------------+-------V(i,j+1)------+---------------------+-- // | | | | // | | | | // | | | | // | | | | // | | i,j | | // | | | | // | | | | // | | | | // | | | | // --+-------V(i-1,j)------+--------V(i,j)-------+-------V(i+1,j)------+-- // | | | | // | | | | // | | | | // | | | | // | | i,j-1 | | // | | | | // | | | | // | | | | // | | | | // --+---------------------+-------V(i,j-1)------+---------------------+-- // | | | | // | | | | // | | | | // | | | | // | | | | // horizontale Geschwindigkeitkomponente setzen // Achtung ! die Geschwindigkeitskomponenten duerfen nur an // den Kanten zwischen zwei Fluidzellen gesetzt werden. Die Werte auf // dem Rand ergeben sich aus den Randbedingungen und werden in den // update-Funktionen berechnet. void CNastStaggeredGrid2d::setV( int i, int j, double v ) { NAST_ASSERT_VALID( this ); NAST_ASSERT_INDEX( i > 0); NAST_ASSERT_INDEX( i < m_gridV.sizeI()-1); NAST_ASSERT_INDEX( j > 0 ); NAST_ASSERT_INDEX( j < m_gridV.sizeJ()-1); NAST_ASSERT( cell(i,j).isFluid( CNastCell2d::THIS | CNastCell2d::SOUTH )); // Die westlichen und oestlichen Werte // V(i-1,j), V(i+1,j) brauchen nicht berechnet werden. // Denn es gibt nur zwei Moeglichkeiten die eintreten koennen. // 1. V(i-1,j), V(i+1,j) liegen am Rand // dann bekommen sie ihre Werte sowieso von // einer anderen Fluidzelle // 2. Die Werte liegen im Hindernis // dann werden sie nicht direkt verwendet, sondern // immer ueber westV, eastV beim Zugriff berechnet // Es brauchen also nur die Werte V(i,j+1) und V(i,j-1) // gesetzt werden, falls es sie auf dem Rand liegen m_gridV(i,j) = v; // den neuen Wert uebernehmen updateV(i,j-1); updateV(i,j+1); } // horizontale Geschwindigkeitkomponente setzen void CNastStaggeredGrid2d::setU( int i, int j, double u ) { NAST_ASSERT_VALID( this ); NAST_ASSERT_INDEX( i > 0); NAST_ASSERT_INDEX( i < m_gridU.sizeI()-1); NAST_ASSERT_INDEX( j > 0 ); NAST_ASSERT_INDEX( j < m_gridU.sizeJ()-1); NAST_ASSERT( cell(i,j).isFluid( CNastCell2d::THIS | CNastCell2d::WEST )); m_gridU( i,j ) = u; // den neuen Wert uebernehmen updateU( i-1, j); updateU( i+1, j); } // Beim Initialiseren und wenn sich die Geometrie veraendert hat // die Randwerte und die Zellen neu berechnen void CNastStaggeredGrid2d::updateAll() { int i,j; // IDEE fuer bewegte Geometrien kann hier festgestellt werden // was sich geandert hat, damit nicht jedes mal fuer // alle U und V Elemtente update aufgerufen werden muss m_geom.setCells( m_cells ); for( i = 1; i <= sizeI()+1; i++) for( j = 1; j <= sizeJ()+1; j++) { updateU(i, j); updateV(i, j); } } // Zeitpunkt abfragen double CNastStaggeredGrid2d::time() const { NAST_ASSERT_VALID( this ); return m_time; } // Zeitpunkt setzen void CNastStaggeredGrid2d::time( double t) { NAST_ASSERT_VALID( this ); NAST_ASSERT( t > m_time ); // die Zeit kann nicht zurueck laufen :-) m_time = t; } // ---------------------------------------------------------------------------- // Debug // ---------------------------------------------------------------------------- void CNastStaggeredGrid2d::debugDump( CNastDumpContext &dumpContext ) const { CNastObject::debugDump( dumpContext ); dumpContext << "CNastStaggeredGrid2d" << "\n"; dumpContext << "\tm_box = " << m_box << "\n"; dumpContext << "\tsizeI() = " << sizeI() << "\n"; dumpContext << "\tsizeJ() = " << sizeJ() << "\n"; dumpContext << "\tm_gridU = " << m_gridU << "\n";; dumpContext << "\tm_gridV = " << m_gridV << "\n";; dumpContext << "\tm_gridP = " << m_gridP << "\n";; dumpContext << "\tm_cells = " << m_cells << "\n";; } void CNastStaggeredGrid2d::assertValid() const { // zuerst einmal AssertValid der Basisklasse aufrufen CNastObject::assertValid(); }