FlowModelSlipJump.h

Go to the documentation of this file.

// This file os part of FVM
// Copyright (c) 2012 FVM Authors
// See LICENSE file for terms.

#ifndef _FLOWMODELSLIPJUMP_H_
#define _FLOWMODELSLIPJUMP_H_

// this file is meant to be included inside the FlowModel::Impl class
// the code is here just because FlowModel_impl.h has grown too big

void slipJumpMomentumBC(const StorageSite& faces,
                        const Mesh& mesh,
                        GenericBCS<VectorT3,DiagTensorT3,T>& gbc,
                        const T accomodationCoefficient,
                        const FloatValEvaluator<VectorT3>& bVelocity)
{
  const StorageSite& cells = mesh.getCells();

  const CRConnectivity& faceCells = mesh.getFaceCells(faces);

  const VectorT3Array& faceArea =
    dynamic_cast<const VectorT3Array&>(_geomFields.area[faces]);

  const VectorT3Array& faceCentroid =
    dynamic_cast<const VectorT3Array&>(_geomFields.coordinate[faces]);
    
  const VectorT3Array& cellCentroid =
    dynamic_cast<const VectorT3Array&>(_geomFields.coordinate[cells]);
    
  const TArray& faceAreaMag =
    dynamic_cast<const TArray&>(_geomFields.areaMag[faces]);

  const VectorT3Array& V = dynamic_cast<const VectorT3Array&>(_flowFields.velocity[cells]);

  const TArray& p = dynamic_cast<const TArray&>(_flowFields.pressure[cells]);
  const TArray& mu = dynamic_cast<const TArray&>(_flowFields.viscosity[cells]);

  const T opPressure( _options["operatingPressure"]);
  const T opTemperature(_options["operatingTemperature"]);
  const T molWt(_options["molecularWeight"]);
  const T Rgas = 8314.472/molWt;

  const bool incompressible = _options.incompressible;
  
  const int nFaces = faces.getCount();

  for(int f=0; f<nFaces; f++)
  {
      const int c0 = faceCells(f,0);

      // face normal

      const VectorT3 en(faceArea[f]/faceAreaMag[f]);
      
      // normal velocity
      const T Vn = dot(V[c0],en);

      // wall parallel component of cell velocity is Vc - Vn * en
      const VectorT3 Vp(V[c0] - Vn*en);

      const VectorT3 ds(faceCentroid[f]-cellCentroid[c0]);

      // normal distance between face and cell centroid
      const T dn = dot(ds,en);

      
      const T pAbs = incompressible ? opPressure : (p[c0]+opPressure);
      
      
      const T muCell = mu[c0];
      const T MeanFreePath = muCell/pAbs*sqrt(0.5*M_PI*Rgas*opTemperature);

      const T coeff = accomodationCoefficient*MeanFreePath/ (dn+(accomodationCoefficient*MeanFreePath));
      const VectorT3 Vwp(Vp*coeff);

      // specified velocity at the boundary
      const VectorT3 bv = bVelocity[f];
      
      // normal component of the specified boundary velocity
      const VectorT3 Vwn(bv-dot(bv,en)*bv);

      // the velocity at the boundary is the parallel component
      // computed above + the normal from the bc specfication
      const VectorT3 Vw(Vwn + Vwp);
      gbc.applyDirichletBC(f,Vw);
     
  }
}

#endif