|
[Sponsors] |
Solving a coupled linear PDE using modified icoFoam solver |
![]() |
|
LinkBack | Thread Tools | Search this Thread | Display Modes |
![]() |
![]() |
#1 |
Member
Uttam
Join Date: May 2020
Location: Southampton, United Kingdom
Posts: 35
Rep Power: 6 ![]() |
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; } 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. |
|
![]() |
![]() |
![]() |
Tags |
adjoint equations, coupled pde, icofoam, simple algorithm |
|
|
![]() |
||||
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 |