From 655e4542400ffb2e3e0ae17cd3d6d8fe23d402a6 Mon Sep 17 00:00:00 2001 From: bangerth Date: Fri, 3 Jun 2011 03:00:38 +0000 Subject: [PATCH] Expand discussion and move possibilities for extensions to the end, where they usually are. git-svn-id: https://svn.dealii.org/trunk@23784 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-18/doc/intro.dox | 125 --------------- deal.II/examples/step-18/doc/results.dox | 187 +++++++++++++++++++++++ 2 files changed, 187 insertions(+), 125 deletions(-) diff --git a/deal.II/examples/step-18/doc/intro.dox b/deal.II/examples/step-18/doc/intro.dox index 10940ee4cd..d912f3b072 100644 --- a/deal.II/examples/step-18/doc/intro.dox +++ b/deal.II/examples/step-18/doc/intro.dox @@ -555,131 +555,6 @@ the results section will show what displacement pattern comes out of this simulation. -

Possible directions for extensions

- -The program as is does not really solve an equation that has many applications -in practice: quasi-static material deformation based on a purely elastic law -is almost boring. However, the program may serve as the starting point for -more interesting experiments, and that indeed was the initial motivation for -writing it. Here are some suggestions of what the program is missing and in -what direction it may be extended: - -
Plasticity models
- - The most obvious extension is to use a more -realistic material model for large-scale quasistatic deformation. The natural -choice for this would be plasticity, in which a nonlinear relationship between -stress and strain replaces equation [stress-strain]. Plasticity -models are usually rather complicated to program since the stress-strain -dependence is generally non-smooth. The material can be thought of being able -to withstand only a maximal stress (the yield stress) after which it starts to -“flow”. A mathematical description to this can be given in the form of a -variational inequality, which alternatively can be treated as minimizing the -elastic energy -@f[ - E(\mathbf{u}) = - (\varepsilon(\mathbf{u}), C\varepsilon(\mathbf{u}))_{\Omega} - - (\mathbf{f}, \mathbf{u})_{\Omega} - (\mathbf{b}, \mathbf{u})_{\Gamma_N}, -@f] -subject to the constraint -@f[ - f(\sigma(\mathbf{u})) \le 0 -@f] -on the stress. This extension makes the problem to be solved in each time step -nonlinear, so we need another loop within each time step. - -Without going into further details of this model, we refer to the excellent -book by Simo and Hughes on “Computational Inelasticity” for a -comprehensive overview of computational strategies for solving plastic -models. Alternatively, a brief but concise description of an algorithm for -plasticity is given in an article by S. Commend, A. Truty, and Th. Zimmermann, -titled “Stabilized finite elements applied to -elastoplasticity: I. Mixed displacement-pressure formulation” -(Computer Methods in Applied Mechanics and Engineering, vol. 193, -pp. 3559-3586, 2004). - - -
Stabilization issues
- -The formulation we have chosen, i.e. using -piecewise (bi-, tri-)linear elements for all components of the displacement -vector, and treating the stress as a variable dependent on the displacement is -appropriate for most materials. However, this so-called displacement-based -formulation becomes unstable and exhibits spurious modes for incompressible or -nearly-incompressible materials. While fluids are usually not elastic (in most -cases, the stress depends on velocity gradients, not displacement gradients, -although there are exceptions such as electro-rheologic fluids), there are a -few solids that are nearly incompressible, for example rubber. Another case is -that many plasticity models ultimately let the material become incompressible, -although this is outside the scope of the present program. - -Incompressibility is characterized by Poisson's ratio -@f[ - \nu = \frac{\lambda}{2(\lambda+\mu)}, -@f] -where $\lambda,\mu$ are the Lam\'e constants of the material. -Physical constraints indicate that $-1\le \nu\le \frac 12$ (the condition -also follows from mathematical stability considerations). If $\nu$ -approaches $\frac 12$, then the material becomes incompressible. In that -case, pure displacement-based formulations are no longer appropriate for the -solution of such problems, and stabilization techniques have to be employed -for a stable and accurate solution. The book and paper cited above give -indications as to how to do this, but there is also a large volume of -literature on this subject; a good start to get an overview of the topic can -be found in the references of the paper by -H.-Y. Duan and Q. Lin on “Mixed finite elements of least-squares type for -elasticity” (Computer Methods in Applied Mechanics and Engineering, vol. 194, -pp. 1093-1112, 2005). - - -
Refinement during timesteps
- -In the present form, the program -only refines the initial mesh a number of times, but then never again. For any -kind of realistic simulation, one would want to extend this so that the mesh -is refined and coarsened every few time steps instead. This is not hard to do, -in fact, but has been left for future tutorial programs or as an exercise, if -you wish. The main complication one has to overcome is that one has to -transfer the data that is stored in the quadrature points of the cells of the -old mesh to the new mesh, preferably by some sort of projection scheme. This -is only slightly messy in the sequential case; in fact, the functions -FETools::get_projection_from_quadrature_points_matrix will do -the projection, and the FiniteElement::get_restriction_matrix and -FiniteElement::get_prolongation_matrix functions will do the -transfer between mother and child cells. Alternatively, one can use -FETools::get_projection_from_quadrature_points_matrix to convert data only -defined in quadrature points to a finite element field, use the -SolutionTransfer class to move the field to a refined/coarsened version of the -mesh, and then use FETools::compute_interpolation_to_quadrature_points_matrix -to get the information back into the quadrature points of the new mesh. - -However, whichever way we do it, it becomes complicated -once we run the program in parallel, since then each process only stores this -data for the cells it owned on the old mesh, and it may need to know the -values of the quadrature point data on other cells if the corresponding cells -on the new mesh are assigned to this process after subdividing the new mesh. A -global communication of these data elements is therefore necessary, making the -entire process a little more unpleasant. The step-28 tutorial -program shows how to work with different meshes at the same time, albeit for a -different kind of problem, giving indications on how to approach the problem -for time-dependent adaptivity as well. - - -
Ensuring mesh regularity
- -At present, the program makes no attempt -to make sure that a cell, after moving its vertices at the end of the time -step, still has a valid geometry (i.e. that its Jacobian determinant is -positive and bounded away from zero everywhere). It is, in fact, not very hard -to set boundary values and forcing terms in such a way that one gets distorted -and inverted cells rather quickly. Certainly, in some cases of large -deformation, this is unavoidable with a mesh of finite mesh size, but in some -other cases this should be preventable by appropriate mesh refinement and/or a -reduction of the time step size. The program does not do that, but a more -sophisticated version definitely should employ some sort of heuristic defining -what amount of deformation of cells is acceptable, and what isn't. - - @note PETSc appears not to co-operate well when using threads and the program crashes when deal.II is compiled in its usual mode supporting multithreading. The program can therefore only be compiled and run diff --git a/deal.II/examples/step-18/doc/results.dox b/deal.II/examples/step-18/doc/results.dox index 4a80e0ad4a..8bb54ed7c8 100644 --- a/deal.II/examples/step-18/doc/results.dox +++ b/deal.II/examples/step-18/doc/results.dox @@ -334,3 +334,190 @@ simulation we would also need to make sure that the body does not intersect itself during deformation. Without such a formulation we cannot expect anything to make physical sense, even if it produces nice pictures! + +

Possible directions for extensions

+ +The program as is does not really solve an equation that has many applications +in practice: quasi-static material deformation based on a purely elastic law +is almost boring. However, the program may serve as the starting point for +more interesting experiments, and that indeed was the initial motivation for +writing it. Here are some suggestions of what the program is missing and in +what direction it may be extended: + +
Plasticity models
+ + The most obvious extension is to use a more +realistic material model for large-scale quasistatic deformation. The natural +choice for this would be plasticity, in which a nonlinear relationship between +stress and strain replaces equation [stress-strain]. Plasticity +models are usually rather complicated to program since the stress-strain +dependence is generally non-smooth. The material can be thought of being able +to withstand only a maximal stress (the yield stress) after which it starts to +“flow”. A mathematical description to this can be given in the form of a +variational inequality, which alternatively can be treated as minimizing the +elastic energy +@f[ + E(\mathbf{u}) = + (\varepsilon(\mathbf{u}), C\varepsilon(\mathbf{u}))_{\Omega} + - (\mathbf{f}, \mathbf{u})_{\Omega} - (\mathbf{b}, \mathbf{u})_{\Gamma_N}, +@f] +subject to the constraint +@f[ + f(\sigma(\mathbf{u})) \le 0 +@f] +on the stress. This extension makes the problem to be solved in each time step +nonlinear, so we need another loop within each time step. + +Without going into further details of this model, we refer to the excellent +book by Simo and Hughes on “Computational Inelasticity” for a +comprehensive overview of computational strategies for solving plastic +models. Alternatively, a brief but concise description of an algorithm for +plasticity is given in an article by S. Commend, A. Truty, and Th. Zimmermann, +titled “Stabilized finite elements applied to +elastoplasticity: I. Mixed displacement-pressure formulation” +(Computer Methods in Applied Mechanics and Engineering, vol. 193, +pp. 3559-3586, 2004). + + +
Stabilization issues
+ +The formulation we have chosen, i.e. using +piecewise (bi-, tri-)linear elements for all components of the displacement +vector, and treating the stress as a variable dependent on the displacement is +appropriate for most materials. However, this so-called displacement-based +formulation becomes unstable and exhibits spurious modes for incompressible or +nearly-incompressible materials. While fluids are usually not elastic (in most +cases, the stress depends on velocity gradients, not displacement gradients, +although there are exceptions such as electro-rheologic fluids), there are a +few solids that are nearly incompressible, for example rubber. Another case is +that many plasticity models ultimately let the material become incompressible, +although this is outside the scope of the present program. + +Incompressibility is characterized by Poisson's ratio +@f[ + \nu = \frac{\lambda}{2(\lambda+\mu)}, +@f] +where $\lambda,\mu$ are the Lam\'e constants of the material. +Physical constraints indicate that $-1\le \nu\le \frac 12$ (the condition +also follows from mathematical stability considerations). If $\nu$ +approaches $\frac 12$, then the material becomes incompressible. In that +case, pure displacement-based formulations are no longer appropriate for the +solution of such problems, and stabilization techniques have to be employed +for a stable and accurate solution. The book and paper cited above give +indications as to how to do this, but there is also a large volume of +literature on this subject; a good start to get an overview of the topic can +be found in the references of the paper by +H.-Y. Duan and Q. Lin on “Mixed finite elements of least-squares type for +elasticity” (Computer Methods in Applied Mechanics and Engineering, vol. 194, +pp. 1093-1112, 2005). + + +
Refinement during timesteps
+ +In the present form, the program +only refines the initial mesh a number of times, but then never again. For any +kind of realistic simulation, one would want to extend this so that the mesh +is refined and coarsened every few time steps instead. This is not hard to do, +in fact, but has been left for future tutorial programs or as an exercise, if +you wish. + +The main complication one has to overcome is that one has to +transfer the data that is stored in the quadrature points of the cells of the +old mesh to the new mesh, preferably by some sort of projection scheme. The +general approach to this would go like this: + +- At the beginning, the data is only available in the quadrature points of + individual cells, not as a finite element field that is defined everywhere. + +- So let us find a finite element field that is defined everywhere so + that we can later interpolate it to the quadrature points of the new + mesh. In general, it will be difficult to find a continuous finite element + field that matches the values in the quadrature points exactly because the + number of degrees of freedom of these fields does not match the number of + quadrature points there are, and the nodal values of this global field will + either be over- or underdetermined. But it is usually not very difficult to + find a discontinuous field that matches the values in the quadrature points; + for example, if you have a QGauss(2) quadrature formula (i.e. 4 points per + cell in 2d, 8 points in 3d), then one would use a finite element of kind + FE_DGQ(1), i.e. bi-/tri-linear functions as these have 4 degrees of freedom + per cell in 2d and 8 in 3d. There are functions that can make this + conversion from individual points to a global field simpler; the following + piece of pseudo-code should help if you use a QGauss(2) quadrature formula + (the prefix history_ indicates that we work with quantities + related to the history variables defined in the quadrature points): + @code + FE_DGQ history_fe (1); + DoFHandler history_dof_handler (triangulation); + history_dof_handler.distribute_dofs (history_fe); + + Vector history_field (history_dof_handler.n_dofs()); + + FullMatrix qpoint_to_dof_matrix (history_fe.dofs_per_cell, + quadrature.size()); + FETools::compute_projection_from_quadrature_points_matrix + (history_fe, + quadrature, quadrature, + qpoint_to_dof_matrix); + + Vector local_history_values_at_qpoints (quadrature.size()); + Vector local_history_fe_values (fe.dofs_per_cell); + for (cell=...) + { + ...collect values from quadrature points into + local_history_values_at_qpoints... + qpoint_to_dof_matrix.vmult (local_history_fe_values, + local_history_values_at_qpoints); + cell->set_dof_values (local_history_fe_values, + history_field); + } + @endcode + +- Now that we have a global field, we can refine the mesh and transfer the + history_field vector as usual using the SolutionTransfer class. This will + interpolate everything from the old to the new mesh. + +- In a final step, we have to get the data back from the now interpolated + global field to the quadrature points on the new mesh. The following code + will do that: + @code + FullMatrix dof_to_qpoint_matrix (quadrature.size(), + history_fe.dofs_per_cell); + FETools::compute_interpolation_to_quadrature_points_matrix + (history_fe, + quadrature, + dof_to_qpoint_matrix); + + Vector local_history_values_at_qpoints (quadrature.size()); + Vector local_history_fe_values (fe.dofs_per_cell); + for (cell=...) + { + cell->set_get_values (history_field, + local_history_fe_values); + dof_to_qpoint_matrix.vmult (local_history_values_at_qpoints, + local_history_fe_values); + ...put values back from local_history_values_at_qpoints + quadrature points into... + } + @endcode + +It becomes a bit more complicated once we run the program in parallel, since +then each process only stores this data for the cells it owned on the old +mesh. That said, using a parallel vector for history_field will +do the trick if you put a call to compress after the transfer +from quadrature points into the global vector. + + +
Ensuring mesh regularity
+ +At present, the program makes no attempt +to make sure that a cell, after moving its vertices at the end of the time +step, still has a valid geometry (i.e. that its Jacobian determinant is +positive and bounded away from zero everywhere). It is, in fact, not very hard +to set boundary values and forcing terms in such a way that one gets distorted +and inverted cells rather quickly. Certainly, in some cases of large +deformation, this is unavoidable with a mesh of finite mesh size, but in some +other cases this should be preventable by appropriate mesh refinement and/or a +reduction of the time step size. The program does not do that, but a more +sophisticated version definitely should employ some sort of heuristic defining +what amount of deformation of cells is acceptable, and what isn't. + -- 2.39.5