Re: [deal.II] Time dependent Navier-Stokes system setup

[Guido Kanschat]

Hi there, here a few comments.

This does not look like Newton's method to me: the equation for the update is

DF(u^m) du = F(u^m)

and

u^m+1 = u^m + du

You seem to be missing parts of the residual on the right hand side. Still, the method seems to be consistent at least for the velocity.

If the pressure does not converge: dies the error converge to a constant?

Your effective Reynolds number is determined by the mesh size. The dynamics of the transient problem get easily destroyed by the backward Euler method if your time step is too large.

Best,
Guido

On Mar 2, 2015 7:16 PM, "Lukas Bystricky" <lukasby...@gmail.com> wrote:

Hello,

I am trying to set up a time dependent Navier-Stokes solver. I'm using Dirichlet boundary conditions on the velocities. Currently testing the code using the method of manufactured solutions gives the correct convergence in the velocity, but not in the pressure. When I try to replicate lid driven cavity flow, the solution approaches a solution that is very similar to what one would expect from Stokes flow, even for high Reynold's number.

Below is my procedure:

Using Gunzburger (1989), given u^0, we can use Newton's method to solve a sequence of linear problems:

a(u^m, v^h )+ c(u^m,u^(m-1), v^h)  + c(u^(m-1),u^m,v^h) + b(v^h, p^m) = (f, v^h) + c(u^(m-1),u^(m-1),v^h)
b(u^m, q^h) = 0

where

a(u,v) = nu*int(grad u : grad v)
b(v,q) = -int(q div(v))
c(w,u,v) = int(w * grad(u) * v)

This is for steady flow, so for time dependent flow, if I discretize velocity using backwards Euler I get the linear problems:

(u^m, v^h) + dt*( a(u^m, v^h )+ c(u^m,u^(m-1), v^h)  + c(u^(m-1),u^m,v^h) + b(v^h, p^m) ) = dt*( (f, v^h) + c(u^(m-1),u^(m-1),v^h) ) + (u^(k-1), v^h)
b(u^m, q^h) = 0

where u^(k-1) is the velocity from the previous time step.

Here is how I am setting up system in Deal.ii:

template<int dim>
void NavierStokesSolver<dim>::assemble_system()
{
    printf("Assembling matrix...");
    Timer timer;
    timer.start ();
   
    QGauss<dim>                 quadrature_formula(fe.degree+1);
    const int                   dofs_per_cell = fe.dofs_per_cell;
    const int                   n_q_points = quadrature_formula.size();   
    Vector<double>              rhs_values(dim+1); 
   
    std::vector<Tensor<1,dim> > previous_newton_velocity_values(n_q_points);
    std::vector<Tensor< 2, dim> > previous_newton_velocity_gradients(n_q_points);
   
    std::vector<Tensor<2,dim> > grad_phi_u (dofs_per_cell);
    std::vector<double>         div_phi_u (dofs_per_cell);
    std::vector<double>         phi_p (dofs_per_cell);
    std::vector<Tensor<1,dim> > phi_u (dofs_per_cell);
   
    Vector<double>              cell_rhs(dofs_per_cell);
    FullMatrix<double>          cell_matrix (dofs_per_cell, dofs_per_cell);
    FullMatrix<double>          cell_mass (dofs_per_cell, dofs_per_cell);
   
    std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);   
    const FEValuesExtractors::Vector velocities (0);
    const FEValuesExtractors::Scalar pressure (dim);
    FEValues<dim> fe_values(fe, quadrature_formula,
            update_values | update_gradients | update_JxW_values | update_quadrature_points);   
      
    typename DoFHandler<dim>::active_cell_iterator 
            cell = dof_handler.begin_active(), endc = dof_handler.end();
   
    system_matrix.reinit(sparsity_pattern);
    system_matrix_copy.reinit(sparsity_pattern);
    mass_matrix.reinit(sparsity_pattern);
    system_rhs.reinit(dof_handler.n_dofs());  
   
    for (; cell!=endc; ++cell)
    {
        fe_values.reinit(cell);
        fe_values[velocities].get_function_values(previous_newton_step, previous_newton_velocity_values);       
        fe_values[velocities].get_function_gradients(previous_newton_step, previous_newton_velocity_gradients);
                       
        cell_matrix = 0;
        cell_mass   = 0;
        cell_rhs    = 0;
       
        //calculate cell contribution to system        
        for (int q = 0; q < n_q_points; q++)
        {
            for (int k=0; k<dofs_per_cell; k++)
            {
                grad_phi_u[k] = fe_values[velocities].gradient (k, q);
                div_phi_u[k]  = fe_values[velocities].divergence (k, q);
                phi_p[k]      = fe_values[pressure].value (k, q);
                phi_u[k]      = fe_values[velocities].value (k, q);
            }
           
            for (int i = 0; i < dofs_per_cell; i++)
            {       
                for (int j = 0; j < dofs_per_cell; j++)
                {
                    cell_matrix(i,j) +=
                                (phi_u[i]*phi_u[j]
                                        + (1/input.re)*input.dt*double_contract(grad_phi_u[i],grad_phi_u[j])
                                        + input.dt*phi_u[j]*previous_newton_velocity_gradients[q]*phi_u[i]
                                        + input.dt*previous_newton_velocity_values[q]*grad_phi_u[j]*phi_u[i]
                                        - input.dt*phi_p[j]*div_phi_u[i]
                                        - input.dt*phi_p[i]*div_phi_u[j])
                                        *fe_values.JxW(q);
                   
                    cell_mass(i,j) += phi_u[i]*phi_u[j]*fe_values.JxW(q);                                                       
                }
               
                int equation_i = fe.system_to_component_index(i).first;   
                forcing_function->vector_value(fe_values.quadrature_point(q), rhs_values);
               
                cell_rhs[i] +=
                            input.dt*(fe_values.shape_value(i,q)*rhs_values[equation_i]
                             + previous_newton_velocity_values[q]*previous_newton_velocity_gradients[q]*phi_u[i])
                            *fe_values.JxW(q);                   
            }
        }
       
        cell->get_dof_indices(local_dof_indices);
       
        for (int i=0; i<dofs_per_cell; i++)
        {
            for (int j=0; j<dofs_per_cell; j++)
            {
                system_matrix.add (local_dof_indices[i], local_dof_indices[j], cell_matrix(i,j));
                mass_matrix.add(local_dof_indices[i], local_dof_indices[j], cell_mass(i,j));
            }
   
            system_rhs(local_dof_indices[i]) += cell_rhs(i);
        }
    }   
   
    timer.stop();
    printf("done (%gs)\n", timer());                       
}

This sets up a system matrix, a right hand side and a mass matrix for the velocity. I then solve as follows at each time step:

while (iter == 0 || (residual > 1e-8 && iter < 10))
        {           
            previous_newton_step = solution;
            assemble_system();
            mass_matrix.vmult(tmp, old_solution);//add velocity from previous time step
           
            system_rhs.add(tmp);   
           
            constraints.condense(system_matrix, system_rhs);
             
            solve_newton_step();
           
            residual = compute_residual();
                           
            iter++;
           
            printf("Residual = %e\n", residual);           
        }


Are there any glaring omissions/errors in my procedure? Thanks in advance.

--
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.