[Aicharem]

Hi Foamers,

I am working on creating a model of an EHD thruster. I have posted about this before, when trying to create a steady state solver, however I now believe that my solver must be a transient one. I made progress with simpleFoam as a base solver for the steady state model, however when trying to use pisoFoam as a base I am not doing so well. Both solvers use code from electrostaticFoam to solve the electrostatic problem. Variables from electrostaticFoam have the letter E attached to them so that they do not overlap with the fluid solver variables.

My aim is to follow a similar approach to the ‘adding temperature to icoFoam’ tutorial, to add in the equations for Gauss’s Law and charge conservation. These results will then be implemented into the UEqn.H file to change the flow parameters.

I currently have a model that will run, but is (a) extremely slow and (b) produces floating point exceptions before any reasonable results can be reached. It would be a huge help if someone could look over the code created in my solver and tell me if the solver is correct and it is the boundary conditions that are causing the problem, or vice versa.

The steps for creating the solver are as follows:

• In a copy of the pisoFoam source code, change the names of pisoFoam.C and alter ‘files’ in the ‘make’ folder to create pisoThrust.

• In ‘createFields.H’, add instruction to read the ‘physicalProperties’ file, create the necessary constants (epsilon0, kE, dens, DE, defined as dimensioned scalars) and create the fields rhoE, phiE and rhoEFlux.

• In ‘pisoThrust.C’ add the electrostatic equations (added terms in red):
  

      while (runTime.loop())
      {
          Info<< "Time = " << runTime.timeName() << nl << endl;
  
          #include "CourantNo.H"
  
          // Pressure-velocity PISO corrector
          {
  
          solve
          (
              fvm::laplacian(phiE) + rhoE/epsilon0  \\Gauss's Law
          );
  
          rhoEFlux = -kE*mesh.magSf()*fvc::snGrad(phiE); \\Charge conservation
  
          solve
          (
              fvm::ddt(rhoE) + fvm::div(rhoEFlux, rhoE) + fvm::div(phi, rhoE) \\Charge conservation pt 2
          );
  
              #include "UEqn.H"
  
              // --- PISO loop
              while (piso.correct())
              {
                  #include "pEqn.H"
              }
          }
  
          laminarTransport.correct();
          turbulence->correct();
  
          runTime.write();
  
          Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
              << "  ClockTime = " << runTime.elapsedClockTime() << " s"
              << nl << endl;
      }

• In UEqn.H add the electrostatic force as a variable body force:
  

  fvVectorMatrix UEqn
  (
      fvm::ddt(U) + fvm::div(phi, U)
    + MRF.DDt(U)
    + turbulence->divDevReff(U)
    + fvc::grad(rhoE/dens*phiE)
   ==
      fvOptions(U)
  );

The solver is then compiled and is run on test cases. I have attached the full source code. If I have made any obvious errors regarding the placement of equations within different loops please can you let me know. I am still very inexperienced with OpenFOAM.

Thanks,
Aicharem