From: wolf
This tutorial program is another one in the series on the elasticity problem that we have already started with step-8 and step-17. It extends it into two @@ -294,23 +295,31 @@ possibly time-varying force function
Note that these equations are posed on a domain that changes with time, with the boundary moving according to the @@ -386,17 +395,17 @@ displacements of the points on the boundary. To -complete this system, we have to specify the relationship between the -stress and the strain, as follows: +complete this system, we have to specify the incremental relationship between +the stress and the strain, as follows:(11) | (12) | @@ -447,7 +456,7 @@ where||||||||||
(13) |
in | on | on @@ -557,7 +574,7 @@ such that | |||||||||
@@ -570,7 +587,7 @@ are no boundary forces, i.e.
, and that the deformation of the
body is driven by body forces
alone. It is also worth noting that when integrating by parts, we
would get terms of the form
, but that we replace it with the term involving the
@@ -598,14 +615,14 @@ symmetric gradient
instead of
. Due to
the symmetry of
The system at time step
, to be solved on the old domain
, has exactly the form of a stationary elastic
problem, and is therefore similar to what we have already implemented
@@ -639,7 +656,7 @@ continuous finite elements.
There are differences, however:
| |||||||||||
@@ -753,16 +770,16 @@ course replace this term by numerical quadrature:
where are the quadrature weights and the quadrature points on cell . This should make clear that what we really need is not the stress in itself, but only the values of the stress in the quadrature points on all cells. This, however, is a simpler task: we only have to provide @@ -782,7 +799,7 @@ end of each time step we then only have to evaluate , multiply it by the stress-strain tensor at quadrature point . @@ -809,14 +826,14 @@ due to the fact that we compute on the domain , and then use this displacement increment to both update the stress as well as @@ -825,7 +842,7 @@ move the mesh nodes around to get to on which the next increment is computed. What we have to make sure, in this context, is that moving the @@ -838,10 +855,11 @@ can be understood as follows: locally, the incremental deformation can be decomposed into three parts, a linear translation (the constant part -of the displacement field in the neighborhood of a point), a dilational +of the displacement increment field in the neighborhood of a point), a +dilational component (that part of the gradient of the displacement field that has a nonzero divergence), and a rotation. A linear translation of the material does not affect the stresses that are frozen into it - the stress values are @@ -853,14 +871,14 @@ situation where , with which ). Nevertheless, if the the material was pre-stressed in a certain direction, then this direction will be rotated along with the material. To @@ -869,20 +887,20 @@ this end, we have to define a rotation matrix that describes, in each point the rotation due to the displacement increments. It is not hard to see that the actual dependence of on can only be through the curl of the displacement, rather than the displacement @@ -890,7 +908,7 @@ itself or its full gradient (as mentioned above, the constant components of the increment describe translations, its divergence the dilational modes, and the curl the rotational modes). Since the exact form of is cumbersome, we only state it in the program code, and note that the correct updating formula @@ -900,7 +918,7 @@ for the stress variable is then | |||||||||||
(20) | |||||||||||
(21) | (22) | @@ -1052,11 +1071,11 @@ In practice, with quadrature points and weights . We have built these contributions before, in step-8 and step-17, but in both of these cases we @@ -1071,7 +1090,7 @@ with quadrature points . This is not really convenient, in particular if we want to consider more complicated elasticity @@ -1086,7 +1105,7 @@ with quadrature points . While we in fact do not use a more complicated @@ -1111,23 +1130,23 @@ What we then need is two things: a way to create the stress-strain rank-4 as well as to create a symmetric tensor of rank 2 (the strain tensor) from the gradients of a shape function at a quadrature point on a given cell. At the top of the implementation of this example program, you will find such functions. The first one, get_stress_strain_tensor, takes two arguments corresponding to the Lamé constants and and returns the stress-strain tensor for the isotropic case corresponding to these constants (in the program, we @@ -1137,11 +1156,11 @@ What we then need is two things: a way to create the stress-strain rank-4 get_strain takes an object of type FEValues and indices and and returns the symmetric gradient, i.e. the strain, corresponding to shape function , evaluated on the cell on which the FEValues object was last reinitialized. @@ -1193,17 +1212,18 @@ Assembling the local contributions||||||||||
(23) |
@@ -1346,7 +1368,8 @@ 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 (11). Plasticity +stress and strain replaces equation (). 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 @@ -1358,7 +1381,7 @@ elastic energy