準(zhǔn)備工作

1. 拷貝求解器源文件：cp -r $FOAM_APP/solvers/multiphase/twoPhaseEulerFoam $WM_PROJECT_USER_DIR/twoPhaseEulerFoam-1

2. 拷貝算例：cp -r $FOAM_TUTORIALS/multiphase/twoPhaseEulerFoam/laminar/fluidisedBed $FOAM_RUN/fluidisedBed

3. 生成算例網(wǎng)格：blockMesh -case $FOAM_RUN/fluidisedBed

4. 修改求解器make文件$WM_PROJECT_USER_DIR/twoPhaseEulerFoam/Make/file:

EXE = ($FOAM_RUN)/twoPhaseEulerFoam

5. 運(yùn)行$WM_PROJECT_USER_DIR/twoPhaseEulerFoam/Allwmake，編譯求解器碉输。

新建曳力模型

1. 進(jìn)入目錄cd $WM_PROJECT_USER_DIR/twoPhaseEulerFoam/interfacialModels/dragModels籽前，復(fù)制其中一個(gè)曳力模型文件夾并命名為EMMScp -r WenYu/ EMMS。然后進(jìn)入到EMMS文件夾下敷钾，把文件名和文件內(nèi)部涉及到模型名稱(chēng)的地方全部替換成EMMS枝哄。

2. 打開(kāi)vim $WM_PROJECT_USER_DIR/twoPhaseEulerFoam/interfacialModels/Make/files，向其中添加一行dragModels/EMMS/EMMS.C阻荒，并將生成目標(biāo)進(jìn)行修改：

//- 修改前
LIB = $(FOAM_LIBBIN)/libcompressibleEulerianInterfacialModels
//- 修改后
LIB = $(FOAM_USER_LIBBIN)/libcompressibleEulerianInterfacialModels_hemm

3. 在目錄$WM_PROJECT_USER_DIR/twoPhaseEulerFoam/interfacialModels下運(yùn)行編譯生成鏈接庫(kù)挠锥，運(yùn)行wmake。生成的鏈接庫(kù)可以在目錄$FOAM_USER_LIBBIN下：

[mmhe@K226 interfacialModels]$ ll $FOAM_USER_LIBBIN
總用量 1.9M
drwxr-xr-x 2 mmhe ipe   61 10月  9 2019 .
drwxr-xr-x 3 mmhe ipe   24 9月  27 15:30 ..
-rwxr-xr-x 1 mmhe ipe 1.9M 10月  9 2019 libcompressibleEulerianInterfacialModels_hemm.so

4. 修改vim $WM_PROJECT_USER_DIR/twoPhaseEulerFoam/Make/options侨赡，在其中將默認(rèn)的鏈接庫(kù)位置換成剛才生成的鏈接庫(kù)：

//- 修改前
-lcompressibleEulerianInterfacialModels
//- 修改后
-L$(FOAM_USER_LIBBIN) \
-lcompressibleEulerianInterfacialModels_hemm

5. 運(yùn)行$WM_PROJECT_USER_DIR/twoPhaseEulerFoam/Allwmake蓖租，編譯求解器粱侣。修改了求解器的Make/options文件后，必須要重新編譯蓖宦，否則無(wú)法更新求解器齐婴。

曳力模型算例測(cè)試

1. 修改vim $FOAM_RUN/fluidisedBed/constant/phaseProperties，將曳力模型設(shè)置成EMMS：

drag
(
    (particles in air)
    {
        type            EMMS;
        residualAlpha   1e-6;
        residualRe      1e-3;
        swarmCorrection
        {
            type        none;
        }
    }
);

2. 修改計(jì)算時(shí)間vim $FOAM_RUN/fluidisedBed/system/controlDict球昨，使算例向前運(yùn)行一步即可：

startTime       0;
stopAt          endTime;
endTime         2e-4;
deltaT          2e-4;

3. 運(yùn)行算例：$FOAM_RUN/twoPhaseEulerFoam -case $FOAM_RUN/fluidisedBed尔店，正常運(yùn)行：

Calculating face flux field phi.air
Selecting diameterModel for phase air: constant
Selecting turbulence model type laminar
Selecting default blending method: none
Selecting dragModel for (particles in air): EMMS
Selecting swarmCorrection for (particles in air): none
Selecting virtualMassModel for (particles in air): constantCoefficient
Selecting heatTransferModel for (particles in air): RanzMarshall
Calculating field DDtU1 and DDtU2

新建固相應(yīng)力模型

1. 拷貝一個(gè)固相應(yīng)力模型：cp -r $WM_PROJECT_USER_DIR/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/phasePressureModel $WM_PROJECT_USER_DIR/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/EMMSStressModel。

2. 將所有模型關(guān)鍵字phasePressure替換為EMMSStress主慰。

3. 修改vim $WM_PROJECT_USER_DIR/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/phaseCompressibleTurbulenceModels.C嚣州，在最后加上EMMSStressModel模型的相關(guān)信息：

#include "EMMSStressModel.H"
#include "EMMSStressModel.H"
makeTurbulenceModel
(phaseModelPhaseCompressibleTurbulenceModel, RAS, EMMSStressModel);

4. 在$WM_PROJECT_USER_DIR/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/Make/files文件中添加一行EMMSStressModel/EMMSStressModel.C，并將目標(biāo)鏈接庫(kù)生成配置改為：

//- 修改前
LIB = $(FOAM_LIBBIN)/libphaseCompressibleTurbulenceModels
//- 修改后
LIB = $(FOAM_USER_LIBBIN)/libphaseCompressibleTurbulenceModels_hemm

5. 在目錄$WM_PROJECT_USER_DIR/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels下運(yùn)行wmake編譯鏈接庫(kù)共螺。

6. 修改vim $WM_PROJECT_USER_DIR/twoPhaseEulerFoam/Make/options该肴，在其中將默認(rèn)的鏈接庫(kù)位置換成剛才生成的鏈接庫(kù)：

-lphaseCompressibleTurbulenceModels_hemm

7. 運(yùn)行$WM_PROJECT_USER_DIR/twoPhaseEulerFoam/Allwmake，編譯求解器藐不，成功匀哄！

應(yīng)力模型算例測(cè)試

1. 修改算例配置文件vim $FOAM_RUN/fluidisedBed/constant/turbulenceProperties.particles，將RAS模型設(shè)置成EMMSStress雏蛮，然后將phasePressure模型的輸入?yún)?shù)拷貝成EMMSStress模型的輸入?yún)?shù)：

RAS
{
    RASModel EMMSStress;
    EMMSStressCoeffs
    {
        preAlphaExp     500;
        expMax          1000;
        alphaMax        0.62;
        g0              1000;
    }
}

2. 運(yùn)行算例：$FOAM_RUN/twoPhaseEulerFoam -case $FOAM_RUN/fluidisedBed涎嚼，正常運(yùn)行：

Calculating face flux field phi.particles
Selecting diameterModel for phase particles: constant
Selecting turbulence model type RAS
Selecting RAS turbulence model EMMSStress
EMMSStressCoeffs
{
    preAlphaExp     500;
    expMax          1000;
    alphaMax        0.62;
    g0              1000;
}
Calculating face flux field phi.air
Selecting diameterModel for phase air: constant
Selecting turbulence model type laminar
Selecting default blending method: none
Selecting dragModel for (particles in air): EMMS
Selecting swarmCorrection for (particles in air): none
Selecting virtualMassModel for (particles in air): constantCoefficient
Selecting heatTransferModel for (particles in air): RanzMarshall
Calculating field DDtU1 and DDtU2

修改曳力模型

1. 打開(kāi)文件vim $WM_PROJECT_USER_DIR/twoPhaseEulerFoam/interfacialModels/dragModels/EMMS/EMMS.C，只需要修改其中的函數(shù)Foam::tmp<Foam::volScalarField> Foam::dragModels::EMMS::CdRe() const即可挑秉。
   CdRe()函數(shù)返回有效曳力系數(shù)與雷諾數(shù)的乘積法梯，其中的變量CdsRe表示標(biāo)準(zhǔn)曳力系數(shù)與雷諾數(shù)的乘積。

Foam::tmp<Foam::volScalarField> Foam::dragModels::EMMS::CdRe() const
{
    Info <<"EMMS曳力模型的CdRe函數(shù)Ｏ拧Ａ⒀啤！"<<endl;
    volScalarField alpha2(max(scalar(1) - pair_.dispersed(), residualAlpha_));
    volScalarField Re(pair_.Re());
    volScalarField CdsRe
    (
        neg(Re - 1000)*24.0*(1.0 + 0.15*pow(Re, 0.687))
      + pos(Re - 1000)*0.44*max(Re, residualRe_)
    );
    //- 獲取EMMS模型輸入?yún)?shù)
    const phaseModel& dispersed = pair_.dispersed();//離散相
    const phaseModel& continuous = pair_.continuous();//連續(xù)相
    const volVectorField Up = dispersed.U();//顆粒速度
    const volVectorField Ug = continuous.U();//流體速度
    const volScalarField Ur = Up.component(vector::Z)-Ug.component(vector::Z);//相對(duì)速度
    const volScalarField epsl_g = pair_.continuous();//平均空隙率

    //- 非均勻結(jié)構(gòu)因子
    volScalarField Hd = Ur*epsl_g;//擬合關(guān)聯(lián)式
    Hd.dimensions().reset(dimensionSet(0,0,0,0,0));//設(shè)置單位為無(wú)量綱數(shù)
    

    return
        Hd
       *CdsRe
       *pow(alpha2, -2.65)
       *max(pair_.continuous(), residualAlpha_);
}

2. 在目錄$WM_PROJECT_USER_DIR/twoPhaseEulerFoam/interfacialModels下運(yùn)行編譯生成鏈接庫(kù)姻灶，運(yùn)行wmake铛绰。

3. 運(yùn)行算例，從輸出信息中可以確定已經(jīng)能夠執(zhí)行該函數(shù)产喉，但是因?yàn)椴皇諗慷崆敖K止捂掰。

修改應(yīng)力模型

devRhoReff()函數(shù)：

divDevRhoReff()函數(shù)：

在phasePressure模型中，是不考慮固相應(yīng)力的曾沈，因此這兩個(gè)函數(shù)返回的都是0值尘颓。

第一步要改的是四個(gè)函數(shù)。首先是pPrime函數(shù)：

1. 修改vim $WM_PROJECT_USER_DIR/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/EMMSStressModel/EMMSStressModel.H：

//- 刪除私有變量-Kenetic Theory Model coefficients
/*
        // Kinetic Theory Model coefficients

            //- Maximum packing phase-fraction
            scalar alphaMax_;

            //- Pre-exponential factor
            scalar preAlphaExp_;

            //- Maximum limit of the exponential
            scalar expMax_;

            //- g0
            dimensionedScalar g0_;
*/

2. 修改vim $WM_PROJECT_USER_DIR/twoPhaseEulerFoam/phaseCompressibleTurbulenceModels/EMMSStressModel/EMMSStressModel.C：

//- 構(gòu)造函數(shù)-刪除沒(méi)有定義的變量初始化
//- read()函數(shù)-刪除文件數(shù)據(jù)讀取
//- pPrime()和pPrimef()函數(shù)-替換掉返回表達(dá)式為擬合關(guān)聯(lián)式
//- devRhoReff()函數(shù)-復(fù)制kineticTheoryModels中對(duì)應(yīng)的部分到此處
//- divDevRhoReff()函數(shù)-復(fù)制kineticTheoryModels中對(duì)應(yīng)的部分到此處

3. 編譯鏈接庫(kù)晦譬，成功！運(yùn)行算例互广，報(bào)錯(cuò)敛腌！

PIMPLE: iteration 1
new cannot satisfy memory request.
This does not necessarily mean you have run out of virtual memory.
It could be due to a stack violation caused by e.g. bad use of pointers or an out of date shared library

4. 經(jīng)過(guò)檢查卧土，是我定義Ur的時(shí)候觸動(dòng)的錯(cuò)誤。問(wèn)題解決像樊，應(yīng)該是我使用了返回的引用尤莺，因此導(dǎo)致了運(yùn)行時(shí)發(fā)生段錯(cuò)誤。