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

Solving a coupled linear PDE using modified icoFoam solver

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

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   June 28, 2022, 05:36
Unhappy Solving a coupled linear PDE using modified icoFoam solver
  #1
New Member
 
Uttam
Join Date: May 2020
Location: Southampton, United Kingdom
Posts: 20
Rep Power: 4
openfoam_aero is on a distinguished road
Dear all,

OpenFOAM version: OpenFOAM 7

I am currently working on data assimilation and I have to solve a set of adjoint equations for a domain consisting of a flow past a cylinder. The thing is, these equations are similar to NS equations but with the exception that they are not non-linear. I present the equattions below

adjoint_equation.png

In these equations, I have the value of \Tilde{u} and the term on the RHS (which is a source term in this case). I need to solve for u_dagger and p_dagger. There is a continuity equation that is also present given by

continuity.png

Now I use the icoFoam solver by removing the time derivative (since it is steady state problem). I am presenting the code below

Code:
#include "fvCFD.H"
#include "simpleControl.H"

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

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

    simpleControl simple(mesh);

    #include "createFields.H"
    #include "initContinuityErrs.H"

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

    Info<< "\nStarting time loop\n" << endl;

    while (simple.loop(runTime))
    {
        Info<< "Time = " << runTime.timeName() << nl << endl;

        //#include "CourantNo.H"

        // Momentum predictor


	tmp<fvVectorMatrix> tUEqn
	(
		 - (fvm::laplacian(nu,UDagger)) +(U1 & T(fvc::grad(UDagger))) - (UDagger & fvc::grad(U1)) - dEdU
	
	);

        
        fvVectorMatrix& UEqn = tUEqn.ref();

    	UEqn.relax();


        solve(UEqn == -fvc::grad(pDagger));
        


        volScalarField rAU(1.0/UEqn.A());
        
        volVectorField HbyA(constrainHbyA(rAU*UEqn.H(),UDagger, pDagger));
        
        surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA));
	
	adjustPhi(phiHbyA, UDagger, pDagger);

	
	tmp<volScalarField> rAtU(rAU);

    	if (simple.consistent())
    		{
		rAtU = 1.0/(1.0/rAU - UEqn.H1());
		phiHbyA +=
		    fvc::interpolate(rAtU() - rAU)*fvc::snGrad(pDagger)*mesh.magSf();
		HbyA -= (rAU - rAtU())*fvc::grad(pDagger);
    		}
        
        tUEqn.clear();
        


        // Update the pressure BCs to ensure flux consistency
        constrainPressure(pDagger, UDagger, phiHbyA, rAU);

        // Non-orthogonal pressure corrector loop
        //while (simple.correctNonOrthogonal())
        //{
             // Pressure corrector

             fvScalarMatrix pEqn
             (
                 fvm::laplacian(rAU, pDagger) == fvc::div(phiHbyA)
             );

 //            pEqn.setReference(pRefCell, pRefValue);

             pEqn.solve();

             //if (simple.finalNonOrthogonalIter())
             //{
                 phi = phiHbyA - pEqn.flux();
             //}
         //}

         #include "continuityErrs.H"
            
         pDagger.relax();

         UDagger = HbyA - rAU*fvc::grad(pDagger);
         UDagger.correctBoundaryConditions();

         runTime.write();

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

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

    return 0;
}
The boundary conditions for U_Dagger and p_dagger are similar to the flow past a cylinder case with symmetry BCs on the top and bottom, empty BCs for front and back, noSlip for the cylinder, inlet BC for velocity and zeroGradient outlet for velocity.

I use a GAMG solver for p_dagger and PBiCG solver for U_Dagger with residuals of 10^-4 for each.

I get a U_Dagger field after about 600 iterations but it is not correct when compared to the reference case from literature. The gradients close to the cylinder wall are not present in my case. I tried debugging a lot and looked through a lot of options.

My questions
1. Am I setting up the equations correctly?
2. Am I using the right method of solving these equations? Is there a simpler alternative?

I suspect that the issue might be with fvVectorMatrix tUEqn. Any suggestions on the problem will be appreciated.
P.S i have included #createPhi1.H at the end of my #createFields.H file and the #createPhi1.H file looks like this

Code:
Info<< "Reading/calculating face flux field phi\n" << endl;

surfaceScalarField phi
(
    IOobject
    (
        "phi",
        runTime.timeName(),
        mesh,
        IOobject::READ_IF_PRESENT,
        IOobject::AUTO_WRITE
    ),
    fvc::flux(UDagger)
);
__________________
Best Regards
Uttam

-----------------------------------------------------------------

“When everything seem to be going against you, remember that the airplane takes off against the wind, not with it.” – Henry Ford.
openfoam_aero is offline   Reply With Quote

Reply

Tags
adjoint equations, coupled pde, icofoam, simple algorithm

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
Help sought on axial compressor simulation jyotir OpenFOAM Running, Solving & CFD 0 November 17, 2021 10:49
Floating point exception error lpz_michele OpenFOAM Running, Solving & CFD 53 October 19, 2015 02:50
SLTS+rhoPisoFoam: what is rDeltaT??? nileshjrane OpenFOAM Running, Solving & CFD 4 February 25, 2013 04:13
Orifice Plate with a fully developed flow - Problems with convergence jonmec OpenFOAM Running, Solving & CFD 3 July 28, 2011 05:24
IcoFoam parallel woes msrinath80 OpenFOAM Running, Solving & CFD 9 July 22, 2007 02:58


All times are GMT -4. The time now is 22:26.