Navigation

    CFD中文网

    CFD中文网

    • Login
    • Search
    • 最新

    OpenFOAM有方法能够使一部分网格不参与计算吗?

    OpenFOAM
    10
    29
    3753
    Loading More Posts
    • Oldest to Newest
    • Newest to Oldest
    • Most Votes
    Reply
    • Reply as topic
    Log in to reply
    This topic has been deleted. Only users with topic management privileges can see it.
    • S
      Samuel-Tu last edited by

      想在计算过程中使得一部分网格不参与计算,有没有啥近似的方法呢?我想的是对那部分网格赋予一个很大的粘度,使他们近似于固体,不知道可不可行。

      D 1 Reply Last reply Reply Quote
      • 七
        七辆战车 last edited by

        把网格抠除?

        S 1 Reply Last reply Reply Quote
        • D
          dzw05 教授 @Samuel-Tu last edited by

          @Samuel-Tu 试过,可行。

          自主匠心,普惠仿真。

          S 林 2 Replies Last reply Reply Quote
          • S
            Samuel-Tu @七辆战车 last edited by

            @七辆战车 对,就是想达到扣除的效果

            1 Reply Last reply Reply Quote
            • S
              Samuel-Tu @dzw05 last edited by

              @dzw05 :qinqin: 感谢

              1 Reply Last reply Reply Quote
              • 我
                我是河滩 last edited by

                请问你的网格是怎么扣除的,嵌套网格有挖洞的操作,洞单元就不参与计算,你了解这个吗?

                动边界

                S 1 Reply Last reply Reply Quote
                • S
                  Samuel-Tu @我是河滩 last edited by

                  @我是河滩 不了解.......

                  1 Reply Last reply Reply Quote
                  • 林
                    林之流风 @dzw05 last edited by

                    @dzw05 请问您具体是如何处理的?:xiexie:

                    D 1 Reply Last reply Reply Quote
                    • D
                      dzw05 教授 @林之流风 last edited by

                      @林之流风 我计算的是两相流问题,把感兴趣区域中某一相的粘度设的特别大(比如1e5),该相就相当于凝固了。

                      自主匠心,普惠仿真。

                      1 Reply Last reply Reply Quote
                      • 李东岳
                        李东岳 管理员 last edited by

                        @Samuel-Tu 在 OpenFOAM有方法能够使一部分网格不参与计算吗? 中说:

                        想在计算过程中使得一部分网格不参与计算,有没有啥近似的方法呢?我想的是对那部分网格赋予一个很大的粘度,使他们近似于固体,不知道可不可行。

                        这是要模拟什么工况?凝固?

                        CFD课程 改成线上了 http://dyfluid.com/class.html
                        CFD高性能服务器 http://dyfluid.com/servers.html

                        S 1 Reply Last reply Reply Quote
                        • S
                          Samuel-Tu @李东岳 last edited by

                          @东岳 类似吧,单向流问题,想随着流场计算,我可以控制流体域内的固体形状发展

                          1 Reply Last reply Reply Quote
                          • 李东岳
                            李东岳 管理员 last edited by

                            看起来你还需要考虑边界条件的问题:shangxue:

                            CFD课程 改成线上了 http://dyfluid.com/class.html
                            CFD高性能服务器 http://dyfluid.com/servers.html

                            S 1 Reply Last reply Reply Quote
                            • S
                              Samuel-Tu @李东岳 last edited by

                              @东岳 嗯嗯,就是考虑怎么在流固交界面上把壁面剪切力修正一下

                              1 Reply Last reply Reply Quote
                              • C
                                cccrrryyy 教授 last edited by

                                借鉴一下沉浸边界的思路?

                                I don't want to survive, I want to thrive.

                                S 1 Reply Last reply Reply Quote
                                • 马乔
                                  马乔 副教授 last edited by

                                  #include "fvCFD.H"
                                  #include "pisoControl.H"
                                  
                                  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
                                  
                                  int main(int argc, char *argv[])
                                  {
                                      argList::addNote
                                      (
                                          "Transient solver for incompressible, laminar flow"
                                          " of Newtonian fluids."
                                      );
                                  
                                      #include "postProcess.H"
                                  
                                      #include "addCheckCaseOptions.H"
                                      #include "setRootCaseLists.H"
                                      #include "createTime.H"
                                      #include "createMesh.H"
                                  
                                      pisoControl piso(mesh);
                                  
                                      #include "createFields.H"
                                      #include "initContinuityErrs.H"
                                  
                                      // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
                                  
                                      Info<< "\nStarting time loop\n" << endl;
                                  
                                      DynamicList<label> lowerFaceCells(mesh.C().size());
                                      DynamicList<label>  upperFaceCells(mesh.C().size());
                                  
                                      //mesh LDU address
                                      const labelUList& lowAddr = mesh.lduAddr().lowerAddr();
                                      const labelUList& upAddr = mesh.lduAddr().upperAddr();
                                  
                                      if(mesh.cellZones().size())
                                      {
                                  	const cellZone& cellBlock = mesh.cellZones()[0];
                                  	Info << "cell block size: " << cellBlock.size() << endl;
                                  	Info << cellBlock << endl;
                                  
                                  	//low efficient
                                  	forAll(lowAddr, faceI)
                                  	{
                                  	    for(auto cellI : cellBlock)
                                  	    {
                                  		if(lowAddr[faceI] == cellI)
                                  		{
                                  		    lowerFaceCells.append(faceI);
                                  		}
                                  
                                  		if(upAddr[faceI] == cellI)
                                  		{
                                  		    upperFaceCells.append(faceI);
                                  		}
                                  	    }
                                  	}
                                      }
                                  
                                      lowerFaceCells.shrink();
                                      upperFaceCells.shrink();
                                  
                                      Info << "low Cells size: " << lowerFaceCells.size() << ", low Cells size: " << upperFaceCells.size() << endl;
                                  
                                      while (runTime.loop())
                                      {
                                          Info<< "Time = " << runTime.timeName() << nl << endl;
                                  
                                          #include "CourantNo.H"
                                  
                                          // Momentum predictor
                                  
                                          fvVectorMatrix UEqn
                                          (
                                              fvm::ddt(U)
                                            + fvm::div(phi, U)
                                            - fvm::laplacian(nu, U)
                                          );
                                  	
                                  	//modify matrix
                                  	if(mesh.cellZones().size())
                                  	{
                                  	    const cellZone& cellBlock = mesh.cellZones()[0]; //should select cellZone accor name
                                  
                                  	    scalarField& lower = UEqn.lower();
                                  	    scalarField& upper = UEqn.upper();
                                  	    scalarField& diag = UEqn.diag();
                                  
                                  
                                  	    for(auto cellI : cellBlock)
                                  	    {
                                  		diag[cellI] = 1.;
                                  		UEqn.source()[cellI] = vector::zero; //vector(0.,0.,0.); you can set your deserved value here, or read from dict, ref. fvoptions
                                  	    }
                                  
                                  	    for(auto faceI : lowerFaceCells)
                                  	    {
                                  	 	lower[faceI] = 0.;	
                                  	    }
                                  
                                  	    for(auto faceI : upperFaceCells)
                                  	    {
                                  	 	upper[faceI] = 0.;	
                                  	    }
                                  	}
                                  	else
                                  	{
                                  	    Info << "There is no cell selected" << endl;
                                  	}
                                  
                                          if (piso.momentumPredictor())
                                          {
                                              solve(UEqn == -fvc::grad(p));
                                          }
                                  
                                          // --- PISO loop
                                          while (piso.correct())
                                          {
                                              volScalarField rAU(1.0/UEqn.A());
                                              volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
                                              surfaceScalarField phiHbyA
                                              (
                                                  "phiHbyA",
                                                  fvc::flux(HbyA)
                                                + fvc::interpolate(rAU)*fvc::ddtCorr(U, phi)
                                              );
                                  
                                              adjustPhi(phiHbyA, U, p);
                                  
                                              // Update the pressure BCs to ensure flux consistency
                                              constrainPressure(p, U, phiHbyA, rAU);
                                  
                                              // Non-orthogonal pressure corrector loop
                                              while (piso.correctNonOrthogonal())
                                              {
                                                  // Pressure corrector
                                  
                                                  fvScalarMatrix pEqn
                                                  (
                                                      fvm::laplacian(rAU, p) == fvc::div(phiHbyA)
                                                  );
                                  
                                                  pEqn.setReference(pRefCell, pRefValue);
                                  
                                  		//modify p matrix
                                  		if(mesh.cellZones().size())
                                  		{
                                  		    const cellZone& cellBlock = mesh.cellZones()[0];
                                  		    //p matrix symmetry
                                  		    scalarField& offDiag = pEqn.hasUpper() ? pEqn.upper() : pEqn.lower();
                                  		    scalarField& diag = pEqn.diag();
                                  
                                  		    for(auto cellI : cellBlock)
                                  		    {
                                  			diag[cellI] = 1.;
                                  			pEqn.source()[cellI] = 0.; 
                                  		    }
                                  
                                  		    //for(auto faceI : upperFaceCells) //should find faces again?
                                  		    //{
                                  		    //    offDiag[faceI] = 0.;	
                                  		    //}
                                  		    const labelUList& offAddr = pEqn.hasUpper() ? pEqn.lduAddr().upperAddr() : pEqn.lduAddr().lowerAddr();
                                  
                                  		    forAll(offAddr, faceI)
                                  		    {
                                  			for(auto cellI : cellBlock)
                                  			{
                                  			    if(offAddr[faceI] == cellI)
                                  			    {
                                  				offDiag[faceI] = 0.;	
                                  			    }
                                  			}
                                  		    }
                                  		}
                                  
                                                  pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
                                  
                                                  if (piso.finalNonOrthogonalIter())
                                                  {
                                                      phi = phiHbyA - pEqn.flux();
                                                  }
                                              }
                                  
                                              #include "continuityErrs.H"
                                  
                                              U = HbyA - rAU*fvc::grad(p);
                                              U.correctBoundaryConditions();
                                          }
                                  
                                          runTime.write();
                                  
                                          runTime.printExecutionTime(Info);
                                      }
                                  
                                      Info<< "End\n" << endl;
                                  
                                      return 0;
                                  }
                                  
                                  toposetDict
                                  // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
                                  
                                  actions
                                  (
                                      // porousBlockage
                                      {
                                          name    porousBlockageCellSet;
                                          type    cellSet;
                                          action  new;
                                          source  boxToCell;
                                          box     (0.04 0.04 -1) (0.06 0.06 1);
                                      }
                                      {
                                          name    porousBlockage;
                                          type    cellZoneSet;
                                          action  new;
                                          source  setToCellZone;
                                          set     porousBlockageCellSet;
                                      }
                                  );
                                  

                                  UEqn.png pEqn.png
                                  :wocao:

                                  装逼没输过,吵架没赢过!

                                  S 1 Reply Last reply Reply Quote
                                  • S
                                    Samuel-Tu @cccrrryyy last edited by

                                    @cccrrryyy 对,现在我用的就是浸入边界法

                                    C 1 Reply Last reply Reply Quote
                                    • S
                                      Samuel-Tu @马乔 last edited by

                                      @马乔 这是建立了一个集合,然后对这个集合赋予不同的粘度吗

                                      1 Reply Last reply Reply Quote
                                      • 马乔
                                        马乔 副教授 last edited by

                                        modify the matrix, set the value of p and U as 0.

                                        装逼没输过,吵架没赢过!

                                        我 1 Reply Last reply Reply Quote
                                        • C
                                          cccrrryyy 教授 @Samuel-Tu last edited by

                                          @Samuel-Tu 沉浸边界法流体域和固体域都计算吧?我记得它的核心思想就是把固体域也囊括到流体域的计算里面,通过修改计算域的控制方程(类似于NS方程的一个东西?),是的计算域内某一部分的表征是固体某一部分的表征是流体。如果你是想让计算域中其中一部分不参与计算,那和普通的只计算流体域不计算固体域(固体域作为流场边界)的区别是什么呢?挺好奇的。

                                          I don't want to survive, I want to thrive.

                                          S 1 Reply Last reply Reply Quote
                                          • S
                                            Samuel-Tu @cccrrryyy last edited by

                                            @cccrrryyy 浸入边界法其中的一种方法是加入附加体积力,确实如你所说是所有计算域都会计算,只是原本的固体域由于附加体积力的影响形成旋涡,不与流体域发生质量交换了,相当于就是固体了。

                                            1 Reply Last reply Reply Quote
                                            • C
                                              cccrrryyy 教授 last edited by

                                              所以说你这个做的应该不是这个体积力的东西,是另外一种了。有什么参考文献可以分享一下么?我很好奇别的处理是什么样子的:chouchou:

                                              I don't want to survive, I want to thrive.

                                              S O 2 Replies Last reply Reply Quote
                                              • S
                                                Samuel-Tu @cccrrryyy last edited by

                                                @cccrrryyy 我现在做的就是体积力的方法,我在想有没有其他办法。。

                                                1 Reply Last reply Reply Quote
                                                • 我
                                                  我是河滩 @马乔 last edited by

                                                  @马乔 在 OpenFOAM有方法能够使一部分网格不参与计算吗? 中说:

                                                  set the value of p and U as 0.

                                                  可以给出一些有关矩阵修改的参考文献吗?目前也需要对一些给定单元赋值(通过周围正常求解网格单元的值进行插值),OpenFoam中的fixedinternalValue边界条件也有类似功能,但不是很理解。

                                                  动边界

                                                  1 Reply Last reply Reply Quote
                                                  • R
                                                    random_ran 副教授 last edited by

                                                    (Sorry, I don't have chinese input right now)

                                                    I once had the exactly same issue!

                                                    Check my thesis Chapter 6.3 :

                                                    https://scholar.uwindsor.ca/cgi/viewcontent.cgi?article=8585&context=etd

                                                    I wrote a paper in this: https://ascelibrary.org/doi/pdf/10.1061/9780784415153.ch06

                                                    Figure 2 is the original plate. (Sorry due to copy-right reason, I cannot put it here.)

                                                    In my thesis, I had to open the plate with some holes in ICEM. So, I maual dig out the block out of the computation domain.

                                                    It's really a painful process.

                                                    But it works!

                                                    This is the original plate:

                                                    69e5a95f-8963-430f-aded-3ebc52e3aed5-image.png

                                                    This is after the "dig":
                                                    97acbf4d-9894-4c21-88e7-34f655ab3fd7-image.png

                                                    Yours in CFD,

                                                    Ran

                                                    S 1 Reply Last reply Reply Quote
                                                    • 李东岳
                                                      李东岳 管理员 last edited by 李东岳

                                                              fvVectorMatrix UEqn
                                                              (   
                                                                  fvm::ddt(U)
                                                                + fvm::div(phi, U)
                                                                - fvm::laplacian(nu, U)
                                                              );
                                                      
                                                              // 固定值
                                                              labelList cells = mesh.cellZones()["null"];
                                                              vectorField value(cells.size(), Zero);
                                                              UEqn.setValues(cells, value);
                                                      
                                                              if (piso.momentumPredictor())
                                                              {
                                                                  solve(UEqn == -fvc::grad(p));
                                                              }
                                                      
                                                      

                                                      可以通过设定一个null区域,然后对区域的网格点强行赋值,但就像我之前说的,总觉得边界处需要处理。也就像 @cccrrryyy 说的,变成了浸没边界法。这个代码只是大体上是个思路,需要更深的研究。可以看看fvOptions里面一些源项的处理,都是操纵矩阵的,比如fixedTempresure,meanVelocityForce等,

                                                      CFD课程 改成线上了 http://dyfluid.com/class.html
                                                      CFD高性能服务器 http://dyfluid.com/servers.html

                                                      S 1 Reply Last reply Reply Quote
                                                      • S
                                                        Samuel-Tu @李东岳 last edited by

                                                        @东岳 嗯嗯,我试一下看看能不能行。我记得setVAlue是切断了相邻面上的联系。我处理下剪切力试试

                                                        1 Reply Last reply Reply Quote
                                                        • S
                                                          Samuel-Tu @random_ran last edited by

                                                          @random_ran 主要的,我的边界随时要改变。这样的话就会一直改网格,有点像自适应网格了。。

                                                          1 Reply Last reply Reply Quote
                                                          • O
                                                            OItoCFD @cccrrryyy last edited by

                                                            @cccrrryyy 他现在这个方法属于cut cell方法的一种吧 我也今天刚要尝试 在固体切割网格的地方要特殊处理 然后当固体是运动时候 为了保持质量守恒 要对连续性方程修正一下减去一个固体相对速度 详情见这篇论文:A fixed-grid model for simulation of a moving body in free surface flows

                                                            1 Reply Last reply Reply Quote
                                                            • 李东岳
                                                              李东岳 管理员 last edited by

                                                              最近把理论这部分弄了一下,

                                                              捕获.PNG

                                                              CFD课程 改成线上了 http://dyfluid.com/class.html
                                                              CFD高性能服务器 http://dyfluid.com/servers.html

                                                              1 Reply Last reply Reply Quote
                                                              • First post
                                                                Last post

                                                              CFD中文网 | 东岳流体 | 京ICP备15017992号-2