The flux correction at the 1st face

[JihooKang-KOR]

Dear @AndreWeiner,

In addition to the 1st issue, I would like to mention about the correction at the 1st face.

In the 1st issue, I showed that the 2D mapping model didn't work well because the slopes at the wall were too large to converge. Meanwhile, the slopes at the final step (3539th) for the finest mesh case were the similar to the range of the slopes for the 1D channel flow case. Therefore, I kept the original 1D mapping model for the simulation, and came up with another approach.

Regarding the notebook "MeshEffect_Spalding.ipynb" (the commit 997d829), the graph in the notebook can be applied only when the velocity profile of the real simulation is similar to the Spalding's function profile. However, when we turn on the face correction, the velocity profile drastically changes, and thus it isn't similar to the Spalding's function anymore. In contrast, the wall correction works well (steadily increasing wall shear stresses), and smoothly changes the velocity field and the velocity gradient at the wall.

Hence, I decided that I fix "magGradUf" (velocity gradient at the 1st face) at a certain time step as the commit 8d809a7. In the code, I fixed the magGradUf at 6th time step with the variable name of "fixedGradUf". I was not sure if the velocity profile at this time step was similar to the Spalding's function, but I thought that this was the time step before the drastic change in the velocity gradient at the face. Therefore, I tried this setting for y+ = 30, 50, and 100, and the result of skin friction is as follows:

The result shows that the simulation at least converged. If we let the face correction happen with changing magGradUf every time step, the magGradUf diminishes very fast, and then the correction ratio of nuEff always exceeds 1. This makes the simulation crashes because magGradUf reduces further and further.

However, when the face correction with the fixed denominator (fixedGradUf) at a certain early time step turns on, the nuEff correction ratio is less than 1 which can lead to convergence. Therefore, my speculation is that we can properly correct fluxes at the 1st face as well if we are able to find which time step has the velocity profile that is similar to the Spalding's function.

I need to think about any methods how to find it if my speculation is reasonable.

Best regards, Jihoo Kang

[AndreWeiner]

Dear Jihoo, very interesting observation. I think what it means is that the problem with the correction at the first cell face is rather a stability issue than a conceptual problem.

One thing I noticed when browsing through the code again: the code below should not really be needed, since later on, you add SMALL when dividing by the gradient's magnitude. Therefore, I would remove it.

forAll (oppFaceIDs, faceI)
{
    // Set an arbitrary value in order to avoid zero denominator at the variable tnuf and 
    // overflow of the exponential term in the function spaldings_law for the first time step
    if (mag(U_sngrad.internalField()[oppFaceIDs[faceI]]) == 0)
    {
        tmp_magGradUf[faceI] = 1e5;
    }        
    else
    {
        tmp_magGradUf[faceI] = mag(U_sngrad.internalField()[oppFaceIDs[faceI]]);
    }        
}

Regarding the stability issue, I suggest the following analysis:

• perform the plate simulation without correcting the first cell face viscosity for yPlus={50, 100}
• for probe location at the beginning, middle, and end of the plate, plot the main quantities of interest at the first cell face over time (the correction factor, the model slope, the numerical slope, the effective viscosity)
• my hope is that we can devise some rules to stabilize the face correction based on these curves

Best, Andre

[AndreWeiner]

Hi Jihoo, I am currently checking the last commit. I already observed the following:

• for both yPlus values, the model slopes at the beginning and the middle of the plate are essentially zero; at the probe location in the middle, the model slopes take on non-zero values after the first couple of iterations; at the beginning of the plate, the model slope remains essentially zero; if we then apply the correction nueff *= model_slope/numerical_slope, the viscosity is set to zero; zero viscosity certainly causes stability issues and might explain the behavior you're observing; so, we have to think about a meaningful lower limit for the effective viscosity; one option might be nueff=max(nueff, nu); another option might be to consider also the cell-centered nueff of the second cell normal to the wall

The second idea is slightly more involved but also more promising I believe. Here is some pseudo-code of how it might be implemented:

// nueff_1 - cell-centered effective viscosity at first cell-center normal to the wall
// nueff_2 - cell-centered effective viscosity at second cell-center normal to the wall
nueff_face_12 *= model_slope/numerical_slope;
nueff_max = max(nueff_1, nueff_2);
nueff_min = min(nueff_1, nueeff_2);
nueff_face_12 = min(max(nueff_min, nueff_face_12), nueff_max);

Here is a code snippet allowing you to find the second cell-center indices (source):

// Get IDs of secondary adjacent cells
labelList secAdjacentCellIDs(adjacentCellIDs.size());

label globFace = -1;
forAll (oppFaceIDs, faceI)
{
    globFace = oppFaceIDs[faceI];

    if (mesh.owner()[globFace] == adjacentCellIDs[faceI])
    {
        secAdjacentCellIDs[faceI] = mesh.neighbour()[globFace];
    }
    else
    {
        secAdjacentCellIDs[faceI] = mesh.owner()[globFace];
    }
}

Let's talk about this idea in more detail. Good work with the analysis! Best, Andre

[JihooKang-KOR]

Dear @AndreWeiner,

I limited the velocity as per the discussion in the last meeting, and the code is uploaded as the commit 3410693.

First of all, the important part of the code is as follows:

forAll (oppFaceIDs, faceI)
{
    // 1st face correction only if the differences at the cell center and the face are negative
    if (diffU_center[faceI] <= 0 && diffU_face[faceI] <= 0)
    {
        nuEff[oppFaceIDs[faceI]] *= (labelAccessor[faceI][1])/(magGradUf[faceI] + ROOTVSMALL);
        actFaceCorr[faceI] = 1;
    }

    // 1st face correction also if the differences at the cell center and the face are positive
    else if (diffU_center[faceI] > 0 && diffU_face[faceI] > 0)
    {
        nuEff[oppFaceIDs[faceI]] *= (labelAccessor[faceI][1])/(magGradUf[faceI] + ROOTVSMALL);
        actFaceCorr[faceI] = 1;
    }

    else
    {
        actFaceCorr[faceI] = 0;
    }
}

We should consider the time derivative here, not the spatial derivative. Therefore, I calculated the difference of the velocities between two nearest time steps. (fvc::ddt didn't work because the case is a steady-state simulation.) When the sign of the difference at the cell center and at the face is the same, the 1st face correction executes, and the case is named "turbulentFlatPlate_velFixPos_500iter" (as the above code). On the other hand, the case is called "turbulentFlatPlate_velFixNeg_500iter" when the correction only executes for the negative sign.

forAll (oppFaceIDs, faceI)
{
    // 1st face correction only if the differences at the cell center and the face are negative
    if (diffU_center[faceI] <= 0 && diffU_face[faceI] <= 0)
    {
        nuEff[oppFaceIDs[faceI]] *= (labelAccessor[faceI][1])/(magGradUf[faceI] + ROOTVSMALL);
        actFaceCorr[faceI] = 1;
    }

    else
    {
        actFaceCorr[faceI] = 0;
    }
}

The skin friction graphs for y+ = 50 and 100 are as follows: (The correction only for the negative sign)

(The correction for the positive and the negative sign)

It seems to be okay, but there is a residual issue as follows: (The correction only for the negative sign)

(The correction for the positive and the negative sign)

They don't converge at all, but the correction only for negative sign shows a bit better situation.

Here, I would like to show the velocity graphs which are given by: (The correction only for the negative sign)

(The correction for the positive and the negative sign)

We can see that the correction only for negative sign gives us the better result. However, as I mentioned above, the residual doesn't converge. Maybe I would execute a full simulation using the negative sign correction.

Best regards, Jihoo Kang

[AndreWeiner]

Dear Jihoo,

We should consider the time derivative here, not the spatial derivative.

That must have been a misunderstanding. The idea is to preserve the monotonicity of the velocity, so spatial derivatives are required. In pseudocode, the correction should look as follows:

// x - direction along the plate
// Ux_wall - velocity at the wall
// Ux_1 - velocity in the first cell center normal to the wall
// Ux_2 - velocity in the second cell center normal to the wall
if ( (Ux_wall - Ux_1)*(Ux_1 - Ux_2) > 0) {
   do correction;
} else {
  do something else;
}

Best, Andre

[JihooKang-KOR]

Dear @AndreWeiner,

I revised the code again as per your comment (the commit 8879835). The code uploaded in Github uses the velocity at the 2nd cell center, whereas the another version uses the velocity at the 1st cell face in the if condition (made and checked internally).

The skin friction results for y+ = 50 and 100 are as follows: (The criterion of 2nd cell center velocity)

(The criterion of 1st cell face velocity)

I simulated with 500 iteration steps, and the simulation didn't diverge, but the skin friction prediction was not that precise. I guess this limitation method might prevent that the velocities at the 1st cell center and at the 1st cell face are exactly the same. However, this method doesn't prevent that these velocities get closer each other. Therefore, magGradUf is anyway very small that leads to higher skin friction values, and the velocities meet together and look almost the same after some iteration steps as follows:

In addition, there is a residual issue again.

All the plots from the Ux plots are the case of 2nd cell center velocity because another case has almost the similar graphs.

Best regards, Jihoo Kang

[JihooKang-KOR]

Dear @AndreWeiner,

I executed the simulation with the code (f7b94f9) of fixing magGradUf and updating it every 50 iterations for y+ = 50 and 100.

The skin friction graph is as follows:

The residual graphs are given by:

For 500 iteration steps, the case for y+ = 100 converged, but the case for y+ = 50 didn't converge. Therefore, I plot several fields for y+ = 100 as follows:

According to these graphs, magGradUf didn't converge as time increased. At the middle and the end of the plate, the value continuously decreases. Even though the simulation for y+ = 100 converged, I cannot think that this simulation yields any meaningful results, which showed that the method just delays the increase of nuEff at the first face.

At this point, I thought there would be any interruption between two corrections at the wall and the face. Thus, I checked the velocity field without any correction that is given by: (No correction) (Wall correction, for comparison)

Compared to the wall correction velocity field, the value is much larger (not approximately 40, but 60 at the end of the plate). Due to the wall correction, the velocity at the 1st cell center a lot decreases, and subsequently the ML model predicts the face slopes much higher than it should be (and the graph shape will be far from the Spalding's function). That is why the model slopes were always larger than the actual magGradUf, I guess.

In order to clarify if the above speculation is reasonable, I tested two cases.

1. Only turn on the 1st face correction (no wall correction), and multiplies nuWall = np.linspace(7.44, 6.25, 449) to add the wall correction ratio at the post-processing stage. (Skin friction Cf = nuEff_wall * magGradUf)

2. Turn on both methods (for face correction, from 31 time steps), but input velocities should be added 25.0 m/s. For example, if the actual velocity at the cell 200 is 40, the model predicts the slope based on 65.

Both skin friction graphs (small-dotted red lines in both graphs) show the reasonable values, even though they are quite flat shape. Therefore, I would like to conclude that the both correction methods are incompatible each other, and thus we need any mitigation methods to apply them together.

-------------------- Additional comment ------------------------- When I tried to apply the above to the Spalding's wall function model, I realized that it cannot be applied to the model because if the velocity at the same location is faster than before, the related slope is always larger. The shape of the Spalding's function is assumed to be always fully developed. Therefore, the above situation can only be applied to the ML model which has both of developing and fully developed state. I think we need to check the slopes in the 1D channel flow case again. If the time step for the flat plate case increases until a fully developed state but the model predicts the slope from the developing state in the 1D case, then it would also not be useful.

Best regards, Jihoo Kang

[AndreWeiner]

Hi Jihoo, thanks for the update. For the tests with both corrections activated, do the velocity and pressure fields look at least physical? If not, can you figure out, when approximately they become unphysical? (visual inspection with paraview) Best, Andre

[AndreWeiner]

I think the main issue is that the velocity value in the second cell normal to the wall takes on values close to the ones in the first cell, which should not be the case. Therefore, we have to understand why the velocity in the second cell decreases. But first, it is important to know if the overall flow field looks physical (we know what it should look like), e.g.,

• does the velocity field show or over or undershoots (values less than zero or larger than the inlet velocity)?
• does the pressure field show any weird patterns like checkerboarding?

If the flow field becomes unphysical early on, then the low velocity in the second cell is likely only a consequence thereof, and we have to understand why the solution is unphysical.

Best, Andre

[JihooKang-KOR]

Dear @AndreWeiner,

Regarding velocity fields, there is a fluctuation at the front of the plate as follows (inlet 69.4m/s): \<Wall correction>

\<Wall & first face correction>

I tested the case of no correction, but it was wrong because the case didn't use any wall function, and thus it couldn't capture the real behavior. Therefore, I turned on the wall function, and executed the simulation without any correction at wall and faces. (And automatically, the velocity plus 25.0 case is not valid anymore.)

As per the figure, there is also a kink at the front of the plate. The whole velocity field distribution is little bit different from the other cases, but we can conclude that the field shape is not that different each other.

With regard to pressure fields, for the case of 1st face correction, there are some strange patterns in the field as follows: But I think they appear due to the consequence of the same velocity at two locations.

I checked the fields in the case of wall function activated as follows:

As seen in the figures, we can find that the slopes with the wall function from a certain time step are the similar to the wall correction case. This means that the idea of wall correction is okay. However, when the 1st face correction is involved, all the problems occur given by:

1. The model slopes at the 1st face are larger than the numerical slopes after a certain time step -> nuEff will be continuously larger
2. These model slopes only depend on the 1st cell center velocity, and therefore when the numerical slopes decrease, the model slopes remain almost the same, then nuEff will always be much larger and larger.

According to the above, I would like to say that we're back to the problem that we faced last year. I think the wall function works because it doesn't touch anything at the 1st face. Hence, I would like to suggest that we can extract any ratio information from the first face and multiply to the wall correction factor, not to the face correction factor.

Best regards, Jihoo Kang

[AndreWeiner]

Hi Jihoo, thanks for the update. The velocity fluctuation at the beginning of the plate in the following picture is exactly what causes the issues further downstream (towards the middle/end of the plate):

The good news is that I believe to know how to fix this issue (for real this time :-)). We've had similar issues with the original implementation for mass transfer problems, but there, the behavior was caused by the correction of the convection term. The velocity boundary layer is different, though, because the momentum equation is non-linear and, instead, species is just transported passively. I checked the original implementation again and noticed the following:

1. it is crucial to ensure monotonicity; what we did in the previous implementation is slightly different from what we tried here, and I think that's where the issue is; have a look at the following code snippet:

   scalar mono = (Ci[adjICellIDs[faceI]] - Ci[secAdjICellIDs[faceI]])
                        *(Ci[secAdjICellIDs[faceI]] - Ci[thirdAdjICellIDs[faceI]]);
   
           if ( (CiNUM[oppIFaceIDs[faceI]] > SMALL) && (mono > 0) )
           {
               phiSgs_i[oppIFaceIDs[faceI]] =  phi[oppIFaceIDs[faceI]]
                                             *CiSGS[oppIFaceIDs[faceI]]
                                             /CiNUM[oppIFaceIDs[faceI]];
           }

   As you can see, monotonicity should be checked by comparing the sign of the slope between the first and second cell center with the one between the second and third cell center. The correction at the first face should be only applied if mono > 0. With this criterion, in a configuration as in the picture above, no correction would be applied because the velocity goes up and then down again within the first three cells. Finding the third cell center is analogous to finding the second cell center. Here is a snippet:

   
   // Get IDs of secondary adjacent cells
   labelList secAdjICellIDs(adjICellIDs.size());

label globFaceI = -1; forAll (oppIFaceIDs, faceI) { globFaceI = oppIFaceIDs[faceI];

if (mesh.owner()[globFaceI] == adjICellIDs[faceI])
{
    secAdjICellIDs[faceI] = mesh.neighbour()[globFaceI];
}
else
{
    secAdjICellIDs[faceI] = mesh.owner()[globFaceI];
}

}

// Get IDs of faces of the cell that are opposite to the opponent face labelList secOppIFaceIDs(oppIFaceIDs.size());

forAll (oppIFaceIDs, faceI) { secOppIFaceIDs[faceI] = mesh.cells()[secAdjICellIDs[faceI]].opposingFaceLabel ( oppIFaceIDs[faceI],mesh.faces() ); }

// Get IDs of third adjacent cells labelList thirdAdjICellIDs(adjICellIDs.size());

label globFaceThI = -1; forAll (secOppIFaceIDs, faceI) { globFaceThI = secOppIFaceIDs[faceI];

if (mesh.owner()[globFaceThI] == secAdjICellIDs[faceI])
{
    thirdAdjICellIDs[faceI] = mesh.neighbour()[globFaceThI];
}
else
{
    thirdAdjICellIDs[faceI] = mesh.owner()[globFaceThI];
}

}

2. If the linear slope at the first cell  becomes too small, don't apply any correction; here is another snippet:

if ( mag(snGradCi[oppIFaceIDs[faceI]]) > SMALL ) { scalar corrFactDopp = mag(snGradCiSGS[oppIFaceIDs[faceI]] /snGradCi[oppIFaceIDs[faceI]]); ... }

I hope this helps! Let me know if I can support you with the implementation.
Best, Andre

[JihooKang-KOR]

Dear @AndreWeiner,

I executed a simulation based on your comment (the related code in the commit 5ad87be).

First of all, the skin friction plot with 500 iteration steps for y+ = 50 and 100 is given by:

As seen in the figure, it at least didn't diverge until 500th step, but the value is not that stable. The residual plot is as follows:

The simulation didn't converge again according to the plots above.

Here are the figures about pressure and velocity fields at 500th step by Paraview: \<Pressure field near the end of the plate> \<Velocity fields at the end of the plate>

There exist several back pressure areas and checker board phenomenon in the velocity field. This implies that the simulation would not be valid.

The field for velocities and slopes are as follows:

As you can see, two velocities meet together again at the middle and the end of the plate. I saw the csv files for this simulation also, and I found that the if statement for monotonicity only worked at the front of the plate because at the middle and the end of the plate, the velocities were mostly monotonic. At those areas, monoton > 0 wasn't activated, but only magGradUf[faceI] > SMALL was activated. It can make the simulation more stable and keep two velocities at the 1st cell center and the face (or 2nd cell center) different, but they are again similar. Therefore, as time increases, nuEff at the face grows over and over again before magGradUf is less than SMALL.

I think the cause of this phenomenon is that the predicted model slopes are always larger than the numerical slopes at the middle and the end of the plate after a certain time, and they are not explicitly connected each other. Furthermore, I guess this velocity field is directly affected by the face correction unlike concentration fields. Then, the predicted slope is calculated in a new equilibrium state which is totally different from what we animated by using Spalding's function.

The current equilibrium state is almost the identical to the simulation with the wall function. Hence, there would be no error regarding the values of magGradUf for the case of wall correction. In a nutshell, I think that the model slopes would be intrinsically larger than the numerical slopes after some steps because the current equilibrium state is totally different from the ideal condition with Spalding's function.

Best regards, Jihoo Kang

[AndreWeiner]

Thanks, Jihoo, good job testing the modification. It's a pity that it didn't resolve all the issues, though. It appears that at least the velocity overshoots at the beginning of the plate are gone. Regarding the species implementation: the flux over the first face is predicted based on the cell-centered concentration value and changing the flux also changes the cell-centered value until an equilibrium is found. Could you send me a snapshot of the current setup including everything that is needed to run the simulation (e.g., the models)? Best, Andre

[JihooKang-KOR]

Dear @AndreWeiner,

I uploaded the whole files related to this case in the cloud folder.

Please find the folder name "turbulentFlatPlate_1stface_monoton". In the folder, there are the pytorch models, csv files for fields, and fields' status for every 100 iteration steps (until the 500th step) which can be visualized by paraview.

If you need more things or I omitted something, please kindly let me know.

Best regards, Jihoo Kang

[JihooKang-KOR]

Dear Andre,

Please find the plots mentioned in the e-mail I sent. For comparison, the plots of no correction for y+ = 1 and 100 are also located here.

The green graphs are the data-driven model with velocity blending with 5000 iterations. For y+ = 10, Cf seemed to be okay for 500 iterations in the notebook, but we can see that it increased after 5000 iterations. This would also be one issue, but when we compare this to the no wall function version, the values are similar. Therefore, the first priority will be at y+ = 5.

Best regards, Jihoo Kang