[danylo]

Hello everyone,


I am trying to model a fluid to solid phase change with mesh motion caused by the phase change. Now the Mesh deforms with the velocity , but only in z direction. Therefore it should be mathematically possible to solve the following continuity equation for :





where

= volume fraction of solid
= density of liquid solution
= density of solid
= fluid velocity without Mesh velocity


Now I can not solve this scalar Equation for a vector, but since only the z-component of can be non-zero, I can write it like this:





Mathematically that should make sense, but I lack the OpenFoam tools and knowledge to solve this in my program. Since there is no fvm::ddz method. I know that I can get the derivation of the z-component in the z direction by using:

volScalarField Umz = Um.component(2);
volVectorField gradUmz = fvc::grad(Umz);
volScalarField gradUmzz = gradUmz.component(2);

But I have no idea how to incorporate this in my fvScalarMatrix, so that it solves for Umz.


What I also tried was defining an alphaz vector field, which is 0 in his x and y component and has the alpha values in the z component. I calculated the dot product with the face normals, the same way phi gets calculated, to get a surfaceScalarField, and tried to solve for Umz. This is how it looks:

forAll(alphaz,i)
{
    alphaz[i].component(2) = alphaice[i];
}

forAll(alphaz.boundaryField(), patchI)
   {
   forAll(alphaz.boundaryField()[patchI], faceI)
      {
          alphaz.boundaryFieldRef()[patchI][faceI].component(2) = alphaice.boundaryField()[patchI][faceI];
      }
   }


surfaceScalarField surfacealphaz
(
    IOobject
    (
        "phialphaz",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    linearInterpolate(alphaz) & mesh.Sf()
);





tmp<fvScalarMatrix> tUmEqn
(
    fvc::ddt(rhomix,alphawater) 
  + rhoice*fvc::ddt(alphaice)
  + fvc::div(rhomixUrel)
  + rhoice*fvm::div(surfacealphaz,Umz)
);

The solver compiles, but I can not get it to run on my case. I get following error:


#0 Foam::error::printStack(Foam::Ostream&) at ??:?
#1 Foam::sigFpe::sigHandler(int) at ??:?
#2 ? in "/lib/x86_64-linux-gnu/libc.so.6"
#3 Foam::symGaussSeidelSmoother::smooth(Foam::word const&, Foam::Field<double>&, Foam::lduMatrix const&, Foam::Field<double> const&, Foam::FieldField<Foam::Fie


I suppose my method is not very physical, hence a lot of problems occurr.


If anyone knows, how to incorporate a derivative in only one direction of space in the solver Matrix, that would be great. Any help would be appreciated.


Thank you

[danylo]

Does anyone have any idea how i could solve this issue? I did not make any progress yet

Thank you

[danylo]

Update:
I can get the code for the alphaz surfaceScalarField to run but I get nonsense values. Problem was that alphaice was 0 at the start, and so it was impossible to solve for Um. I set alphaice to a small starting value, but the value Um just spirals out of control up to extremely high values.