[jeroendemto]

Hello everyone,

I'm trying to simulate a 1-d stefan problem, but I'm running into some trouble.
The formula I use to implement phase change can be seen below.

fvScalarMatrix TEqn
(
fvc::ddt(rho,h)
+ fvm::Sp((heatCap*rho)/runTime.deltaT(),T)
- rho*heatCap*TOld/runTime.deltaT()
- fvm::laplacian(thermalCond, T)
);

This works fine if I use a constant density defined as a property i.e. rho = 1000.


However if I use a volScalarField as rho as seen below, the solution diverges.

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Version:  7
     \\/     M anipulation  |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       volScalarField;
    object      rho;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

dimensions       [1 -3 0 0 0 0 0];

internalField    uniform 1000;

boundaryField
{
    inlet
    {
        type        	fixedValue;
        value       	$internalField;
    }
    
    outlet
    {
        type            fixedValue;
        value           $internalField;
    }

    walls
    {
    	type		fixedValue;
    	value		$internalField;
    }

    sides
    {
        type        empty;
    }

}

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

I would think this would work the same but it does not, how can I fix this?

kind regards,
Jeroen