はじめに

OpenFOAMのソルバーinterFoamで使用する表面張力モデルはCSF(Continuum Surface Force)モデルである。このモデルではスプリアス・カレントの生成により表面張力の精度が低下する問題が知られている。ESI版OpenFOAM v2212からはSmoothed CSFモデルが使用できるようになり、このモデルを使用すればスプリアス・カレントの生成を抑制することができる。本稿ではこれらのモデルに加えて、SSF(Sharp Surface Force)モデルの実装について説明する。
モデル式などは次の文献を参考にした。

Vachaparambil, K. J., Einarsrud, K. E. 2019a. Comparison of surface tension models for the volume of fluid method. Processes, 7: 542.

なお、実装したモデルの動作確認は、次のページに記載する。

作業環境

OpenFOAMのバージョンは文献に合わせてFoundation版の6とした。OpenFOAM6の動作が保証されているUbuntu 18.04.6 LTSをWindows11のWSL2にインストールした。そしてOpenFOAM6のソースコードをUbuntu上に展開し、コードの改造およびコンパイルができるようにした。

対象となるソースコード

表面張力モデルが実装されているソースコードは、$FOAM_SRC/transportModels/interfacePropertiesディレクトリにあるinterfaceProperties.HおよびinterfaceProperties.Cである。

表面張力モデル

概要

表面張力は次の量の積で計算される。

表面張力
曲率
界面の水の体積分率の勾配

表面張力は、計算を実行するディレクトリの$constant/transportProperties$で$sigma$として設定された値を用いる。
曲率は、対象となるソースコードの$interfaceProperties$クラスのプライベートメンバ関数である$calculateK()$で計算される。
界面の水の体積分率の勾配は、界面に垂直な勾配を用いる。
以下では各表面張力モデルでの計算方法を記す。

曲率の計算

セル中心の水の体積分率の勾配の計算

CSFモデル

与えられた水の体積分率$\alpha_{1}$($alpha.water$のこと)より、セル中心の水の体積分率の勾配を計算する。

(\nabla\alpha)_c=\nabla\alpha_{1}\tag{1}

ここで、添字のcはセル中心を表す。
実装は次の通り。

csf_tgradAlpha

        tmp<volVectorField> tgradAlpha;
        tgradAlpha = fvc::grad(alpha1_, "nHat");

Smoothed CSFモデル

与えられた水の体積分率$\alpha_{1}$をsmoothingして、セル中心の水の体積分率の勾配を計算する。
あるセルの水の体積分率と隣接するセルのセル中心の水の体積分率の補間値をフェイス中心値とした後、セル周りのフェイス中心の水の体積分率をフェイスの面積で加重平均した値を新たなセル中心値とする。この操作を$nAlphaSmoothCurvature$回行う。

\tilde{\alpha_{1}}^{(i+1)}=\frac{\sum_{f=1}^{N} \langle\tilde{\alpha_{1}}^{(i)}\rangle_{c \rightarrow f} S_{f}}{\sum_{f=1}^{N} S_{f}}=\langle\langle\tilde{\alpha_{1}}^{(i)}\rangle_{c \rightarrow f}\rangle_{f \rightarrow c}\tag{2}

ここで、$\tilde{\alpha_{1}}^{(0)}=\alpha_{1}$とする。添字の$f$はフェイスを表し、$S_{f}$はフェイス$f$の面積を表す。
$\langle \phi_{c} \rangle_{c->f}$はセル中心からフェイス中心へ補間値を設定する操作を表す。
$\langle \phi_{f} \rangle_{f->c}$はフェイス中心の面積の加重平均値をセル中心に設定する操作を表す。

この計算の回数は、$system/fvSolution$の$alpha.water$ディクショナリ内の$nAlphaSmoothCurvature$で設定できるようにする（デフォルト値は0で、この場合はCSFモデルを用いる）。
(2)式で計算した水の体積分率よりセル中心の水の体積分率の勾配を計算する。

(\nabla\alpha)_c=\nabla\tilde{\alpha_{1}}^{(nAlphaSmoothCurvature-1)}\tag{3}

実装は次の通り。

smootheCsf_tgradAlpha

        tmp<volVectorField> tgradAlpha;
        // Smooth interface curvature to reduce spurious currents
        volScalarField alpha1L(alpha1_);
        for (int i = 0; i < nAlphaSmoothCurvature_; ++i)
        {
            alpha1L = fvc::average(fvc::interpolate(alpha1L));
        }
        tgradAlpha = fvc::grad(alpha1L, "nHat");

SSFモデル

与えられた水の体積分率$\alpha_{1}$を3回smoothingする。

\alpha_{s}^{(i+1)}=C\langle\langle \alpha_{s}^{(i)}\rangle_{c \rightarrow f}\rangle_{f \rightarrow c} + (1-C)\alpha_{s}^{(i)}\tag{4}

ここで、$\alpha_{s}^{(0)}=\alpha_{1}, i=0,1,2$とする。
smoothingの重み付け係数$C$は、$system/fvSolution$の$alpha.water$ディクショナリ内の$cSmooth$で設定できるようにする（デフォルト値は$0.5$）。
(4)式で計算した水の体積分率よりセル中心の水の体積分率の勾配を計算する。

(\nabla\alpha)_c=\nabla\alpha_{s}^{(3)}\tag{５}

実装は次の通り。

ssf_tgradAlpha

        tmp<volVectorField> tgradAlpha;
        // SSF Model: Smooth alpha1
        volScalarField alpha1S(alpha1_);
        for (int i = 0; i < 3; ++i)
        {
            alpha1S = cSmooth_ * fvc::average(fvc::interpolate(alpha1S))
                    + (scalar(1) - cSmooth_) * alpha1S;
        }
        tgradAlpha = fvc::grad(alpha1S, "nHat");

フェイス中心の水の体積分率の勾配の計算

セル中心の水の体積分率の勾配よりフェイス中心の水の体積分率の勾配を計算する。
全モデル共通で、次式で計算する。

\vec{(\nabla\alpha)_f}=\langle (\nabla\alpha)_{c}\rangle_{c \rightarrow f}\tag{6}

実装は次の通り。

gradAlphaf

    surfaceVectorField gradAlphaf(fvc::interpolate(tgradAlpha));

フェイス単位の界面法線ベクトルの計算

フェイス中心の水の体積分率の勾配からフェイス単位の界面法線ベクトルを計算する。
全モデル共通で、次式で計算する。

\vec{\hat{n}}_{fv}=\frac{\vec{(\nabla\alpha)_{f}}}{|\vec{(\nabla\alpha)_{f}}|+\delta_{N}}\tag{7}

実装は次の通り。

nHatfv

    surfaceVectorField nHatfv(gradAlphaf/(mag(gradAlphaf) + deltaN_));

$\delta_{N}$は次式で与えられ、実装ではコンストラクタで設定される。

\delta_{N}=\frac{10^{-8}}{(\sum \frac{V_{i}}{N})^{1/3}}\tag{8}

deltaN_

    deltaN_
    (
        "deltaN",
        1e-8/pow(average(alpha1.mesh().V()), 1.0/3.0)
    ),

接触角の補正

今回の改造とは無関係のため省略する。

フェイス単位の界面法線ベクトル束の計算

全モデル共通で、次式でフェイス単位の界面法線ベクトル束を計算する。

\hat{n}_f=\vec{\hat{n}}_{fv} \cdot \vec S_f\tag{9}

ここで、$\vec S_f$はフェイス$f$の面積を法線ベクトルに乗じたベクトルである。
実装は次の通り。

nHatf_

    nHatf_ = nHatfv & Sf;

曲率の計算

次式で曲率を計算する。

\kappa = -\nabla \hat{n}_f\tag{10}

実装は次の通り。

K_

    K_ = -fvc::div(nHatf_);

$fvc::div$の引数$nHatf \_ $は$surfaceScalarField$クラスである。
この場合、$fvc::div(nHatf \_ )$は$fvc::surfaceIntegrate(nHatf \_ )$を呼び出す。
これはそれぞれのセルを区切る面（界面）での値$nHatf \_ $の総和をとり、セルの体積で割るという操作を行う。戻り値は$volScalarField$クラスとなる。曲率$K \_ $は$volscalarField$クラスで定義されている。
CSFモデルおよびSmoothed CSFモデルはこの曲率$K \_ $を使用する。
SSFモデルではこの曲率$K \_ $を基にして平滑化(smoothing)・鋭敏化(sharpening)した曲率$Kfinal \_ $を計算して使用する。この曲率$Kfinal \_ $は$surfaceScalarField$クラスである。

曲率の平滑化・鋭敏化(SSFモデルのみ)

次式で曲率の平滑化・鋭敏化を行う。

{\kappa_s}^{(i+1)}=A\kappa+(1-A)\frac{{\langle {\langle w {\kappa_s}^{(i)} \rangle}_{c \rightarrow f} \rangle}_{f \rightarrow c}}{{\langle {\langle w \rangle}_{c \rightarrow f} \rangle}_{f \rightarrow c}}\tag{11}

ここで、

{\kappa_s}^{(0)}=\kappa\tag{12}

\alpha_c = min(1, max(\alpha_1, 0))\tag{13}

A = 2\sqrt{\alpha_c (1-\alpha_c)}\tag{14}

w = \sqrt{\alpha_c(1-\alpha_c)+10^{-3}}\tag{15}

(11)式を2回計算した結果より、SSFモデルの曲率$\kappa_{final}$を算出する。

\kappa_{final}=\frac{{\langle w {\kappa_s}^{(2)} \rangle}_{c \rightarrow f}}{{\langle w \rangle}_{c \rightarrow f}}\tag{16}

実装は次の通り。

Kfinal_

        const volScalarField alphaC
        (
            "alphaC",
            min(max(alpha1_, scalar(0)), scalar(1))
        );

        const volScalarField A
        (
            "A",
            alphaC * (scalar(1) - alphaC)
        );

        const volScalarField w
        (
            "w",
            sqrt(A + scalar(0.001))
        );
        
        // Smoothing curvature
        volScalarField Ks(K_);
        for (int i = 0; i < 2; ++i)
        {
            Ks = scalar(2) * sqrt(A) * K_ + (scalar(1) - scalar(2) * sqrt(A))
               * fvc::average(fvc::interpolate(w * Ks)) / fvc::average(fvc::interpolate(w));
        }
        
        Kfinal_ = fvc::interpolate(w * Ks) / fvc::interpolate(w);

$Kfinal \_ $はヘッダファイルおよびソースファイルのコンストラクタに追加する。

表面張力の計算

CSFモデルおよびSmoothed CSFモデル

(10)式で計算した曲率および水の体積分率の勾配より表面張力を計算する。

F_{st}={\langle \sigma \kappa \rangle}_{c \rightarrow f}(\nabla \alpha_1)_f\tag{17}

ここで、水の体積分率の勾配は、界面に垂直な勾配である。
実装は次の通りである。

Fst_csf

Foam::tmp<Foam::surfaceScalarField>
Foam::interfaceProperties::surfaceTensionForce() const
{
    return fvc::interpolate(sigmaK())*fvc::snGrad(alpha1_);
}

SSFモデル

SSFモデルでも(17)式と同様に表面張力の計算を行う。だたし、水の体積分率はsharpeningしたものを使用する。

F_{st}={\langle \sigma \rangle}_{c \rightarrow f}\kappa_{final}(\nabla \alpha_{sh})_f\tag{18}

\alpha_{sh}=\frac{1}{1-C_{sh}}\Biggl(min\Biggl[max\Biggl(\alpha_1, \frac{C_{sh}}{2} \Biggr), 1-\frac{C_{sh}}{2} \Biggr] -\frac{C_{sh}}{2}\Biggr)\tag{19}

sharpening係数$C_{sh}$は、$system/fvSolutionのalpha.water$ディクショナリ内の$cSharp$で設定できるようにする（デフォルト値は$0.5$）。
実装は次の通り。

Fst_ssf

Foam::tmp<Foam::surfaceScalarField>
Foam::interfaceProperties::surfaceTensionForce() const
{
    const volScalarField alphaSh
    (
        "alphaSh",
        (
            min(
                max(alpha1_, cSharp_ / scalar(2)),
                scalar(1) - cSharp_ / scalar(2)
            )
            - cSharp_ / scalar(2)
        ) / (scalar(1) - cSharp_)
    );
    return sigmaKfinal()*fvc::snGrad(alphaSh);
}

CSFおよびSmoothed CSF使用するメンバ関数$sigmaK()$に代わり、SSFモデルではメンバ関数$sigmaKFinal()$を使用するようにする。SSFモデルで使用する曲率$Kfinal \_ $は$surfaceScalarField$クラスであるためである。

sigmaKFinal

Foam::tmp<Foam::surfaceScalarField>
Foam::interfaceProperties::sigmaKfinal() const
{
    return fvc::interpolate(sigmaPtr_->sigma())*Kfinal_;
}

まとめ

冒頭にしるした文献に沿ってSSFモデルの実装を行った。Smoothed CSFモデルについてはESI版に導入されているモデルを流用して実装した。最後にソースコードの全体を提示する。

interfaceProperties.H

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

Class
    Foam::interfaceProperties

Description
    Contains the interface properties.

    Properties to aid interFoam:
    -# Correct the alpha boundary condition for dynamic contact angle.
    -# Calculate interface curvature.

SourceFiles
    interfaceProperties.C

\*---------------------------------------------------------------------------*/

#ifndef interfaceProperties_H
#define interfaceProperties_H

#include "IOdictionary.H"
#include "surfaceTensionModel.H"
#include "volFields.H"
#include "surfaceFields.H"

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

namespace Foam
{

/*---------------------------------------------------------------------------*\
                           Class interfaceProperties Declaration
\*---------------------------------------------------------------------------*/

class interfaceProperties
{
    // Private data

        //- Keep a reference to the transportProperties dictionary
        const dictionary& transportPropertiesDict_;

        // 2025/5/19 add start
        //- Number of iterations to smooth interfacial curvatures
        scalar nAlphaSmoothCurvature_;
        //- Smoothen alpha coefficient for SSF model
        scalar cSmooth_;
        //- Sharpen alpha coefficient for SSF model
        scalar cSharp_;
        //- Using SSF Model flag
        bool SSFModel_;
        // 2025/5/19 add end

        //- Compression coefficient
        scalar cAlpha_;

        //- Surface tension
        autoPtr<surfaceTensionModel> sigmaPtr_;

        //- Stabilisation for normalisation of the interface normal
        const dimensionedScalar deltaN_;

        const volScalarField& alpha1_;
        const volVectorField& U_;
        surfaceScalarField nHatf_;
        volScalarField K_;
        // 2025/5/29 add start
        surfaceScalarField Kfinal_;
        // 2025/5/29 add end

    // Private Member Functions

        //- Disallow default bitwise copy construct and assignment
        interfaceProperties(const interfaceProperties&);
        void operator=(const interfaceProperties&);

        //- Correction for the boundary condition on the unit normal nHat on
        //  walls to produce the correct contact dynamic angle
        //  calculated from the component of U parallel to the wall
        void correctContactAngle
        (
            surfaceVectorField::Boundary& nHat,
            const surfaceVectorField::Boundary& gradAlphaf
        ) const;

        //- Re-calculate the interface curvature
        void calculateK();


public:

    //- Conversion factor for degrees into radians
    static const scalar convertToRad;


    // Constructors

        //- Construct from volume fraction field gamma and IOdictionary
        interfaceProperties
        (
            const volScalarField& alpha1,
            const volVectorField& U,
            const IOdictionary&
        );


    // Member Functions

        scalar cAlpha() const
        {
            return cAlpha_;
        }

        const dimensionedScalar& deltaN() const
        {
            return deltaN_;
        }

        const surfaceScalarField& nHatf() const
        {
            return nHatf_;
        }

        tmp<volScalarField> sigmaK() const;
        // 2025/5/29 add start
        tmp<surfaceScalarField> sigmaKfinal() const;
        // 2025/5/29 add end
        tmp<surfaceScalarField> surfaceTensionForce() const;

        //- Indicator of the proximity of the interface
        //  Field values are 1 near and 0 away for the interface.
        tmp<volScalarField> nearInterface() const;

        void correct();

        //- Read transportProperties dictionary
        bool read();
};


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

} // End namespace Foam

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

#endif

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

interfaceProperties.C

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "interfaceProperties.H"
#include "alphaContactAngleFvPatchScalarField.H"
#include "mathematicalConstants.H"
#include "surfaceInterpolate.H"
#include "fvcDiv.H"
#include "fvcGrad.H"
#include "fvcSnGrad.H"
// 2025/5/17 add start
#include "fvcAverage.H"
// 2025/5/17 add en

// * * * * * * * * * * * * * * * Static Member Data  * * * * * * * * * * * * //

const Foam::scalar Foam::interfaceProperties::convertToRad =
    Foam::constant::mathematical::pi/180.0;


// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

// Correction for the boundary condition on the unit normal nHat on
// walls to produce the correct contact angle.

// The dynamic contact angle is calculated from the component of the
// velocity on the direction of the interface, parallel to the wall.

void Foam::interfaceProperties::correctContactAngle
(
    surfaceVectorField::Boundary& nHatb,
    const surfaceVectorField::Boundary& gradAlphaf
) const
{
    const fvMesh& mesh = alpha1_.mesh();
    const volScalarField::Boundary& abf = alpha1_.boundaryField();

    const fvBoundaryMesh& boundary = mesh.boundary();

    forAll(boundary, patchi)
    {
        if (isA<alphaContactAngleFvPatchScalarField>(abf[patchi]))
        {
            alphaContactAngleFvPatchScalarField& acap =
                const_cast<alphaContactAngleFvPatchScalarField&>
                (
                    refCast<const alphaContactAngleFvPatchScalarField>
                    (
                        abf[patchi]
                    )
                );

            fvsPatchVectorField& nHatp = nHatb[patchi];
            const scalarField theta
            (
                convertToRad*acap.theta(U_.boundaryField()[patchi], nHatp)
            );

            const vectorField nf
            (
                boundary[patchi].nf()
            );

            // Reset nHatp to correspond to the contact angle

            const scalarField a12(nHatp & nf);
            const scalarField b1(cos(theta));

            scalarField b2(nHatp.size());
            forAll(b2, facei)
            {
                b2[facei] = cos(acos(a12[facei]) - theta[facei]);
            }

            const scalarField det(1.0 - a12*a12);

            scalarField a((b1 - a12*b2)/det);
            scalarField b((b2 - a12*b1)/det);

            nHatp = a*nf + b*nHatp;
            nHatp /= (mag(nHatp) + deltaN_.value());

            acap.gradient() = (nf & nHatp)*mag(gradAlphaf[patchi]);
            acap.evaluate();
        }
    }
}


void Foam::interfaceProperties::calculateK()
{
    const fvMesh& mesh = alpha1_.mesh();
    const surfaceVectorField& Sf = mesh.Sf();

    // Cell gradient of alpha
    // 2025/5/19 mod start
    // const volVectorField gradAlpha(fvc::grad(alpha1_, "nHat"));
    tmp<volVectorField> tgradAlpha;
    //if (nAlphaSmoothCurvature_ < 1)
    if (SSFModel_)
    {
        // SSF Model: Smooth alpha1
        volScalarField alpha1S(alpha1_);
        for (int i = 0; i < 3; ++i)
        {
            alpha1S = cSmooth_ * fvc::average(fvc::interpolate(alpha1S))
                    + (scalar(1) - cSmooth_) * alpha1S;
        }
        tgradAlpha = fvc::grad(alpha1S, "nHat");
    }
    else if (nAlphaSmoothCurvature_ < 1)
    {
        tgradAlpha = fvc::grad(alpha1_, "nHat");
    }
    else
    {
        // Smooth interface curvature to reduce spurious currents
        volScalarField alpha1L(alpha1_);
        for (int i = 0; i < nAlphaSmoothCurvature_; ++i)
        {
            alpha1L = fvc::average(fvc::interpolate(alpha1L));
        }
	tgradAlpha = fvc::grad(alpha1L, "nHat");
    }
    // 2025/5/19 mod end

    // Interpolated face-gradient of alpha
    // 2025/5/19 mod start
    // surfaceVectorField gradAlphaf(fvc::interpolate(gradAlpha));
    surfaceVectorField gradAlphaf(fvc::interpolate(tgradAlpha));
    // 2025/5/19 mod end

    // gradAlphaf -=
    //    (mesh.Sf()/mesh.magSf())
    //   *(fvc::snGrad(alpha1_) - (mesh.Sf() & gradAlphaf)/mesh.magSf());

    // Face unit interface normal
    surfaceVectorField nHatfv(gradAlphaf/(mag(gradAlphaf) + deltaN_));
    // surfaceVectorField nHatfv
    // (
    //     (gradAlphaf + deltaN_*vector(0, 0, 1)
    //    *sign(gradAlphaf.component(vector::Z)))/(mag(gradAlphaf) + deltaN_)
    // );
    correctContactAngle(nHatfv.boundaryFieldRef(), gradAlphaf.boundaryField());

    // Face unit interface normal flux
    nHatf_ = nHatfv & Sf;

    // Simple expression for curvature
    K_ = -fvc::div(nHatf_);

    // Complex expression for curvature.
    // Correction is formally zero but numerically non-zero.
    /*
    volVectorField nHat(gradAlpha/(mag(gradAlpha) + deltaN_));
    forAll(nHat.boundaryField(), patchi)
    {
        nHat.boundaryField()[patchi] = nHatfv.boundaryField()[patchi];
    }

    K_ = -fvc::div(nHatf_) + (nHat & fvc::grad(nHatfv) & nHat);
    */
    // 2025/5/19 add start
    // SSF Model: Smooth curvature
    if (SSFModel_)
    {
        const volScalarField alphaC
        (
            "alphaC",
            min(max(alpha1_, scalar(0)), scalar(1))
        );

        const volScalarField A
        (
            "A",
            alphaC * (scalar(1) - alphaC)
        );

        const volScalarField w
        (
            "w",
            sqrt(A + scalar(0.001))
        );

        volScalarField Ks(K_);
        for (int i = 0; i < 2; ++i)
        {
            Ks = scalar(2) * sqrt(A) * K_ + (scalar(1) - scalar(2) * sqrt(A))
               * fvc::average(fvc::interpolate(w * Ks)) / fvc::average(fvc::interpolate(w));
        }

        Kfinal_ = fvc::interpolate(w * Ks) / fvc::interpolate(w);
    }
    // 2025/5/19 add end
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::interfaceProperties::interfaceProperties
(
    const volScalarField& alpha1,
    const volVectorField& U,
    const IOdictionary& dict
)
:
    transportPropertiesDict_(dict),
    // 2025/5/19 add start
    nAlphaSmoothCurvature_
    (
        alpha1.mesh().solverDict(alpha1.name()).lookupOrDefault
        (
            "nAlphaSmoothCurvature", 0
        )
    ),
    cSmooth_
    (
        alpha1.mesh().solverDict(alpha1.name()).lookupOrDefault
        (
            "cSmooth", 0.5
        )
    ),
    cSharp_
    (
        alpha1.mesh().solverDict(alpha1.name()).lookupOrDefault
        (
            "cSharp", 0.5
        )
    ),
    SSFModel_
    (
        alpha1.mesh().solverDict(alpha1.name()).lookupOrDefault
        (
            "SSFModel", false
        )
    ),
    // 2025/5/19 add end
    cAlpha_
    (
        readScalar
        (
            alpha1.mesh().solverDict(alpha1.name()).lookup("cAlpha")
        )
    ),

    sigmaPtr_(surfaceTensionModel::New(dict, alpha1.mesh())),

    deltaN_
    (
        "deltaN",
        1e-8/pow(average(alpha1.mesh().V()), 1.0/3.0)
    ),

    alpha1_(alpha1),
    U_(U),

    nHatf_
    (
        IOobject
        (
            "nHatf",
            alpha1_.time().timeName(),
            alpha1_.mesh()
        ),
        alpha1_.mesh(),
        dimensionedScalar("nHatf", dimArea, 0.0)
    ),

    K_
    (
        IOobject
        (
            "interfaceProperties:K",
            alpha1_.time().timeName(),
            alpha1_.mesh()
        ),
        alpha1_.mesh(),
        dimensionedScalar("K", dimless/dimLength, 0.0)
    ),
    // 2025/5/29 add start
    Kfinal_
    (
        IOobject
        (
            "interfaceProperties:Kfinal",
            alpha1_.time().timeName(),
            alpha1_.mesh()
        ),
        alpha1_.mesh(),
        dimensionedScalar("Kfinal", dimless/dimLength, 0.0)
    )
    // 2025/5/29 add end
{
    // 2025/5/19 add start
    Info << endl;
    Info << "interface properties parameters" << endl;
    if (SSFModel_)
    {
        Info << "Selecting SSF Surface Tension Model" << endl;
        Info << "cSmooth = " << cSmooth_ << endl;
        Info << "cSharp = " << cSharp_ << endl;
    }
    else if (nAlphaSmoothCurvature_ < 1)
    {
        Info << "Selecting CSF Surface Tension Model" << endl;
    }
    else
    {
        Info << "Selecting Smoothed CSF Surface Tension Model" << endl;
        Info << "nAlphaSmoothCurvature = " << nAlphaSmoothCurvature_ << endl;
    }
    Info << endl;
    // 2025/5/19 add end
    calculateK();
}


// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //

Foam::tmp<Foam::volScalarField>
Foam::interfaceProperties::sigmaK() const
{
    return sigmaPtr_->sigma()*K_;
}

// 2025/5/29 add start
Foam::tmp<Foam::surfaceScalarField>
Foam::interfaceProperties::sigmaKfinal() const
{
    return fvc::interpolate(sigmaPtr_->sigma())*Kfinal_;
}
// 2025/5/29 add end

Foam::tmp<Foam::surfaceScalarField>
Foam::interfaceProperties::surfaceTensionForce() const
{
    // 2025/5/19 mod start
    if (SSFModel_)
    {
        const volScalarField alphaSh
        (
            "alphaSh",
            (
                min(
                    max(alpha1_, cSharp_ / scalar(2)),
                    scalar(1) - cSharp_ / scalar(2)
                )
                - cSharp_ / scalar(2)
            ) / (scalar(1) - cSharp_)
        );
    // 2025/5/29 mod start
        return sigmaKfinal()*fvc::snGrad(alphaSh);
    // 2025/5/29 mod end
    }
    else
    {
        return fvc::interpolate(sigmaK())*fvc::snGrad(alpha1_);
    }
    // 2025/5/19 mod start
}


Foam::tmp<Foam::volScalarField>
Foam::interfaceProperties::nearInterface() const
{
    return pos0(alpha1_ - 0.01)*pos0(0.99 - alpha1_);
}


void Foam::interfaceProperties::correct()
{
    calculateK();
}


bool Foam::interfaceProperties::read()
{
    alpha1_.mesh().solverDict(alpha1_.name()).lookup("cAlpha") >> cAlpha_;
    sigmaPtr_->readDict(transportPropertiesDict_);

    return true;
}


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