CFD中文网

AdjointOptimisationFoam user Mannual 链接：
https://www.openfoam.com/documentation/files/adjointOptimisationFoamManual.pdf
建议使用之前看一下连续伴随方法中的相关公式推导、伴随变量的边界条件推导以及敏感度分析，目前仅有伴随S-A模型，适用于气动外形优化。目前OpenFOAM中实现连续伴随方法的主要难点在于动量伴随方程中ATC（adjoint transpose convection）项的处理，由于该项会导致矩阵刚度较大，如果处理不当会导致伴随方程求解发散。运用伴随方程可以方便的进行各种目标函数和约束的敏感度分析，包括表面敏度和体积敏度，分别对应形状优化和拓扑优化。欧洲相关学者已经研究得非常透彻了，未来应该会在openfoam新版本慢慢放出。已经商业化的是engys公司旗下的Helyx-adjoint产品。

@王昭君 他们封装的是openfoam的cartesianMesh，即cfMesh，还是跟snappyHexMesh一样依赖参数调节。

foam-extend 里已经实现了IBM相关的求解器，可以参考。

可压缩文章里面一般叫通量差分分裂（FDS）和通量矢量分裂(FVS)。FDS中最著名当属Roe格式，FVS中比较著名的是Steger-Warming格式和Van Leer格式。个人一点见解。

目前，我主要侧重于植入我们课题组的可压缩通量格式，主要是一些基于气体动理学的一些工作。替换dbnsFlux的通量模板参数，可以切换到新的通量格式进行计算。OpenFOAM的可压缩求解器平台不够完善，总是遇到大大小小的问题。

我也是用dbnsFoam以及dbnsTurbFoam，植入了一些新的无粘通量的格式，加速收敛技术只用了LTS，湍流的case依旧算的 很慢。shenchun的那篇文章LU-SGS实现纯粹是自己另写了一个求解器，跟OpenFOAM基本没太大关系。采用fvblockMatrix确实是个突破口。 之前思考过把vector3D改为vector5D但底层代码需要改太多就放弃了。
不知楼主可否分享一下代码？考虑粘性通量和湍流的计算确实是个难点。好久没上论坛，突然看见个做可压缩隐式求解器的同道中人，幸会。

@程迪 终于找到和我研究方向相同的FOAMer了！这个代码我也研究过，虽然思路很清晰明了，与我之前Fortran编写可压缩流求解器相同。但似乎没那么FOAM了？OF里面很多东西都没有用起来，似乎只是当成了一个编程语言在用。还有一个AeroFOAM也挺久没更新了，似乎面向对象的思想更多一些？不知是否研究过？　希望以后多多交流。:cheeky:

@youmengtian 目前openfoam的可压缩求解器确实没有不可压缩求解器那么强大。我也在参考densityBased求解器做一些算法植入工作。目前有人已经实现了LU-SGS方法，不知他的实现方法和blockMarix有何联系？
附上它的文章：
Shen C, Sun F, Xia X. Implementation of density-based solver for all speeds in the framework of OpenFOAM[J]. Computer Physics Communications, 2014, 185(10): 2730-2741.
Shen C, Xia X, Wang Y, et al. Implementation of density-based implicit LU-SGS solver in the framework of OpenFOAM[J]. Advances in Engineering Software, 2016, 91: 80-88.
期待进一步交流！

各位：
　　最近在研究这个求解器：
　　compressibleFoam，发现里面的求解思路和我之前采用fortran语言是一样的。思路就是按边来循环，分别对internalFace和boundaryFace进行处理，求得单元交界面上的通量。
　　控制方程如下
　　
其中一些量如下
　　
　　
　　
　　下面这个是求内部面无粘通量的函数

/// Loop over internal faces and get face flux from L and R cell center
label leftCell , rightCell;
forAll( mesh.owner() , iface ) {
/// Get the left and right cell index
leftCell = mesh.owner()[iface];
rightCell = mesh.neighbour()[iface]; 
/// Approximate Riemann solver at interface 
scalar lambda = (*fluxSolver)( &rho[leftCell]  , &U[leftCell]  , &p[leftCell] , /// Left
     &rho[rightCell] , &U[rightCell] , &p[rightCell] ,          /// Right
     &massFlux[iface] , &momFlux[iface] , &energyFlux[iface] ,
     &nf[iface] );
/// Multiply with face area to get face flux
massFlux[iface] *= mesh.magSf()[iface];
momFlux[iface] *= mesh.magSf()[iface];
energyFlux[iface] *= mesh.magSf()[iface];
localDt[ leftCell ]  += lambda * mesh.magSf()[iface];
localDt[ rightCell ] += lambda * mesh.magSf()[iface];
}

下面这个是求物面边界条件的边界面上的无粘通量的函数

\*---------------------------------------------------------------------------*/
   /////////////////////////////////////////////////////////////////////////////////
   /// Remember BC's for Gas Dynamics solvers
   /// are based on characteristics and cannot be applied
   /// to individual variables. Rather it should be applied
   /// to the state vector as a whole. OpenFOAM framework
   /// does not allow us to do this without breaking the
   /// convention.
   /////////////////////////////////////////////////////////////////////////////////  
   /// Loop over boundary patches and determine the
   /// type of boundary condition of the patch
   /////////////////////////////////////////////////////////////////////////////////
   forAll ( mesh.boundaryMesh() , ipatch ) {
     scalar bflux[5];
     word BCTypePhysical = mesh.boundaryMesh().physicalTypes()[ipatch];
     word BCType         = mesh.boundaryMesh().types()[ipatch];
     word BCName         = mesh.boundaryMesh().names()[ipatch];
     const UList<label> &bfaceCells = mesh.boundaryMesh()[ipatch].faceCells();
   /////////////////////////////////////////////////////////////////////////////////
                                /// Slip Wall BC ///
   /////////////////////////////////////////////////////////////////////////////////
     if( BCTypePhysical == "slip" || BCTypePhysical == "symmetry" ) {
       forAll( bfaceCells  , iface ) {
         /// Extrapolate wall pressure
         scalar p_e = p[ bfaceCells[iface] ];
         vector normal = nf.boundaryField()[ipatch][iface];
         scalar face_area = mesh.magSf().boundaryField()[ipatch][iface];
         scalar lambda = std::fabs( U[ bfaceCells[iface] ] & normal ) +
                         std::sqrt ( gama * p_e / rho[ bfaceCells[iface] ] );
         momResidue[ bfaceCells[iface] ] += p_e * normal * face_area;
         localDt[ bfaceCells[iface] ]    += lambda * face_area;
       }
     }

之后我们将边界上的通量从owner单元减去，加到nighbour单元去

   /// Clear all residues
   massResidue = dimensionedScalar( "", massResidue.dimensions() , 0.0 );
   momResidue = dimensionedVector( "" , momResidue.dimensions() , vector( 0.0 , 0.0 , 0.0 ) );
   energyResidue = dimensionedScalar( "" , energyResidue.dimensions() ,  0.0 );
   /// Loop over each face and add the flux to left and subtract from the right
   forAll( mesh.owner() , iface ) {
     /// Store left and right cell reference
     leftCell = mesh.owner()[iface];
     rightCell = mesh.neighbour()[iface];
     /// Note that the normal vector to a face will point
     /// from the owner cell to the neighbour cell (L to R)
     /// Add to left cell
     massResidue[leftCell]    += massFlux[iface];
     momResidue[leftCell]     += momFlux[iface];
     energyResidue[leftCell]  += energyFlux[iface];
     /// Subtract from right cell
     massResidue[rightCell]   -= massFlux[iface];
     momResidue[rightCell]    -= momFlux[iface];
     energyResidue[rightCell] -= energyFlux[iface];
   }

原来编程就是一直按照这种思路来的，可是转到OpenFOAM 下一时适应不过来。看到官方的求解器，根本摸不清楚边界条件和有限体积是如何交互的。有人能解释一下这和我上述思路的差别和联系吗？
大家编写自己的求解器的时候都是怎么实现的？希望看透这一切联系的牛人能给予一点提示和点拨。跪谢！

暂时还没有，对OpenFOAM的边界条件处理还不是很懂，请问这个是OpenFOAM自带的吗？

SOWFA 是美国能源实验室基于OpenFOAM开发的软件，主要用于计算风力机的有关问题https://nwtc.nrel.gov/SOWFA

我的问题如下：已经在能量控制方程当中加入了地表因太阳辐射对大气热量影响的源项，其边界条件需模拟一天当中太阳的照射情况，目前打算采用正弦曲线，请问如何实现边界条件或源项场按照已知规律随时间变化呢？恳请高手不吝赐教。