请教：第三代liutex涡识别方法如何在paraview后处理中实现？？
-
在paraview里面实现不容易，你可以用openfoam编程实现，比如下面的方程33
https://link.springer.com/article/10.1007/s42241-019-0048-7
-
@李东岳 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; } // ************************************************************************* //
