Calculating lift and drag

[Lukas Bystricky]

I am trying to replicate the results here of the benchmark problem of flow around a cylinder (circle in 2D). Here is my code so far:

template<int dim>
void NavierStokesSolver<dim>::calculate_lift_and_drag()
{
    QGauss<dim-1>               face_quadrature_formula(3);
    const int                   n_q_points = face_quadrature_formula.size();
   
    const FEValuesExtractors::Vector velocities (0);
    const FEValuesExtractors::Scalar pressure (dim);
   
    std::vector<double>     pressure_values(n_q_points);
    std::vector<Tensor< 2, dim> > velocity_gradients(n_q_points);
   
    Tensor<1,dim>    grad_u_tau;     
    Tensor<1,dim>    normal_vector;
    Tensor<1,dim>    tangent_vector(1);
   
    FEFaceValues<dim> fe_face_values (fe, face_quadrature_formula, update_values
            | update_quadrature_points | update_gradients | update_JxW_values | update_normal_vectors );
   
    typename DoFHandler<dim>::active_cell_iterator 
            cell = dof_handler.begin_active(), endc = dof_handler.end();
   
    double drag = 0;
    double lift = 0;
           
    for (; cell!=endc; ++cell)
    {
        for (int face=0; face < GeometryInfo<dim>::faces_per_cell; face++)
        {                       
            if (cell->face(face)->at_boundary())
            {
                fe_face_values.reinit (cell, face);
                std::vector<Point<dim> > q_points = fe_face_values.get_quadrature_points();
               
                if (cell->face(face)->boundary_indicator() == 4)//on the cirlce
                {           
                    for (int q = 0; q < n_q_points; q++)
                    {
                        fe_face_values[velocities].get_function_gradients(solution, velocity_gradients);
                        fe_face_values[pressure].get_function_values(solution, pressure_values);
                        normal_vector = fe_face_values.normal_vector(q);
                       
                        tangent_vector[0] = normal_vector[1];
                        tangent_vector[1] = -normal_vector[0];
                       
                        //calculate gradient of tangential component of velocity
                        grad_u_tau[0] = velocity_gradients[q][0][0]*tangent_vector[0] + velocity_gradients[q][1][0]*tangent_vector[1];
                        grad_u_tau[1] = velocity_gradients[q][0][1]*tangent_vector[0] + velocity_gradients[q][1][1]*tangent_vector[1];                   
                       
                        //dot product of tangential component of velocity with normal vector
                        drag += (input.nu*(grad_u_tau[0]*normal_vector[0] + grad_u_tau[1]*normal_vector[1])*normal_vector[1]
                                 - pressure_values[q]*normal_vector[0])*fe_face_values.JxW(q);
                                 
                        lift -= (input.nu*(grad_u_tau[0]*normal_vector[0] + grad_u_tau[1]*normal_vector[1])*normal_vector[0]
                                 - pressure_values[q]*normal_vector[1])*fe_face_values.JxW(q);
                    }
                }
            }
        }
    }
   
    //calculate pressure drop
    Point<dim> p1, p2;
    p1[0] = 0.15;
    p1[1] = 0.2;
    p2[0] = 0.25;
    p2[1] = 0.2;
    Vector<double> solution_values1(dim+1);
    Vector<double> solution_values2(dim+1);
   
    VectorTools::point_value(dof_handler, solution, p1, solution_values1);
    VectorTools::point_value(dof_handler, solution, p2, solution_values2);
   
    double p_diff = solution_values1(dim) - solution_values2(dim);
   
    fprintf (pFile, "%f,  %f, %f\n", 50*drag, 50*lift, p_diff);
    fflush (pFile);   
}

The pressure drop is correct, so I think my solver is working, however my lift and drag values are off. To evaluate the lift and drag you need to know the normal derivative of the tangential velocity at the boundary. I've been assuming the since I am doing quadrature, that the gradient of the tangential vector is 0 along each edge so the gradient of the tangential velocity is just the gradient of the velocity times the tangential vector. Is this valid? Or have I made other mistakes?

[Thomas Wick]

I did such computations in the past - maybe the error is
really in your assumptions. The nice thing in deal.II is that
you can just work the compact form. Here are pieces of
code that I used ( I am not posting the entire code because
it actually an ALE fluid code, which can be used for FSI -
but if you have no FSI deformation, then you arrive at a
standard fluid code):


  for (; cell!=endc; ++cell)
     {

  for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
     if (cell->face(face)->at_boundary() &&
         cell->face(face)->boundary_indicator()==80)
       {
         fe_face_values.reinit (cell, face);
         fe_face_values.get_function_values (solution, face_solution_values);
         fe_face_values.get_function_grads (solution, face_solution_grads);
              
         for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
           {

....

     drag_lift_value -= 500.0 * (stress_fluid + fluid_pressure) *
           fe_face_values.normal_vector(q_point)* fe_face_values.JxW(q_point);


}

  } // end boundary 80 for fluid
      
     }
  
   std::cout << "Drag: " << drag_lift_value[0] << std::endl;
   std::cout << "Lift: " << drag_lift_value[1] << std::endl;





Best Thomas

++--------------------------------------------++
Dr. Thomas Wick
Research Scientist at RICAM Linz, Austria

Email:  thoma...@ricam.oeaw.ac.at
www:    http://people.ricam.oeaw.ac.at/t.wick/
++--------------------------------------------++
--

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Lukas Bystricky]

Thank you Thomas. I tried as you suggested, and now the drag calculation seems to work, however the lift is still off. Here's my new code:

template<int dim>
void NavierStokesSolver<dim>::

calculate_lift_and_drag(double t)

{
    QGauss<dim-1>               face_quadrature_formula(3);
    const int                   n_q_points = face_quadrature_formula.size();
   
    const FEValuesExtractors::Vector velocities (0);
    const FEValuesExtractors::Scalar pressure (dim);
   
    std::vector<double>     pressure_values(n_q_points);
    std::vector<Tensor< 2, dim> > velocity_gradients(n_q_points);
   
   
    Tensor<1,dim>    grad_u_tau;     
    Tensor<1,dim>    normal_vector;

    Tensor<1,dim>    tangent_vector;
   
    Tensor<2,dim>    fluid_stress;
    Tensor<2,dim>    fluid_pressure;
    Tensor<1,dim>    forces;

   
    FEFaceValues<dim> fe_face_values (fe, face_quadrature_formula, update_values
            | update_quadrature_points | update_gradients | update_JxW_values | update_normal_vectors );
   
    typename DoFHandler<dim>::active_cell_iterator 
            cell = dof_handler.begin_active(), endc = dof_handler.end();
   
    double drag = 0;
    double lift = 0;
           
    for (; cell!=endc; ++cell)
    {
        for (int face=0; face < GeometryInfo<dim>::faces_per_cell; face++)
        {                       
            if (cell->face(face)->at_boundary())
            {
                fe_face_values.reinit (cell, face);
                std::vector<Point<dim> > q_points = fe_face_values.get_quadrature_points();
               
                if (cell->face(face)->boundary_indicator() == 4)

                {   

                    fe_face_values[velocities].get_function_gradients(solution, velocity_gradients);
                    fe_face_values[pressure].get_function_values(solution, pressure_values);
                           
                    for (int q = 0; q < n_q_points; q++)
                    {                       

                        normal_vector = -fe_face_values.normal_vector(q);
                       
                        fluid_pressure[0][0] = pressure_values[q];
                        fluid_pressure[1][1] = pressure_values[q];
                         
                        fluid_stress = input.nu*velocity_gradients[q] - fluid_pressure;
                       
                        forces = fluid_stress*normal_vector*fe_face_values.JxW(q);
                       
                        drag += forces[0];
                        lift += forces[1];                       

                    }
                }
            }
        }
    }
   
    //calculate pressure drop
    Point<dim> p1, p2;
    p1[0] = 0.15;
    p1[1] = 0.2;
    p2[0] = 0.25;
    p2[1] = 0.2;
    Vector<double> solution_values1(dim+1);
    Vector<double> solution_values2(dim+1);
   
    VectorTools::point_value(dof_handler, solution, p1, solution_values1);
    VectorTools::point_value(dof_handler, solution, p2, solution_values2);
   
    double p_diff = solution_values1(dim) - solution_values2(dim);
   

    fprintf (pFile, "%f, %f,  %f, %f\n", t, 500*drag, 500*lift, p_diff);
    fflush (pFile);   
}

Is there some obvious mistake I'm making?

[Thomas Wick]

How far is the lift off?

The lift is very, very sensitive. Maybe you try mesh refinement.

Best Thomas

++--------------------------------------------++
Dr. Thomas Wick
Research Scientist at RICAM Linz, Austria

Email: thoma...@ricam.oeaw.ac.at
www: http://people.ricam.oeaw.ac.at/t.wick/
++--------------------------------------------++
--

On 03/04/2016 05:44 PM, Lukas Bystricky wrote:
> Thank you Thomas. I tried as you suggested, and now the drag
> calculation seems to work, however the lift is still off. Here's my
> new code:
>

> template<int dim>
> void NavierStokesSolver<dim>::calculate_lift_and_drag(double t)
> {

> QGauss<dim-1> face_quadrature_formula(3);
> const int n_q_points =
> face_quadrature_formula.size();
>
> const FEValuesExtractors::Vector velocities (0);
> const FEValuesExtractors::Scalar pressure (dim);
>
> std::vector<double> pressure_values(n_q_points);
> std::vector<Tensor< 2, dim> > velocity_gradients(n_q_points);
>
>
> Tensor<1,dim> grad_u_tau;
> Tensor<1,dim> normal_vector;

> Tensor<1,dim> tangent_vector;
>
> Tensor<2,dim> fluid_stress;
> Tensor<2,dim> fluid_pressure;
> Tensor<1,dim> forces;
>

> FEFaceValues<dim> fe_face_values (fe, face_quadrature_formula,
> update_values
> | update_quadrature_points | update_gradients |
> update_JxW_values | update_normal_vectors );
>
> typename DoFHandler<dim>::active_cell_iterator
> cell = dof_handler.begin_active(), endc = dof_handler.end();
>
> double drag = 0;
> double lift = 0;
>
> for (; cell!=endc; ++cell)
> {
> for (int face=0; face < GeometryInfo<dim>::faces_per_cell; face++)
> {
> if (cell->face(face)->at_boundary())
> {
> fe_face_values.reinit (cell, face);
> std::vector<Point<dim> > q_points =
> fe_face_values.get_quadrature_points();
>
> if (cell->face(face)->boundary_indicator() == 4)

> {

> fe_face_values[velocities].get_function_gradients(solution,
> velocity_gradients);
> fe_face_values[pressure].get_function_values(solution, pressure_values);
>
> for (int q = 0; q < n_q_points; q++)
> {

> normal_vector = -fe_face_values.normal_vector(q);
>
> fluid_pressure[0][0] = pressure_values[q];
> fluid_pressure[1][1] = pressure_values[q];
>
> fluid_stress = input.nu*velocity_gradients[q]
> - fluid_pressure;
>
> forces =
> fluid_stress*normal_vector*fe_face_values.JxW(q);
>
> drag += forces[0];
> lift += forces[1];
> }
> }
> }
> }
> }
>

> //calculate pressure drop
> Point<dim> p1, p2;
> p1[0] = 0.15;
> p1[1] = 0.2;
> p2[0] = 0.25;
> p2[1] = 0.2;
> Vector<double> solution_values1(dim+1);
> Vector<double> solution_values2(dim+1);
>
> VectorTools::point_value(dof_handler, solution, p1, solution_values1);
> VectorTools::point_value(dof_handler, solution, p2, solution_values2);
>
> double p_diff = solution_values1(dim) - solution_values2(dim);
>

> fprintf (pFile, "%f, %f, %f, %f\n", t, 500*drag, 500*lift, p_diff);
> fflush (pFile);
> }
>

[Thomas Wick]

How far is the lift off?

The lift is very, very sensitive. Maybe you try mesh refinement.

Best Thomas

++--------------------------------------------++
Dr. Thomas Wick
Research Scientist at RICAM Linz, Austria

Email:  thoma...@ricam.oeaw.ac.at
www:    http://people.ricam.oeaw.ac.at/t.wick/
++--------------------------------------------++
--

[Lukas Bystricky]

For the steady case, with Re=20, it's not off by much. It should be in the range 0.0104 to 0.0110 and it is 0.0093. When I go to higher Reynolds numbers though, for example the periodic Re=100, it is off by quite a lot (it should vary between around -1 and 1 and instead it varies between around -20 and 20). This does seem to be very sensitive to mesh and timestep sizes though. My timestepping scheme is only first order, could that be the problem?

Thanks again,

Lukas

[Thomas Wick]

Yes. This might be very likely that for higher Reynolds number
the timestepping and time step size is a problem.

Please use a finer time step size
and please also test with another time stepping scheme,
for example A-stable second order Crank-Nicolson.

Best Thomas

++--------------------------------------------++
Dr. Thomas Wick
Research Scientist at RICAM Linz, Austria

Email:  thoma...@ricam.oeaw.ac.at
www:    http://people.ricam.oeaw.ac.at/t.wick/
++--------------------------------------------++
--

On 03/04/2016 06:00 PM, Lukas Bystricky wrote:

For the steady case, with Re=20, it's not off by much. It should be in the range 0.0104 to 0.0110 and it is 0.0093. When I go to higher Reynolds numbers though, for example the periodic Re=100, it is off by quite a lot (it should vary between around -1 and 1 and instead it's around 20). This does seem to be very sensitive to mesh and timestep sizes though. My timestepping scheme is only first order, could that be the problem?

Thanks again,

Lukas

[Jie Cheng]

Hi Lucas

   

I am trying to replicate the results here of the benchmark problem of flow around a cylinder (circle in 2D). Here is my code so far:

  It's been a long time since you posted this question. But if you are still reading the threads, could you please tell me how you managed to generate a non-symmetric triangulation in the benchmark case? I tried to use one hypercube_with_cylindrical_hole and two hyper_rectangles to construct it, but the result looks ugly upon global refinement.

   Thank you.

Jie

[Martin Kronbichler]

Hi Jie,

One needs to move some of the vertices around. I have some code that does generate this geometry in one of my projects:

https://github.com/kronbichler/adaflo/blob/master/tests/flow_past_cylinder.cc

(check the create_triangulation function)

The mesh that comes out from this code is accurate but not optimal yet it is not fine enough around the cylinder as compared to be bulk of the domain, at least for some of the benchmark test cases. You can find a picture which gives better convergence results in our recent publication (Fig 11):

http://www.sciencedirect.com/science/article/pii/S0021999117306915

However, I do not have that code online right now and it takes too much time to look it up. But it should be easy enough to compose it of several triangulation, adding the transfinite interpolation manifold to the transition from polar->straight elements that we recently introduced in the deal.II and that is mentioned in that paper.

Best,
Martin