CFD中文网

    CFD中文网

    • 登录
    • 搜索
    • 最新

    请教:第三代liutex涡识别方法如何在paraview后处理中实现??

    OpenFOAM
    2
    3
    1146
    正在加载更多帖子
    • 从旧到新
    • 从新到旧
    • 最多赞同
    回复
    • 在新帖中回复
    登录后回复
    此主题已被删除。只有拥有主题管理权限的用户可以查看。
    • 波
      波流力 最后由 编辑

      看到相关的文献都说第三代liutex涡识别方法具有极大的优势,但是如何在paraview后处理中去实现呢??有没有尝试过的大佬分享一下如何实现?

      1 条回复 最后回复 回复 引用
      • 李东岳
        李东岳 管理员 最后由 李东岳 编辑

        在paraview里面实现不容易,你可以用openfoam编程实现,比如下面的方程33

        https://link.springer.com/article/10.1007/s42241-019-0048-7

        dab85dc8-dc56-4597-91b8-643c7e70e9c9-image.png

        线上CFD课程开始报名:http://www.dyfluid.com/class.html

        CFD高性能服务器 http://dyfluid.com/servers.html

        李东岳 1 条回复 最后回复 回复 引用
        • 李东岳
          李东岳 管理员 @李东岳 最后由 编辑

          @李东岳 http://www.jhydrodynamics.com/en/download-of-liutex-code/ 这个已经有了个开源的代码

          /*---------------------------------------------------------------------------*\
            =========                 |
            \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
             \\    /   O peration     |
              \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
               \\/     M anipulation  |
          -------------------------------------------------------------------------------
          License
              This file is part of OpenFOAM.
          
              OpenFOAM is free software: you can redistribute it and/or modify it
              under the terms of the GNU General Public License as published by
              the Free Software Foundation, either version 3 of the License, or
              (at your option) any later version.
          
              OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
              ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
              FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
              for more details.
          
              You should have received a copy of the GNU General Public License
              along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
          
          Application
              stressComponents
          
          Description
              Calculates and writes the scalar fields of the six components of the stress
              tensor sigma for each time.
          
          \*---------------------------------------------------------------------------*/
          
          #include "fvCFD.H"
          #include "OFstream.H"
          
          // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
          tensor rotation(vector u, vector v)
          {
              scalar eps = 1.0e-10;
              tensor q(tensor::zero); 
              vector a(
                         u.y()*v.z()-u.z()*v.y(),
                         u.z()*v.x()-u.x()*v.z(),
                         u.x()*v.y()-u.y()*v.x()
                       );
          
             if(Foam::mag(a)<eps)
             {
                   q.xx()= 1.0;
                   q.yy()= 1.0;
                   q.zz()= 1.0;
             }
             else
             {
                   a=a/Foam::mag(a);
                   scalar t=Foam::sign( u & v);
                   scalar alpha = Foam::acos(t);
                   scalar c = Foam::cos(alpha);
                   scalar s = Foam::sin(alpha);
          
                   q.xx()=Foam::pow(a.x(),2)*(1-c)+c;
                   q.xy()=a.x()*a.y()*(1-c)-a.z()*s;
                   q.xz()=a.x()*a.z()*(1-c)+a.y()*s;
           
                   q.yx() = a.y()*a.x()*(1-c)+a.z()*s;
                   q.yy() = Foam::pow(a.y(),2)*(1-c)+c;
                   q.yz() = a.y()*a.z()*(1-c)-a.x()*s;
          
                   q.zx() = a.z()*a.x()*(1-c)-a.y()*s;
                   q.zy() = a.z()*a.y()*(1-c)+a.x()*s;
                   q.zz() = Foam::pow(a.z(),2)*(1-c)+c;
             }
             return q;
          }
          
          
          int main(int argc, char *argv[])
          {
              timeSelector::addOptions();
          
              #include "setRootCase.H"
              #include "createTime.H"
          
              instantList timeDirs = timeSelector::select0(runTime, args);
          
              #include "createMesh.H"
          
              forAll(timeDirs, timeI)
              {
                  runTime.setTime(timeDirs[timeI], timeI);
          
                  Info<< "Time = " << runTime.timeName() << endl;
          
                  IOobject Uheader
                  (
                      "U",
                      runTime.timeName(),
                      mesh,
                      IOobject::MUST_READ
                  );
                  
                  vector z0(0,0,1);
          
                  // Check U exists
                  if (Uheader.headerOk())
                  {
                      mesh.readUpdate();
          
                      Info<< "    Reading U" << endl;
                      volVectorField U(Uheader, mesh);
                      volVectorField Rotex
                      (
                          IOobject
                          (
                              "Rotex",
                              runTime.timeName(),
                              mesh,
                              IOobject::NO_READ,
                              IOobject::AUTO_WRITE
                          ),
                          mesh,
                          dimensionedVector("zero", dimless, vector::zero)
                      );
                      volTensorField gradU = fvc::grad(U);
                      
                      forAll(U,cellI)
                      {
                         scalar aa=-(  gradU.component(tensor::XX)()[cellI]
                                     + gradU.component(tensor::YY)()[cellI]
                                     + gradU.component(tensor::ZZ)()[cellI]
                                    );
                         tensor a=gradU[cellI];
                         tensor tt =a & a;
                         scalar bb = -0.5*(tt.xx()+tt.yy()+tt.zz()-pow((tt.xx()+tt.yy()+tt.zz()),2));
                         scalar cc = -( a.xx()*(a.yy()*a.zz()-a.yz()*a.zy())
                                       -a.xy()*(a.yx()*a.zz()-a.yz()*a.zx())
                                       +a.xz()*(a.yx()*a.zy()-a.yy()*a.zx())
                                      );
                         scalar delta = 18.0*aa*bb*cc-4.0*Foam::pow(aa,3)*cc+Foam::pow(aa,2)*pow(bb,2)-4.0*Foam::pow(bb,3)-27.0*Foam::pow(cc,2);
                         scalar qq=(Foam::pow(aa,2),3*bb)/9.0;
                         scalar rr=(2.0*Foam::pow(aa,3)-9.0*aa*bb+27.0*cc)/54.0;
                         delta=-delta/108.0;
                        
                         if(delta > 0.0)
                         {
                           
                            scalar aaaa=-sign(rr)*Foam::pow(abs(rr)+Foam::sqrt(delta),scalar(1.0/3.0));
                            scalar bbbb=0;
           
                            if( aaaa !=0.0)
                               bbbb=qq/aaaa;
           
          		  scalar eig1c_r=-0.5*(aaaa+bbbb)-aa/3.0;
          		  scalar eig1c_i=0.5*Foam::sqrt(scalar(3.0))*(aaaa-bbbb);
          		  scalar eig2c_r=-0.5*(aaaa+bbbb)-aa/3.0;
          		  scalar eig2c_i=-0.5*Foam::sqrt(scalar(3.0))*(aaaa-bbbb);
                            scalar eig3r = aaaa+bbbb-aa/3.0;
                            scalar delta1(0.0),delta2(0.0),delta3(0.0);
                            delta1 = (a.xx()-eig3r)*(a.yy()-eig3r) -a.yx()*a.xy();
                            delta2 = (a.yy()-eig3r)*(a.zz()-eig3r) -a.yz()*a.zy();
                            delta3 = (a.xx()-eig3r)*(a.zz()-eig3r) -a.zx()*a.xz();
                           
                            if(delta1==0.0 && delta2==0.0 && delta3==0.0)
          		  {
                                FatalErrorInFunction
          	                    << "delta1-3 are:" << delta1 << ","
                  	            << delta2 <<","<< delta3
                          	    << exit(FatalError);
                            }
                            vector vr(0,0,0);
                            
                            if(    Foam::mag(delta1)>=Foam::mag(delta2) 
                                &&  Foam::mag(delta1)>=Foam::mag(delta3)
                               ) 
                            {
                                  vr.x()=(-(a.yy()-eig3r)*a.xz()+a.xy()*a.yz())/delta1;
                                  vr.y()= (a.yx()*a.xz() - (a.xx()-eig3r)*a.yz())/delta1;
                                  vr.z()=1.0;
                            }
                            else if(    Foam::mag(delta2)>=Foam::mag(delta1) 
                                   &&  Foam::mag(delta2)>=Foam::mag(delta3)
                                 )
                            {
                                  vr.x()=1.0;
                                  vr.y()=(-(a.zz()-eig3r)*a.yx()+a.yz()*a.zx())/delta2;
                                  vr.z()=(a.zy()*a.yx() - (a.yy()-eig3r)*a.zx())/delta2;
                            }   
                            else if(    Foam::mag(delta3)>=Foam::mag(delta1) 
                                   &&  Foam::mag(delta3)>=Foam::mag(delta2)
                                 )
                            {
                                  vr.x()=(-(a.zz()-eig3r)*a.xy()+a.xz()*a.zy())/delta3;
                                  vr.y()= 1.0 ;
                                  vr.z()=(a.zx()*a.xy() - (a.xx()-eig3r)*a.zy())/delta3;
                            }   
                            else
                                 FatalErrorInFunction<< "vr error"<< exit(FatalError);
                          
                            vr = vr/Foam::sqrt(vr & vr);
           
                            tensor qqq=rotation(z0,vr);
                            tensor vg(qqq.T() & a);
                         	  vg=vg  & qqq; 
                                      
                            scalar alpha = 0.5*Foam::sqrt(Foam::pow(vg.yy()-vg.xx(),2)+Foam::pow(vg.yx()+vg.xy(),2)); 
                            scalar beta = 0.5*(vg.yx()-vg.xy());
                          
                            if(Foam::magSqr(beta) > Foam::magSqr(alpha))
                            {  
                                 scalar rm=0.0;
                                 if(beta > 0.0)
                                    rm=2.0*(beta-alpha);
                                 else
                                    rm=2.0*(beta+alpha);
                                    
                                 Rotex[cellI]=rm*vr; 
                            }
                         }
                     }
                      Rotex.write();
                      Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
                          << "  ClockTime = " << runTime.elapsedClockTime() << " s"
                          << nl << endl;
          
                  }
                  else
                  {
                      Info<< "    No U" << endl;
                  }
          
              }
          
              Info<< "End" << endl;
          
              return 0;
          }
          
          
          // ************************************************************************* //
          
          

          线上CFD课程开始报名:http://www.dyfluid.com/class.html

          CFD高性能服务器 http://dyfluid.com/servers.html

          1 条回复 最后回复 回复 引用
          • First post
            Last post