Hello everybody,
I would be really gratefully if you could help me.
I have to solve the incompressible Navier-Stokes equations in a square
domain [0,1]x[0x1].
I calculate velocity and pressure fields with a decoupled scheme, using
Q2Q1 pair of elements as spatial discretization. As temporal
discretization I use Backward Euler Scheme, so it is a scheme of first
order.
In order to do an error analysis, I simulate an analytical case where
the solution, both velocity and pressure,depends on time, something like
"sin(a*pi*x)*cos(b*pi*y)*exp(-c*t)".
I use a grid of 2^4*2^4 elements, i.e. a spatial step=h=1/16, another
one of 2^5*2^5 elements, i.e. h=1/32, another one of 2^6*2^6, i.e.
h=1/64 and another grid of 2^7*2^7, i.e. h=1/128.
I use the pair element Q2Q1.
I use four different temporal steps: incrT=1/64,1/256,1/512 and 1/1024.
************************************************************************************************************************
1) Then, I calculate the error in velocity and pressure with L^2 norm in
each iteration in the following way:
Vector<float> difference_per_cell_u (triangulation.n_active_cells());
Vector<float> difference_per_cell_v (triangulation.n_active_cells());
Vector<float> difference_per_cell_p (triangulation.n_active_cells());
ComponentSelectFunction<DIMENSION> seleccionu (0,1,3), seleccionv
(1,1,3), seleccionp (2,1,3);
//Bvector_solution is the numerical solution in the present iteration.
//Kim_Moin2d_final... is the analytical solution.
VectorTools::integrate_difference(dof_handler_total, Bvector_solution,
mi_Functions::Kim_Moin2d_final<DIMENSION>
(DIMENSION+1,viscosity,(num)*time_step_value),
difference_per_cell_u,
QGauss<DIMENSION>(3),VectorTools::L2_norm,
& seleccionu);
VectorTools::integrate_difference(dof_handler_total, Bvector_solution,
mi_Functions::Kim_Moin2d_final<DIMENSION>
(DIMENSION+1,viscosity,(num)*time_step_value),
difference_per_cell_v,
QGauss<DIMENSION>(3),VectorTools::L2_norm,
& seleccionv);
VectorTools::integrate_difference(dof_handler_total,
Bvector_solution,mi_Functions::Kim_Moin2d_final
<DIMENSION>(DIMENSION+1,viscosity,(num)*time_step_value),
difference_per_cell_p,
QGauss<DIMENSION>(3),VectorTools::L2_norm,
& seleccionp);
error_u(iteration)=sqrt((difference_per_cell_u.l2_norm()*difference_per_cell_u.l2_norm())
+
(difference_per_cell_v.l2_norm()*difference_per_cell_v.l2_norm()));
error_p(iteration)=difference_per_cell_p.l2_norm();
2) Once I have all errors for all time steps during the total time to
simulate, I also want to obtain the error in time with L² norm, so I do:
for (unsigned int i=0;i<num_iterations+1;i++)
error_velocidad+=(
error_u(i)*error_u(i)
);//L2(t)
error_velocidad/=(num_iterations+1);
error_velocidad=sqrt(error_velocidad);
for (unsigned int i=0;i<num_iterations+1;i++)
error_presion+=(
error_p(i)*error_p(i)
);
error_presion/=(num_iterations+1);
error_presion=sqrt(error_presion);
****************************************************************************************************************
I think I have programmed everything in the right way. But the problem
is that I don't obtain a good result.
Well, the pressure error behaves in a good way since less refined the
mesh, higher is the error.
But, the velocity error doesn't follow this pattern, for example, I obtain:
Always with Q2Q1:
a)With incrT=1/64.
(error with h=1/32) is the lowest.
(error with h=1/128)>(error with h=1/64)>(error with h=1/32)
b)With incrT=1/256.
(error with h=1/32) is the lowest.
(error with h=1/128)>(error with h=1/64)>(error with h=1/32)
c)With incrT=1/512.
(error with h=1/64) is the lowest
(error with h=1/128)>(error with h=1/64)
d)With incrT=1/1024.
(error with h=1/64) is the lowest
(error with h=1/128)>(error with h=1/64)
Looking at the final solution and comparing it to the exact one, all the
previous situations seem to reach a good agreement.
i) I have thought that perhaps I should try to implement a temporal
discretization scheme of order two, such as Runge-Kutta 2, instead a
Backward Euler, but I don't know if this would solve the problem.
ii) As I am using Q2 as element for the velocity, I had also thought to
use a MappingQ of order 2, but as I read in tutorial and some mails in
the mailing list, it is only used when we are working with curved
domain, and this is not the case.
iii) Besides, the most surprising thing is that the pressure (Q1
element) error has the properly behaviour.
Any advice will be welcome.
Thanks in advance.
Isa
_______________________________________________
