CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Running, Solving & CFD (http://www.cfd-online.com/Forums/openfoam-solving/)
-   -   Changing patch temperature directly from the solver (buoyantPimpleFoam) (http://www.cfd-online.com/Forums/openfoam-solving/91109-changing-patch-temperature-directly-solver-buoyantpimplefoam.html)

kriztof August 1, 2011 04:22

Changing patch temperature directly from the solver (buoyantPimpleFoam)
 
Hi,

I would like to be able to manually change the temperature of a patch at each time step directly from the solver. I use buoyantPimpleFoam.
I can't use any available type of a boundary condition, because at each time step I have to solve the heat equation (I do it with the finite difference method).
Unfortunately, I do not know how to force the solver to overwrite its temperature with my own. The BC type of a patch is set to fixedValue at the moment (I had not any better idea).
I know how to change the temperature of the internal field - I can do it by changing the enthalpy. However when I am changing the enthalpy of my patch with the following code:

Code:

forAll(h.boundaryField()[PatchID], i)
{
    h.boundaryField()[PatchID][i] = 295000; // new enthalpy magnitude
}

... it has no effect in the temperature, which is still fixed.
On the other hand, when I try to change directly the patch temperature, by:

Code:

forAll(T.boundaryField()[PatchID], i)
{
    T.boundaryField()[PatchID][i] = 300; // new temperature magnitude
}

... I get the following compilation error:

Code:

error: assignment of read-only location ‘((const Foam::fvPatchField<double>*)
((const Foam::GeometricField<double, Foam::fvPatchField,
Foam::volMesh>::GeometricBoundaryField*)((const Foam::volScalarField*)T)-
>Foam::GeometricField<Type, PatchField, GeoMesh>::boundaryField [with Type =
double, PatchField = Foam::fvPatchField, GeoMesh = Foam::volMesh]())-
>Foam::GeometricField<double, Foam::fvPatchField,
Foam::volMesh>::GeometricBoundaryField::
<anonymous>.Foam::FieldField<Foam::fvPatchField, double>::
<anonymous>.Foam::PtrList<T>::operator[] [with T = Foam::fvPatchField<double>]
(window_S_ID))->Foam::fvPatchField<double>::<anonymous>.Foam::Field<double>::
<anonymous>.Foam::List<double>::<anonymous>.Foam::UList<T>::operator[] [with T =
double](i)’

Any ideas would be appreciated...

newOFuser January 3, 2012 14:32

Hi Krzysztof

Were you able to resolve this?

-Amit



kriztof January 16, 2012 08:26

Dear Amit,

yes, I solved the problem.

What I had to do was to directly change the patch temperature in the following way:

Code:

       
// create an alias for the temperature of patchi
fvPatchScalarField& Tnew = const_cast<fvPatchScalarField&>(thermo.T().boundaryField()[patchi]);

forAll(Tnew, facei)
{
    Tnew[facei] = your_desired_temperature;
}

It is important to use the const_cast, since the temperature is declared as a "const" field.

Regards,
Krzysztof A.

newOFuser January 16, 2012 12:34

Hi Krzysztof,

Thanks so much for the reply! I tried using it (by using a boundary patch), and it works!

However, I'd like to now define the patch inside the mesh region, so that I can specify the temperature on the patch (inside the solver code) for a few time steps, to cause ignition.

Is it possible to define this internal patch in the blockMesh file? If so, how can this be done? (This patch is not a boundary)

Many thanks

kriztof January 16, 2012 14:39

I do not know how to add a patch inside the mesh. I thought, that patches can be defined only at the boundaries.
However if you will find the solution, you should be able to change its temperature in a similar way.

But isn't it better for you to change the temperature of some certain cells in your case?

If so, then just somehow get the IDs of cells you want to manipulate,
and create a rule similar to the following one:

Code:

volScalarField& Tnew = const_cast<volScalarField&>(thermo.T());

forAll(Tnew, celli)
{
    if(celli == proper_cell_ID)
    {
        Tnew[celli] = your_desired_temperature;
    }
}

Of course you should define your own condition for the cell IDs.
I have not tested this, but it should work.

Does it solve your problem?

newOFuser January 23, 2012 17:17

Thanks so much Krzysztof,

I was able to run the code, and it seems to work such that the code indicates the temperature Tnew in the chosen cells (inside the mesh too) has been changed to T_ignition.

But when I plot the results in paraview for T, I notice that the temperature has been affected only at the boundary patches. Any idea why this could be occuring?

Thanks again,
ak

Code:

volScalarField& Tnew = const_cast<volScalarField&>(thermo.T());
forAll(ign.sites(), i)
{
const ignitionSite& ignSite = ign.sites()[i];
if (ignSite.igniting())
{
forAll(ignSite.cells(), icelli)
{
label ignCell = ignSite.cells()[icelli];
Info << "Igniting cell " << ignCell << endl;
Tnew[ignCell] = ignitionT;
Info << "Temperature in cell " << Tnew[ignCell] << endl;
}
}
}

Sample simulation result:

Found ignition cell 8897
location 0.0504707 0.0182951 0.028575 // This is inside region of the mesh

later..

Igniting cell 8897
Temperature in cell 1500 //Tnew changed

But max T shows up as less than 1500 (which could be due to subsequent computations), and only the temperature in the boundary patch is affected.
T gas min/max = 293, 881.926
T max is seen only at boundary patch (inside it is still 293K)

newOFuser January 23, 2012 18:11

I looked into it further. The temperature in the cells does change but subsequent enthalpy calculations seem to bring it down immediately.
Changing the enthalpy instead seems to result in higher temperatures that sustain for longer time duration, as required.

kriztof January 23, 2012 18:19

Oh...

I would guess the problem is that the solver solves the energy equation for enthalpy, instead for temperature (file: "hEqn.H"). The temperature must be subsequently derived from the enthalpy at each time step.

Maybe you should try to change the enthalpy instead of the temperature?
You should be able to do it using thermo.h() instead of thermo.T().

However I don't know, how it will affect the stability of the solution... (this doubt regards any manual changes of the solution during the simulation, not only the enthalpy :) )


--- EDIT ---

I see that you were faster... ;)

newOFuser January 23, 2012 18:35

Hi Krzysztof,

That is right. The reason I wanted to avoid changing enthalpy is that it makes the temperature rise very fast (and I get janaf error).

Previously I was setting the internal field temperature at start of simulation (in 0/T directory) as T_ignition, and it worked better (than changing enthalpy), and so I was looking at ways to specify temperature in specific region, for causing ignition.
(I want ignition at later time and in a specific region now, rather than specify T in entire internal domain at start for ignition)

Per your suggestions, I am able to do it now, but the enthalpy equation is causing issues again! Not sure why it worked fine when 0/T file was modified though. Any thoughts?

Regards,
ak

kriztof January 23, 2012 18:49

Quote:

Originally Posted by newOFuser (Post 340744)
and so I was looking at ways to specify temperature in specific region, for causing ignition.

Maybe the setFields utility?
http://www.openfoam.org/docs/user/da...#x7-530002.3.3

If you do not want to start the ignition at the beginning, you can calculate some preliminary time without the ignition, i.e. 5s.
Then you can manually change the temperature field in the "5" directory.
Then continue your simulation using startTime=5 in the controlDict.

Does it make any sense? I hope so!

newOFuser January 23, 2012 18:52

Thanks! I'll look into the possible ways, and post once I've acceptable results!

Cheers,
ak

kriztof January 23, 2012 18:55

Good luck :)


All times are GMT -4. The time now is 03:11.