 |
deal.II version GIT relicensing-2250-g88cb8ba3cb 2024-12-13 12:20:00+00:00
|
Loading...
Searching...
No Matches
|
deal.II version GIT relicensing-2250-g88cb8ba3cb 2024-12-13 12:20:00+00:00
|
This tutorial depends on step-22.
This program grew out of a student project by Abner Salgado at Texas A&M University. Most of the work for this program is by him.
The purpose of this program is to show how to effectively solve the incompressible time-dependent Navier-Stokes equations. These equations describe the flow of a viscous incompressible fluid and read
\begin{align*} u_t + u \cdot \nabla u - \nu \Delta u + \nabla p = f, \\ \nabla \cdot u = 0, \end{align*}
where \(u\) represents the velocity of the flow and \(p\) the pressure. This system of equations is supplemented by the initial condition
\[ u |_{t=0} = u_0, \]
with \(u_0\) sufficiently smooth and solenoidal, and suitable boundary conditions. For instance, an admissible boundary condition, is
\[ u|_{\partial\Omega} = u_b. \]
It is possible to prescribe other boundary conditions as well. In the test case that we solve here the boundary is partitioned into two disjoint subsets \(\partial\Omega = \Gamma_1 \cup \Gamma_2\) and we have
\[ u|_{\Gamma_1} = u_b, \]
and
\[ u\times n|_{\Gamma_2} = 0, \quad p|_{\Gamma_2} = 0 \]
where \(n\) is the outer unit normal. The boundary conditions on \(\Gamma_2\) are often used to model outflow conditions.
In previous tutorial programs (see for instance step-20 and step-22) we have seen how to solve the time-independent Stokes equations using a Schur complement approach. For the time-dependent case, after time discretization, we would arrive at a system like
\begin{align*} \frac1\tau u^k - \nu \Delta u^k + \nabla p^k = F^k, \\ \nabla \cdot u^k = 0, \end{align*}
where \(\tau\) is the time-step. Although the structure of this system is similar to the Stokes system and thus it could be solved using a Schur complement approach, it turns out that the condition number of the Schur complement is proportional to \(\tau^{-2}\). This makes the system very difficult to solve, and means that for the Navier-Stokes equations, this is not a useful avenue to the solution.
Rather, we need to come up with a different approach to solve the time-dependent Navier-Stokes equations. The difficulty in their solution comes from the fact that the velocity and the pressure are coupled through the constraint
\[ \nabla \cdot u = 0, \]
for which the pressure is the Lagrange multiplier. Projection methods aim at decoupling this constraint from the diffusion (Laplace) operator.
Let us shortly describe how the projection methods look like in a semi-discrete setting. The objective is to obtain a sequence of velocities \(\{u^k\}\) and pressures \(\{p^k\}\). We will also obtain a sequence \(\{\phi^k\}\) of auxiliary variables. Suppose that from the initial conditions, and an application of a first order method we have found \((u^0,p^0,\phi^0=0)\) and \((u^1,p^1,\phi^1=p^1-p^0)\). Then the projection method consists of the following steps:
\[ u^\star = 2u^k - u^{k-1}, \quad p^\sharp = p^k + \frac43 \phi^k - \frac13 \phi^{k-1}. \]
Step 1: Diffusion step. We find \(u^{k+1}\) that solves the single linear equation
\[ \frac1{2\tau}\left( 3u^{k+1} - 4u^k + u^{k-1} \right) + u^\star \cdot\nabla u^{k+1} + \frac12 \left( \nabla \cdot u^\star \right) u^{k+1} -\nu \Delta u^{k+1} + \nabla p^\sharp = f^{k+1}, \quad u^{k+1}|_{\Gamma_1} = u_b, \quad u^{k+1} \times n|_{\Gamma_2} = 0. \]
\[ \Delta \phi^{k+1} = \frac3{2\tau} \nabla \cdot u^{k+1}, \quad \partial_n \phi^{k+1}|_{\Gamma_1} = 0, \quad \phi^{k+1}|_{\Gamma_2} = 0 \]
\[ p^{k+1} = p^k + \phi^{k+1}. \]
\[ p^{k+1} = p^k + \phi^{k+1} - \nu \nabla \cdot u^{k+1}. \]
Without going into details, let us remark a few things about the projection methods that we have just described:
\[ u \cdot \nabla u + \frac12 \left( \nabla\cdot u \right) u. \]
This is consistent with the continuous equation (because \(\nabla\cdot u = 0\), though this is not true pointwise for the discrete solution) and it is needed to guarantee unconditional stability of the time-stepping scheme. Moreover, to linearize the term we use the second order extrapolation \(u^\star\) of \(u^{k+1}\).\[ L^2(\Omega)^d = H \oplus \nabla H^1_{\Gamma_2}(\Omega), \]
where\[ H = \left\{ v \in L^2(\Omega)^d:\ \nabla\cdot v =0, \ v\cdot n|_{\Gamma_1} = 0 \right\}, \]
and\[ H^1_{\Gamma_2}(\Omega) = \left\{ q \in H^1(\Omega):\ q|_{\Gamma_2} = 0 \right\}. \]
Indeed, if we use this decomposition on \(u^{k+1}\) we obtain\[ u^{k+1} = v^{k+1} + \nabla \left( \frac{2\tau}{3} \phi^{k+1} \right), \]
with \(v^{k+1}\in H\). Taking the divergence of this equation we arrive at the projection equation.The standard incremental scheme and the rotational incremental scheme were first considered by van Kan in
and is analyzed by Guermond in
for the case \(\nu = 1\). It turns out that this technique suffers from unphysical boundary conditions for the kinematic pressure that lead to reduced rates of convergence. To prevent this, Timmermans et al. proposed in
the rotational pressure-correction projection method that uses a divergence correction for the kinematic pressure. A thorough analysis for scheme has first been performed in
for the Stokes problem.
To obtain a fully discrete setting of the method we, as always, need a variational formulation. There is one subtle issue here given the nature of the boundary conditions. When we multiply the equation by a suitable test function one of the term that arises is
\[ -\nu \int_\Omega \Delta u \cdot v. \]
If we, say, had Dirichlet boundary conditions on the whole boundary then after integration by parts we would obtain
\[ -\nu \int_\Omega \Delta u \cdot v = \nu \int_\Omega \nabla u : \nabla v - \int_{\partial\Omega} \partial_n u \cdot v = \nu \int_\Omega \nabla u : \nabla v. \]
One of the advantages of this formulation is that it fully decouples the components of the velocity. Moreover, they all share the same system matrix. This can be exploited in the program.
However, given the nonstandard boundary conditions, to be able to take them into account we need to use the following identity
\[ \Delta u = \nabla\nabla\cdot u - \nabla\times\nabla\times u, \]
so that when we integrate by parts and take into account the boundary conditions we obtain
\[ -\nu \int_\Omega \Delta u \cdot v = \nu \int_\Omega \left[ \nabla \cdot u \nabla \cdot v + \nabla \times u \nabla \times v \right], \]
which is the form that we would have to use. Notice that this couples the components of the velocity. Moreover, to enforce the boundary condition on the pressure, we need to rewrite
\[ \int_\Omega \nabla p \cdot v = -\int_\Omega p \nabla \cdot v + \int_{\Gamma_1} p v\cdot n + \int_{\Gamma_2} p v\cdot n = -\int_\Omega p \nabla \cdot v, \]
where the boundary integral in \(\Gamma_1\) equals zero given the boundary conditions for the velocity, and the one in \(\Gamma_2\) given the boundary conditions for the pressure.
In the simplified case where the boundary \(\Gamma_2\) is parallel to a coordinate axis, which holds for the testcase that we carry out below, it can actually be shown that
\[ \nu \int_\Omega \nabla u : \nabla v = \nu \int_\Omega \left[ \nabla \cdot u \nabla \cdot v + \nabla \times u \nabla \times v \right]. \]
This issue is not very often addressed in the literature. For more information the reader can consult, for instance,
Our implementation of the projection methods follows verbatim the description given above. We must note, however, that as opposed to most other problems that have several solution components, we do not use vector-valued finite elements. Instead, we use separate finite elements the components of the velocity and the pressure, respectively, and use different DoFHandler's for those as well. The main reason for doing this is that, as we see from the description of the scheme, the dim components of the velocity and the pressure are decoupled. As a consequence, the equations for all the velocity components look all the same, have the same system matrix, and can be solved in parallel. Obviously, this approach has also its disadvantages. For instance, we need to keep several DoFHandlers and iterators synchronized when assembling matrices and right hand sides; obtaining quantities that are inherent to vector-valued functions (e.g. divergences) becomes a little awkward, and others.
The testcase that we use for this program consists of the flow around a square obstacle. The geometry is as follows:
with \(H=4.1\), making the geometry slightly non-symmetric.
We impose no-slip boundary conditions on both the top and bottom walls and the obstacle. On the left side we have the inflow boundary condition
\[ u = \left( \begin{array}{c} 4 U_m y (H-y)/H^2 \\ 0 \end{array} \right), \]
with \(U_m = 1.5\), i.e. the inflow boundary conditions correspond to Poiseuille flow for this configuration. Finally, on the right vertical wall we impose the condition that the vertical component of the velocity and the pressure should both be zero. The final time \(T=10\).
We start by including all the necessary deal.II header files and some C++ related ones. Each one of them has been discussed in previous tutorial programs, so we will not get into details here.
Finally this is as in all previous programs:
Since our method has several parameters that can be fine-tuned we put them into an external file, so that they can be determined at run-time.
This includes, in particular, the formulation of the equation for the auxiliary variable \(\phi\), for which we declare an enum. Next, we declare a class that is going to read and store all the parameters that our program needs to run.
In the constructor of this class we declare all the parameters. The details of how this works have been discussed elsewhere, for example in step-29.
In the next namespace, we declare the initial and boundary conditions:
As we have chosen a completely decoupled formulation, we will not take advantage of deal.II's capabilities to handle vector valued problems. We do, however, want to use an interface for the equation data that is somehow dimension independent. To be able to do that, our functions should be able to know on which spatial component we are currently working, and we should be able to have a common interface to do that. The following class is an attempt in that direction.
With this class defined, we declare classes that describe the boundary conditions for velocity and pressure:
NavierStokesProjection classNow for the main class of the program. It implements the various versions of the projection method for Navier-Stokes equations. The names for all the methods and member variables should be self-explanatory, taking into account the implementation details given in the introduction.
The next few structures and functions are for doing various things in parallel. They follow the scheme laid out in Parallel computing with multiple processors accessing shared memory, using the WorkStream class. As explained there, this requires us to declare two structures for each of the assemblers, a per-task data and a scratch data structure. These are then handed over to functions that assemble local contributions and that copy these local contributions to the global objects.
One of the things that are specific to this program is that we don't just have a single DoFHandler object that represents both the velocities and the pressure, but we use individual DoFHandler objects for these two kinds of variables. We pay for this optimization when we want to assemble terms that involve both variables, such as the divergence of the velocity and the gradient of the pressure, times the respective test functions. When doing so, we can't just anymore use a single FEValues object, but rather we need two, and they need to be initialized with cell iterators that point to the same cell in the triangulation but different DoFHandlers.
To do this in practice, we declare a "synchronous" iterator – an object that internally consists of several (in our case two) iterators, and each time the synchronous iteration is moved forward one step, each of the iterators stored internally is moved forward one step as well, thereby always staying in sync. As it so happens, there is a deal.II class that facilitates this sort of thing. (What is important here is to know that two DoFHandler objects built on the same triangulation will walk over the cells of the triangulation in the same order.)
The same general layout also applies to the following classes and functions implementing the assembly of the advection term:
The final few functions implement the diffusion solve as well as postprocessing the output, including computing the curl of the velocity:
NavierStokesProjection::NavierStokesProjection In the constructor, we just read all the data from the Data_Storage object that is passed as an argument, verify that the data we read is reasonable and, finally, create the triangulation and load the initial data.
NavierStokesProjection::create_triangulation_and_dofsThe method that creates the triangulation and refines it the needed number of times. After creating the triangulation, it creates the mesh dependent data, i.e. it distributes degrees of freedom and renumbers them, and initializes the matrices and vectors that we will use.
NavierStokesProjection::initialize This method creates the constant matrices and loads the initial data
NavierStokesProjection::initialize_*_matrices In this set of methods we initialize the sparsity patterns, the constraints (if any) and assemble the matrices that do not depend on the timestep dt. Note that for the Laplace and mass matrices, we can use functions in the library that do this. Because the expensive operations of this function – creating the two matrices – are entirely independent, we could in principle mark them as tasks that can be worked on in parallel using the Threads::new_task functions. We won't do that here since these functions internally already are parallelized, and in particular because the current function is only called once per program run and so does not incur a cost in each time step. The necessary modifications would be quite straightforward, however.
The initialization of the matrices that act on the pressure space is similar to the ones that act on the velocity space.
For the gradient operator, we start by initializing the sparsity pattern and compressing it. It is important to notice here that the gradient operator acts from the pressure space into the velocity space, so we have to deal with two different finite element spaces. To keep the loops synchronized, we use the alias that we have defined before, namely PairedIterators and IteratorPair.
NavierStokesProjection::run This is the time marching function, which starting at t_0 advances in time using the projection method with time step dt until T.
Its second parameter, verbose indicates whether the function should output information what it is doing at any given moment: for example, it will say whether we are working on the diffusion, projection substep; updating preconditioners etc. Rather than implementing this output using code like
we use the ConditionalOStream class to do that for us. That class takes an output stream and a condition that indicates whether the things you pass to it should be passed through to the given output stream, or should just be ignored. This way, above code simply becomes
and does the right thing in either case.
NavierStokesProjection::diffusion_stepThe implementation of a diffusion step. Note that the expensive operation is the diffusion solve at the end of the function, which we have to do once for each velocity component. To accelerate things a bit, we allow to do this in parallel, using the Threads::new_task function which makes sure that the dim solves are all taken care of and are scheduled to available processors: if your machine has more than one processor core and no other parts of this program are using resources currently, then the diffusion solves will run in parallel. On the other hand, if your system has only one processor core then running things in parallel would be inefficient (since it leads, for example, to cache congestion) and things will be executed sequentially.
NavierStokesProjection::assemble_advection_term The following few functions deal with assembling the advection terms, which is the part of the system matrix for the diffusion step that changes at every time step. As mentioned above, we will run the assembly loop over all cells in parallel, using the WorkStream class and other facilities as described in the documentation topic on Parallel computing with multiple processors accessing shared memory.
NavierStokesProjection::projection_stepThis implements the projection step:
NavierStokesProjection::update_pressure This is the pressure update step of the projection method. It implements the standard formulation of the method, that is
\[ p^{n+1} = p^n + \phi^{n+1}, \]
or the rotational form, which is
\[ p^{n+1} = p^n + \phi^{n+1} - \frac{1}{Re} \nabla\cdot u^{n+1}. \]
NavierStokesProjection::output_results This method plots the current solution. The main difficulty is that we want to create a single output file that contains the data for all velocity components, the pressure, and also the vorticity of the flow. On the other hand, velocities and the pressure live on separate DoFHandler objects, and so can't be written to the same file using a single DataOut object. As a consequence, we have to work a bit harder to get the various pieces of data into a single DoFHandler object, and then use that to drive graphical output.
We will not elaborate on this process here, but rather refer to step-32, where a similar procedure is used (and is documented) to create a joint DoFHandler object for all variables.
Let us also note that we here compute the vorticity as a scalar quantity in a separate function, using the \(L^2\) projection of the quantity \(\text{curl} u\) onto the finite element space used for the components of the velocity. In principle, however, we could also have computed it as a pointwise quantity from the velocity, and do so through the DataPostprocessor mechanism discussed in step-29 and step-33.
Following is the helper function that computes the vorticity by projecting the term \(\text{curl} u\) onto the finite element space used for the components of the velocity. The function is only called whenever we generate graphical output, so not very often, and as a consequence we didn't bother parallelizing it using the WorkStream concept as we do for the other assembly functions. That should not be overly complicated, however, if needed. Moreover, the implementation that we have here only works for 2d, so we bail if that is not the case.
The main function looks very much like in all the other tutorial programs, so there is little to comment on here:
We run the code with the following parameter-file.prm, which can be found in the same directory as the source:
# First a global definition
# the type of method we want to use
set Method_Form = rotational
subsection Physical data
# In this subsection we declare the physical data
# The initial and final time, and the Reynolds number
set initial_time = 0.
set final_time = 25.
set Reynolds = 100
end
subsection Time step data
# In this subsection we declare the data that is to be used for time discretization,
# i.e. the time step dt
set dt = 5e-3
end
subsection Space discretization
# In this subsection we declare the data that is relevant to the space discretization
# we set the number of global refines the triangulation must have
# and the degree k of the pair Q_(k+1)--Q_k of velocity--pressure finite element spaces
set n_of_refines = 3
set pressure_fe_degree = 1
end
subsection Data solve velocity
# In this section we declare the parameters that are going to control the solution process
# for the velocity.
set max_iterations = 1000 # maximal number of iterations that GMRES must make
set eps = 1e-6 # stopping criterion
set Krylov_size = 30 # size of the Krylov subspace to be used in GMRES
set off_diagonals = 60 # number of off diagonals that ILU must compute
set diag_strength = 0.01 # diagonal strengthening value
set update_prec = 10 # this number indicates how often the preconditioner must be updated
end
#The output frequency
set output_interval = 50
#Finally we set the verbosity level
set verbose = false
Since the verbose parameter is set to false, we do not get any kind of output besides the number of the time step the program is currently working on. If we were to set it to true we would get information on what the program is doing and how many steps each iterative process had to make to converge, etc.
Let us plot the obtained results for \(t=1,5,12,20,25\) (i.e. time steps 200, 1000, 2400, 4000, and 5000), where in the left column we show the vorticity and in the right the velocity field:
 |  |
 |  |
 |  |
 |  |
 |  |
The images show nicely the development and extension of a vortex chain behind the obstacles, with the sign of the vorticity indicating whether this is a left or right turning vortex.
We can change the Reynolds number, \(Re\), in the parameter file to a value of \(500\). Doing so, and reducing the time step somewhat as well, yields the following images at times \(t=20,40\):
 |  |
 |  |
For this larger Reynolds number, we observe unphysical oscillations, especially for the vorticity. The discretization scheme has now difficulties in correctly resolving the flow, which should still be laminar and well-organized. These phenomena are typical of discretization schemes that lack robustness in under-resolved scenarios, where under-resolved means that the Reynolds number computed with the mesh size instead of the physical dimensions of the geometry is large. We look at a zoom at the region behind the obstacle, and the mesh size we have there:
We can easily test our hypothesis by re-running the simulation with one more mesh refinement set in the parameter file :
Indeed, the vorticity field now looks much smoother. While we can expect that further refining the mesh will suppress the remaining oscillations as well, one should take measures to obtain a robust scheme in the limit of coarse resolutions, as described below.
This program can be extended in the following directions:
Adaptive mesh refinement: As we have seen, we computed everything on a single fixed mesh. Using adaptive mesh refinement can lead to increased accuracy while not significantly increasing the computational time.
Adaptive time-stepping: Although there apparently is currently no theory about projection methods with variable time step, practice shows that they perform very well.
High Reynolds numbers: As we can see from the results, increasing the Reynolds number changes significantly the behavior of the discretization scheme. Using well-known stabilization techniques we could be able to compute the flow in this, or many other problems, when the Reynolds number is very large and where computational costs demand spatial resolutions for which the flow is only marginally resolved, especially for 3D turbulent flows.
Variable density incompressible flows: There are projection-like methods for the case of incompressible flows with variable density. Such flows play a role if fluids of different density mix, for example fresh water and salt water, or alcohol and water.
1.9.8