CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   Adding an internal heat source for compressibleTwoPhaseEulerFoam (https://www.cfd-online.com/Forums/openfoam-programming-development/123155-adding-internal-heat-source-compressibletwophaseeulerfoam.html)

Kenna September 5, 2013 11:10
Adding an internal heat source for compressibleTwoPhaseEulerFoam
 
1 Attachment(s)
Hello dear Foamers,

I am working on the simulation of a plasma-arc-welding process and I am using the compressibleTwoPhaseEulerFoam Solver now, because it calculates the two phases (I use Argon and Air) and the temperature which I am looking for mainly.
But to be exact enough I need a heat source, where the electrical current is flowing through the plasma.
I would try to create a new Field "q" which would describe the increasing energy in the cells relative to the geometrical size of the cell.

I found already an interesting post like this one:
http://www.cfd-online.com/Forums/ope...ionfoam-2.html
And this one:
http://www.cfd-online.com/Forums/ope...ew-solver.html
But the solver is still different. I don't have a TEqn.H. I guess I need to change the EEqn.H, but I don't know how to do it in detail.
I don't know much about C++ and the "Programmers Guide" didn't enlighten me much too.

My plan would be:
- a readable heat source file for the Mesh [Joule/Second/m³]
- accessing the volumetric size of the cells [m³]
- calculating the total additional energy of the cell [Joule/Second/m³]*[m³]*[StepTimeSeconds]
- subtracting it in the right part of the following equation

Code:

    fvScalarMatrix he1Eqn
    (
        fvm::ddt(alpha1, he1) + fvm::div(alphaPhi1, he1)
      + fvc::ddt(alpha1, K1) + fvc::div(alphaPhi1, K1)

        // Compressibity correction
      - fvm::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), he1)
      - (fvc::ddt(alpha1) + fvc::div(alphaPhi1))*K1

      + (
            he1.name() == "e"
          ? fvc::div(alphaPhi1, p)
          : -dalpha1pdt
        )/rho1

      - fvm::laplacian(k1, he1)
    ==
        heatTransferCoeff*(thermo2.T() - thermo1.T())/rho1
      + heatTransferCoeff*he1/Cpv1/rho1
      - fvm::Sp(heatTransferCoeff/Cpv1/rho1, he1)
    );
Do I assume right, if I say this is the important part of the EEqn.H?

Also I attached a slice of my 3D-model.



I would appreciate your help a lot.
Sincerely Kenna

Kenna September 16, 2013 10:08
Current try
 
Well, so far I could create a Qp (heat flux) field with this dimensions (0, 2, -3, 0, 0, 0, 0).

I put a Qp file in the latest Time folder and put where the electric arc should be, via setFieldsDict in a cylindrical Form the value 2.14748e+09. The rest stays zero.

in the createFields.H from the solver I entered this:

Code:
    Info<< "Reading field Qp\n" << endl;
    volScalarField Qp
    (
        IOobject
        (
            "Qp",
            mesh.time().timeName(),
            //runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh,
    dimensionedScalar("Qp", dimensionSet(0, 2, -3, 0, 0, 0, 0), 0)
Qp is Energy/kg/second or m²/s³

I added Qp very simple in the EEqn.
Code:
    fvScalarMatrix he1Eqn
    (
        fvm::ddt(alpha1, he1) + fvm::div(alphaPhi1, he1)
      + fvc::ddt(alpha1, K1) + fvc::div(alphaPhi1, K1)

        // Compressibity correction
      - fvm::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), he1)
      - (fvc::ddt(alpha1) + fvc::div(alphaPhi1))*K1

      + (
            he1.name() == "e"
          ? fvc::div(alphaPhi1, p)
          : -dalpha1pdt
        )/rho1

      - fvm::laplacian(k1, he1)
    ==
        heatTransferCoeff*(thermo2.T() - thermo1.T())/rho1
      + heatTransferCoeff*he1/Cpv1/rho1
      - fvm::Sp(heatTransferCoeff/Cpv1/rho1, he1)
      + Qp// neu
    );
The Qp file looks like this:
Code:
dimensions      [0 2 -3 0 0 0 0];

internalField  nonuniformFixedValue //List<scalar>
78072
(
2.14748e+09
2.14748e+09
2.14748e+09
2.14748e+09
2.14748e+09
2.14748e+09
//etc
0
0
0
);
}
    item
    {
        type            zeroGradient;
    }
//etc
But it seems like the Qp file isn't even considered in the calculation and runs without the file (I used IOobject::MUST_READ ?) the solver writes every time step a new Qp file. But the new generated file contains as values only zeros like this:

Code:
dimensions      [0 2 -3 0 0 0 0];

internalField  uniform 0;

boundaryField
{
    inlet1
    {
        type            calculated;
        value          uniform 0;
    }
//etc
I guess the main problem is, Qp is always zero, because it doesn't get read in.

jherb September 16, 2013 11:30
Could you try without the red line (and the comma at the line above after "mesh").

Quote:
Originally Posted by Kenna (Post 452033)
in the createFields.H from the solver I entered this:

Code:
    Info<< "Reading field Qp\n" << endl;
    volScalarField Qp
    (
        IOobject
        (
            "Qp",
            mesh.time().timeName(),
            //runTime.timeName(),
            mesh,
            IOobject::MUST_READ,
            IOobject::AUTO_WRITE
        ),
        mesh,
        dimensionedScalar("Qp", dimensionSet(0, 2, -3, 0, 0, 0, 0), 0)
  )

Kenna September 17, 2013 03:51
Quote:
Originally Posted by jherb (Post 452047)
Could you try without the red line (and the comma at the line above after "mesh").
Thanks, this works!

The temperature is rising in the right area. The pressure doesn't change as expected, because it is an open system.

So far I have no further question :)


All times are GMT -4. The time now is 07:36.