How to calculate the frictional and form drags in Navier_solver?

[Watchmojo]

Hi,

I am interested to add two functions in Navier_solver that calculate the skin friction drag and form drag over the walls, which require two integrations as follows:

Can I get some hint?

Thank you very much.

[jandrej]

Both can be done in the application code (meaning no change to Navier is necessary). I would do the following for skin friction

1. Retrieve the velocity GridFunction. (You should have that anyways because of initial conditions etc.)
2. Create a GridFunctionCoefficient from it
3. Create a new LinearForm
4. Add the BoundaryNormalLFIntegrator with the GridFunctionCoefficient you just created
5. Now you can (at any time you like, because the coefficients are automatically updated!) assemble this form and retrieve the quantity.

See https://github.com/mfem/mfem/blob/master/miniapps/navier/navier_solver.cpp#L197 for similar usage and a hint how to use attributes for your "wall" that you want to integrate on.

This can be done more straightforward, see the QoI computation in the TGV example. I recommend using the above mentioned way first.

[Watchmojo]

Thank you for your help.

Here is what I have done for the form drag

`void NavierSolver::ComputeCp()
{
if (Cp_gfcoeff == nullptr)
{
Cp_gfcoeff = new VectorGridFunctionCoefficient(&pn_gf);
Cp_form = new ParLinearForm(pfes);
Cp_form->AddBoundaryIntegrator(new BoundaryNormalLFIntegrator(
*Cp_gfcoeff), vel_ess_attr);
}
Cp_form->Assemble();
Cp_form->ParallelAssemble(Cp_bdr);

}`

If I understand correctly, at the end of this function, Cp_bdr contains the integration of pn_gf in the normal direction (the interior grids/modes are marked zero). If the mesh is perfectly orthogonal and there is only one wall normal direction, summing Cp_bdr over all the processors will give the form drag. Right? If yes, what if the mesh is not orthogonal? In that case, Cp_bdr should be integrated while considering arbitrary wall normal directions for different locations. That is, the form drag in x, y, z directions should be calculated separated. How can this be done?

---

I have done the following for the skin friction:
`void NavierSolver::ComputeCf()
{
if (Cf_gfcoeff == nullptr)
{
Cf_gfcoeff = new VectorGridFunctionCoefficient(&Cf_gf);
Cf_form = new ParLinearForm(vfes);
Cf_form->AddBoundaryIntegrator(new BoundaryNormalLFIntegrator(
*Cf_gfcoeff), vel_ess_attr);
}
G->Mult(un, Cfn); // this has to be wrong
Cf_gf.GetTrueDofs(Cfn);
Cf_form->Assemble();
Cf_form->ParallelAssemble(Cf_bdr);

}
`
The way of calculating velocity gradient above has to be wrong. Another function returning the velocity gradient tensor probably needs to be added, in analogy to ComputeCurl2D/ComputeCurl3D. The remaining problem is then: how to integrate Cf_bdr while considering arbitrary wall normal directions for different locations?

---

By the way, in ComputeCFL, there might be two minor issues: (1) hmin is the minimal dimension of element. Should it be the minimal grid size for each integration point? (2) Should CFL be calculated as follows:
` cflmax = fmax(cflmax, cflx);
cflmax = fmax(cflmax, cfly);
cflmax = fmax(cflmax, cflz);

`----
A rudimental question: in the case of Gauss-Lobatto bases, the nodal value is equivalent to the modal coefficient. Right?

[jandrej]

• I'm not sure what you mean by your statement that you need the different directional components. Your mentioned integral is fulfilled on any boundary.

• You are right regarding the gradient of the velocity. You can either implement something similar to https://github.com/mfem/mfem/blob/master/fem/coefficient.hpp#L519 for vectors, or go the custom function/integrator route like in ComputeCurl2D or the QoI computation in the TGV example.

• The element length is a choice in the CFL computation. The choice here is consistent with the one in Nek5000 for example.

• I'm not familiar with modal methods.

[Watchmojo]

hi Jandrej

I meant that the different directional components over walls should be integrated globally to get the overall drag. Therefore, the variation of surface norm needs to be accounted for, for example, the pressure drag on the circular cylinder.

In miniapp Navier_solver, gradient and curl operators are implemented. Why did not you use the functions provided in GridFunction?

In the first substep of time splitting, it appears that the nonlinear convective term is projected to the test function space and then the nodal velocity field is updated by dividing the inverse of mass matrix. Why not calculate the nonlinear convective term at each node, and then the nodal velocity field can be advanced without using the mass matrix?

I have tried to implement NS solver using MFEM 5 years ago, although it failed. I will read more of your implementation and absorb for future research.
NS_.zip