CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Programming & Development

Compressible Dynamic Mesh Solver in Foam-Extend 3.1

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree1Likes
  • 1 Post By matthewshorter

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   June 21, 2016, 10:09
Default Compressible Dynamic Mesh Solver in Foam-Extend 3.1
  #1
New Member
 
Matthew Shorter
Join Date: Feb 2016
Posts: 10
Rep Power: 10
matthewshorter is on a distinguished road
Hello,
I am looking to implement a moving mesh solver in foam-extend 3.1. The problem we are studying is a shock tunnel with an actuated valve connecting each section. We have successfully created a moving mesh based on the movingConeTopo case under icoDyMFoam, but we are having considerable difficulty implementing mesh motion into a compressible flow solver. Our goal is to effectively recreate "rhoCentralDyMFoam" for foam-extend. Using the SOD shock tube case, we have confirmed that the custom solver returns accurate results for a stationary mesh, but it crashes with a floating point error due to negative pressure occurring at the trailing side of the wall when a new mesh layer is added.
Here is the code for the relative flux we are currently using and is placed after the directed interpolation of primitive fields after faces:
Code:
if (mesh.moving()){

Info<<"Mesh flux is being accounted for\n";  
 phiv_pos -= mesh.phi();// fvc::meshPhi(rho,U);
        

 phiv_neg -= mesh.phi();//fvc::meshPhi(rho,U);



        Info <<endl <<phiv_neg <<endl<<"printed phiv_neg";
         }
and later here is the pressure made absolute
Code:
  phiEp +=mesh.phi()*(a_pos*p_pos + a_neg*p_neg);
I've attached the mesh at the point of crash. This is right after the dynamic mesh added a new layer. Note the negative temperature at the wall, which obviously blows up the solver....

Below is the solver in its entirety:
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
    rhoCentralDyMFoam

Description
    Density-based compressible flow solver based on central-upwind schemes of
    Kurganov and Tadmor with support for mesh-motion and topology changes

\*---------------------------------------------------------------------------*/

#include "fvCFD.H"

#include "dynamicFvMesh.H"

#include "basicPsiThermo.H"

#include "zeroGradientFvPatchFields.H"
#include "fixedRhoFvPatchScalarField.H"

#include "motionSolver.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
    #include "setRootCase.H"
    #include "createTime.H"
    #include "createDynamicFvMesh.H"

    #include "createFields.H"
    #include "readTimeControls.H"
#include "readThermophysicalProperties.H"   //?xz
    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    #include "readFluxScheme.H"

    dimensionedScalar v_zero("v_zero", dimVolume/dimTime, 0.0);
    scalar CoNum = 0.0;
    scalar meanCoNum = 0.0;
    Info<< "\nStarting time loop\n" << endl;
    while (runTime.run())
      {
    runTime++;
        Info<< "Time = " << runTime.timeName() << nl << endl;
    bool meshChanged =        mesh.update();
    reduce(meshChanged, orOp<bool>());
        // --- Directed interpolation of primitive fields onto faces
        surfaceScalarField rho_pos=fvc::interpolate(rho, pos,"reconstruct(rho)");
        surfaceScalarField rho_neg=fvc::interpolate(rho, neg,"reconstruct(rho)");

        surfaceVectorField rhoU_pos=fvc::interpolate(rhoU, pos,"reconstruct(U)");
        surfaceVectorField rhoU_neg=fvc::interpolate(rhoU, neg, "reconstruct(U)");

        volScalarField rPsi= 1.0/psi;
        surfaceScalarField rPsi_pos=fvc::interpolate(rPsi, pos,"reconstruct(T)");
        surfaceScalarField rPsi_neg=fvc::interpolate(rPsi, neg,"reconstruct(T)");

        surfaceScalarField e_pos=fvc::interpolate(e, pos,"reconstruct(T)");
        surfaceScalarField e_neg=fvc::interpolate(e, neg,"reconstruct(T)");
        surfaceVectorField U_pos= rhoU_pos/rho_pos;
        surfaceVectorField U_neg = rhoU_neg/rho_neg;
        surfaceScalarField p_pos= rho_pos*rPsi_pos;
        surfaceScalarField p_neg = rho_neg*rPsi_neg;
        surfaceScalarField phiv_pos= U_pos & mesh.Sf();
        surfaceScalarField phiv_neg= U_neg & mesh.Sf();






    if (mesh.moving()){

Info<<"Mesh flux is being accounted for\n";  
 phiv_pos -= mesh.phi();// fvc::meshPhi(rho,U);
        

 phiv_neg -= mesh.phi();//fvc::meshPhi(rho,U);



        Info <<endl <<phiv_neg <<endl<<"printed phiv_neg";
         }

    Info<<"  adjusted for moving mesh ";
        volScalarField c= sqrt(thermo.Cp()/thermo.Cv()*rPsi);

    surfaceScalarField cSf_pos =         fvc::interpolate(c, pos,"reconstruct(T)")*mesh.magSf() ;
        surfaceScalarField cSf_neg=    fvc::interpolate(c, neg,"reconstruct(T)")*mesh.magSf();

        surfaceScalarField ap =   max(max(phiv_pos + cSf_pos, phiv_neg + cSf_neg), v_zero);
        surfaceScalarField am =  min(min(phiv_pos - cSf_pos, phiv_neg - cSf_neg), v_zero);

        surfaceScalarField a_pos= ap/(ap - am);

        surfaceScalarField amaxSf("amaxSf", max(mag(am), mag(ap)));
    
        surfaceScalarField aSf= am*a_pos;

        if (fluxScheme == "Tadmor")
        {
            aSf = -0.5*amaxSf;
            a_pos = 0.5;
        }

        surfaceScalarField a_neg= 1.0 - a_pos;

        phiv_pos *= a_pos;
        phiv_neg *= a_neg;

        surfaceScalarField aphiv_pos= phiv_pos - aSf;
        surfaceScalarField aphiv_neg= phiv_neg + aSf;



        amaxSf = max(mag(aphiv_pos), mag(aphiv_neg));
    
        #include "compressibleCourantNo.H"
    #include "readTimeControls.H"
    #include "setDeltaT.H"

       surfaceScalarField phi("phi",aphiv_pos*rho_pos + aphiv_neg*rho_neg);


        surfaceVectorField phiUp=
            (aphiv_pos*rhoU_pos + aphiv_neg*rhoU_neg)
          + (a_pos*p_pos + a_neg*p_neg)*mesh.Sf();

        surfaceScalarField phiEp=
            aphiv_pos*(rho_pos*(e_pos + 0.5*magSqr(U_pos)) + p_pos)
          + aphiv_neg*(rho_neg*(e_neg + 0.5*magSqr(U_neg)) + p_neg)
          + aSf*p_pos - aSf*p_neg;




      phiEp +=mesh.phi()*(a_pos*p_pos + a_neg*p_neg);

 Info << "\n print out check phiEp\n";

    Info << "\n beginning to solve the density and momentum equations\n";

//        volScalarField muEff("muEff", turbulence->muEff());
        volTensorField tauMC("tauMC", mu*dev2(fvc::grad(U)().T()));



    // solve density

        solve(fvm::ddt(rho) + fvc::div(phi));
    //Info << "Solving Momentum\n";
        // --- Solve momentum
        solve(fvm::ddt(rhoU) + fvc::div(phiUp));
    
        U.dimensionedInternalField() =
            rhoU.dimensionedInternalField()
           /rho.dimensionedInternalField();
        U.correctBoundaryConditions();
        rhoU.boundaryField() = rho.boundaryField()*U.boundaryField();
    volScalarField rhoBydt(rho/runTime.deltaT());
        if (!inviscid)
        {
            solve
            (
                fvm::ddt(rho, U) - fvc::ddt(rho, U)
              - fvm::laplacian(mu, U)
              - fvc::div(tauMC)
            );
      Info << "Flow is not inviscid";
          rhoU = rho*U;
        }
    //Info <<"Define sigmaDotU\n";
        // --- Solve energy
        surfaceScalarField sigmaDotU=(
                fvc::interpolate(mu)*mesh.magSf()*fvc::snGrad(U)
              + (mesh.Sf() & fvc::interpolate(tauMC))
            )
            & (a_pos*U_pos + a_neg*U_neg);
    Info<<"Solving Energy Eqn\n";
        solve
        (
            fvm::ddt(rhoE)
          + fvc::div(phiEp)
      -fvc::div(sigmaDotU)         
        );

        e = rhoE/rho - 0.5*magSqr(U);
        e.correctBoundaryConditions();
        thermo.correct();
    //Info <<"\n corrected thermo and using the solve functions \n";
        rhoE.boundaryField() =
            rho.boundaryField()*
            (
                e.boundaryField() + 0.5*magSqr(U.boundaryField())
            );

        if (!inviscid)
        {
        //turned off turbulence for now MS 6/7/16
            solve
            (
                fvm::ddt(rho, e) - fvc::ddt(rho, e)
        - fvm::laplacian(thermo.alpha(), e)
            );
            thermo.correct();
            rhoE = rho*(e + 0.5*magSqr(U));
        }

        p.dimensionedInternalField() =
            rho.dimensionedInternalField()
           /psi.dimensionedInternalField();
        p.correctBoundaryConditions();
        rho.boundaryField() = psi.boundaryField()*p.boundaryField();

    // turbulence->correct(); //see above; turned off 6/7/16 MS

        runTime.write();

        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
    }

    Info<< "End\n" << endl;

    return 0;
}

// ************************************************************************* //
Attached Images
File Type: jpg tempError.jpg (93.2 KB, 46 views)
OlegSutyrin likes this.
matthewshorter is offline   Reply With Quote

Old   September 19, 2023, 23:14
Default
  #2
Member
 
Michael Sukham
Join Date: Mar 2020
Location: India
Posts: 79
Rep Power: 6
2538sukham is on a distinguished road
Have you solved this particular issue? When layer addition or removal occurs, the simulation crashed. I am using foamExtend 3.2 and the LAR is based on linearValveLayersFvMesh. Is it to do with mapping when topoChange occurs?
With regards
2538sukham is offline   Reply With Quote

Reply

Tags
compressible, dynamic mesh, foam extend 3.1

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
decomposePar problem: Cell 0contains face labels out of range vaina74 OpenFOAM Pre-Processing 37 July 20, 2020 05:38
[mesh manipulation] Importing Multiple Meshes thomasnwalshiii OpenFOAM Meshing & Mesh Conversion 18 December 19, 2015 18:57
[blockMesh] BlockMesh FOAM warning gaottino OpenFOAM Meshing & Mesh Conversion 7 July 19, 2010 14:11
[Gmsh] Import problem ARC OpenFOAM Meshing & Mesh Conversion 0 February 27, 2010 10:56


All times are GMT -4. The time now is 18:28.