[KarthickRajkumar]

Dear all,

Greetings!

I am trying to do FSI simulation for a musical application in Foam-extend 3.2. Since the application involves compressible flow solver, I had to implement a compressible solver in FSI and run it first with HronTurekfsi tutorial case (tutorial case under fsiFoam solver) and then use it for my application.

1. rhoPisoFlow: So initially I learnt how equations are coded, did some prework and implemented rhoPisoFoam solver in FSI which I call it rhoPisoFlow. When I run it with HronTurekcase, the simulation stops at 2 s due to High courant number (time step set to be 1e-3).

2. rhoPimpleFlow: I learnt then that Pimple solver is good for unsteady compressible flow simulations because of relaxation factors and nCorrectors. So I implemented rhoPimpleFoam, which I call it rhoPimpleFlow. This one when I run with HronTurek case, the simulation stops very early around 0.4 s, again due to high courant number (time step 1e-3)

In the original HronTurek case which is an incompressible FSI case, the time step set to be 1e-3. So I kept the same. Shouldn't the same time step work for compressible case too? or will it fix this problem by reducing it further? If there should be any problem with implementation, I could explain how I implemented too.

I am facing trouble with finding where lies the problem actually? If atleast I know where lies the problem, life would not have been this difficlult! So I thought of creating a thread to get suggestion from people who can help me out on this!

Should you require more information, I can share it too!

Thanks,
Karthick

[KarthickRajkumar]

Hai,

I tried to fix the error and some how managed to get the rhoPimpleflow code work, but this time the simulation doesn't converge at all. I had set the number of fsi loop to be 45, played with relaxation factors, but that couldn't fix the convergence problem either!

I will share the code of the implemented rhoPimpleflow model here, if anyone has already implemented this or have some idea about it, please let me know if I am making any obvious mistakes!

So, the evolve function in rhoPimpleFlow.C looks this way,

void rhoPimpleFlow::evolve()
{
    Info << "Evolving flow model" << endl;

    const fvMesh& mesh = flowModel::mesh();

	// ReadPISOControls.H begins

	    int nOuterCorr(readInt(flowProperties().lookup("nOuterCorrectors")));
	    int nCorr(readInt(flowProperties().lookup("nCorrectors")));

	    int nNonOrthCorr =
		flowProperties().lookupOrDefault<int>("nNonOrthogonalCorrectors", 0);

	    bool momentumPredictor =
		flowProperties().lookupOrDefault<Switch>("momentumPredictor", true);

	    bool transonic =
		flowProperties().lookupOrDefault<Switch>("transonic", false);

	// ReadPISOControls.H ends

		if (nOuterCorr != 1)
		{
		    p_.storePrevIter();
		    rho_.storePrevIter();
		}

	   if(mesh.moving())
	    {
		// Make the fluxes relative
		phi_ -= fvc::interpolate(rho_)*fvc::meshPhi(rho_, U_);
	    }

		// CompressibleCourantNo.H begins
		scalar CoNum = 0.0;
		scalar meanCoNum = 0.0;
		scalar velMag = 0.0;

		if (mesh.nInternalFaces())
		{
		    surfaceScalarField phiOverRho = mag(phi_)/fvc::interpolate(rho_);

		    surfaceScalarField SfUfbyDelta =
			mesh.surfaceInterpolation::deltaCoeffs()*phiOverRho;

		    CoNum = max(SfUfbyDelta/mesh.magSf()).value()*runTime().deltaT().value();

		    meanCoNum = (sum(SfUfbyDelta)/sum(mesh.magSf())).value()*
			runTime().deltaT().value();

		    velMag = max(phiOverRho/mesh.magSf()).value();
		}

		Info<< "Courant Number mean: " << meanCoNum
		    << " max: " << CoNum
		    << " velocity magnitude: " << velMag
		    << endl;
		// CompressibleCourantNo.H ends


		// rhoEqn.H begins
		{
		    solve(fvm::ddt(rho_) + fvc::div(phi_));
		}
		// rhoEqn.H ends

        // --- Pressure-velocity PIMPLE corrector loop
        for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
        {
            // "UEqn.H"

		// Solve the Momentum equation

		tmp<fvVectorMatrix> UEqn
		(
		    fvm::ddt(rho_, U_)
		  + fvm::div(phi_, U_)
		  + turbulence_->divDevRhoReff(U_)
		);

		if (oCorr == nOuterCorr - 1)
		{
		    if (mesh.solutionDict().relax("UFinal"))
		    {
			UEqn().relax(mesh.solutionDict().relaxationFactor("UFinal"));
		    }
		    else
		    {
			UEqn().relax(1);
		    }
		}
		else
		{
		    UEqn().relax();
		}

		volScalarField rUA = 1.0/UEqn().A();

		if (momentumPredictor)
		{
		    if (oCorr == nOuterCorr-1)
		    {
			solve(UEqn() == -fvc::grad(p_), mesh.solutionDict().solver("UFinal"));
		    }
		    else
		    {
			solve(UEqn() == -fvc::grad(p_));
		    }
		}
		else
		{
		    U_ = rUA*(UEqn().H() - fvc::grad(p_));
		    U_.correctBoundaryConditions();
		}
		
		//UEqn ends            
		
		// hEqn.H begins
		{
		    fvScalarMatrix hEqn
		    (
			fvm::ddt(rho_, h_)
		      + fvm::div(phi_, h_)
		      - fvm::laplacian(turbulence_->alphaEff(), h_)
		     ==
			DpDt_
		    );

		    if (oCorr == nOuterCorr-1)
		    {
			hEqn.relax();
			hEqn.solve(mesh.solutionDict().solver("hFinal"));
		    }
		    else
		    {
			hEqn.relax();
			hEqn.solve();
		    }

		    thermo_.correct();
		}

		// hEqn.H ends

            // --- PISO loop
            for (int corr = 0; corr < nCorr; corr++)
            {
				//PEqn Begins

		rho_ = thermo_.rho();

		volScalarField rUA = 1.0/UEqn().A();
		U_ = rUA*UEqn().H();

		if (nCorr <= 1)
		{
		    UEqn.clear();
		}

		if (transonic)
		{
		    surfaceScalarField phid
		    (
			"phid",
			fvc::interpolate(psi_)
		       *(
			    (fvc::interpolate(U_) & mesh.Sf())
			  + fvc::ddtPhiCorr(rUA, rho_, U_, phi_)
			)
		    );

		    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
		    {
			fvScalarMatrix pEqn
			(
			    fvm::ddt(psi_, p_)
			  + fvm::div(phid, p_)
			  - fvm::laplacian(rho_*rUA, p_)
			);

			if
			(
			    oCorr == nOuterCorr-1
			    && corr == nCorr-1
			    && nonOrth == nNonOrthCorr
			)
			{
			    pEqn.solve(mesh.solutionDict().solver("pFinal"));
			}
			else
			{
			    pEqn.solve();
			}

			if (nonOrth == nNonOrthCorr)
			{
			    phi_ == pEqn.flux();
			}
		    }
		}
		else
		{
    phi_ = fvc::interpolate(rho_)
        *((fvc::interpolate(U_) & mesh.Sf()) - fvc::meshPhi(rho_, U_));

		        adjustPhi(phi_, U_, p_);

		    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
		    {
			// Pressure corrector
			fvScalarMatrix pEqn
			(
			    fvm::ddt(psi_, p_)
			  + fvc::div(phi_)
			  - fvm::laplacian(rho_*rUA, p_)
			);

			if
			(
			    oCorr == nOuterCorr-1
			 && corr == nCorr-1
			 && nonOrth == nNonOrthCorr
			)
			{
			    pEqn.solve(mesh.solutionDict().solver("pFinal"));
			}
			else
			{
			    pEqn.solve();
			}

			if (nonOrth == nNonOrthCorr)
			{
			    phi_ += pEqn.flux();
			}
		    }
		}

				//rhoEqn.H begins
			{
			    solve(fvm::ddt(rho_) + fvc::div(phi_));
			}
				// rhoEqn.H ends
	
			// compressibleContinuityErrs.H begins
			{
			    dimensionedScalar totalMass = fvc::domainIntegrate(rho_);

			    sumLocalContErr =
				(fvc::domainIntegrate(mag(rho_ - thermo_.rho()))/totalMass).value();

			    globalContErr =
				(fvc::domainIntegrate(rho_ - thermo_.rho())/totalMass).value();

			    cumulativeContErr += globalContErr;

			    Info<< "time step continuity errors : sum local = " << sumLocalContErr
				<< ", global = " << globalContErr
				<< ", cumulative = " << cumulativeContErr
				<< endl;
			}
			// compressibleContinuityErrs.H ends

			//if (oCorr != nOuterCorr-1)
			{
			    // Explicitly relax pressure for momentum corrector
			    p_.relax();

			    rho_ = thermo_.rho();
			    Info<< "rho max/min : " << max(rho_).value()
				<< " " << min(rho_).value() << endl;
			    rho_.relax();
			    Info<< "rho max/min : " << max(rho_).value()
				<< " " << min(rho_).value() << endl;
			}

			U_ -= rUA*fvc::grad(p_);
			U_.correctBoundaryConditions();

			DpDt_ = fvc::DDt(surfaceScalarField("phiU", phi_/fvc::interpolate(rho_)), p_);

			bound(p_, pMin_);
				//PEqn Ends

            }

            turbulence_->correct();
        }

}

The code which I doubt is

	   if(mesh.moving())
	    {
		// Make the fluxes relative
		phi_ -= fvc::interpolate(rho_)*fvc::meshPhi(rho_, U_);
	    }

Normally in solvers which has mesh motion option, correctPhi.H will be included in the code to correct the fluxes. But I think here that is done internally in the code where I am not clear enough. Could anyone enlighten me on this to check whether I have implemented it correctly or not! I just followed the already implemented icoFlow model in order to create my model.

Any help would be appreciated.

Thanks,
Karthick