Linear statics – small strain¶
In this first example we will look at a linear elastic, 2-D plane strain problem, with imposed boundary displacements. The problem is sketched in the figure below. From which it is observed that:
- The sample is homogeneous.
- Symmetry conditions are assumed in - and -direction.
- The outer (right) edge is displaced by in -direction (whereby the displacement is constant in -direction).
- The boundary tractions are zero everywhere, except for where the displacement is prescribed. There a non-zero reaction force may appear.
- (Not shown) We assume that , and thus . By doing this we can make the problem truly linear.
The geometry is discretized in 2-D linear quadrilateral elements, which have four nodes per element. To keep things simple we use only three, equi-sized, elements in each direction. The mesh is shown below, whereby the element numbers have a regular font while the node numbers are italic.
At this point we focus our attention on the internal force. Thereby we first consider the constitutive response:
I.e the double contraction between the fourth-order stiffness
(with and the bulk and the shear modulus), and the linear strain
(more information in the documentation of GooseSolid). Because of the symmetries in we can simplify the constitutive expression as follows
The displacement of the final configuration, , is now decomposed some known pre-strain plus an unknown update :
For our simple problem, which is initially stress and strain free, we find that
Some of the above expressions could thus be simplified further, while also part of the implementation can be omitted. We keep it here to build on this problem later on.
For the internal force this implies that
Since is known we easily evaluate the original (in principle non-linear) expression of the internal force:
For the update we use the explicit relation for the stress
We then again apply our discretization scheme to obtain
The tractions are fully given by the boundary conditions, which also can be used to show that – or at least for the relevant components (see below). The final force balance then reads:
We continue, by writing the problem in terms of scalar degrees-of-freedoms (DOFs). Each node has two DOFs, the two vector directions. For our mesh we define the DOFs as follows (whereby the cyan colored DOF-numbers correspond to the -direction, while the magenta ones correspond to the -direction).
By doing this, a little bit of book-keeping allows us to write our balance equation as the following system of equations
Next, we have to impose the boundary conditions. Given the fact that we impose displacements on part of the boundary here, a part of will be prescribed. More specifically, our mesh looks as follows
where the yellow nodes are prescribed DOFs while the blue ones are yet unknown. We are now ready to solve the balance equation in two steps. First we will partition the system is unknown and prescribed DOFs:
We are now ready to solve the former part:
which gives us that
and finally that
which can be reassembled as displacement vectors per node, .
Should you be interested, one can compute the reaction force, i.e. the boundary tractions there where the displacement has been prescribed. One has to do the following:
For this, one thus has to compute the new . Because this specific model is linear we can however obtain the reaction forces without having to re-evaluate . Specifically
Our first attempt of an implementation literately follows the steps above: it constructs and , which are then partitioned. The prescribed displacements are then set. Thereafter the problem is solved, and the displacements are reconstructed to nodal vector for easy post-processing.
One of the things that made the previous examples not very suitable for becoming a production code is the fact that the stiffness matrix was first fully assembled and afterwards partitioned. Besides costing a lot of memory for storing the matrix twice, it might cost a lot of time since partitioning might become a costly operation in the case that sparse matrices are used. To avoid this, the system may be pre-partioned. In that case we renumber the DOFs such that we end up with first all the unknown DOFs (denoted
iiu in the code), and then all the known DOFs (denoted
iip). For our example this results in:
This allows us to consider four different matrices (denoted using
_pp) and two different columns (denoted using
_p) to which the internal force and the stiffness are directly assembled. The rows (and columns) of these matrices and columns follow from introducing separate indices for
Some additional notes on the theory discussed on a simplified scalar system, for the same mesh as presented here, are included in a separate document. One is invited to study this document before continuing.
In our first example we will consider the same material and mesh as above. However, now we will assume periodicity is both spatial directions and prescribe a change in the macroscopic deformation gradient, equal to
First of all we will specify periodicity for our mesh. It applies that the following equalities hold (in terms of node numbers:
where are the microscopic fluctuations, that do not affect the macroscopic affine deformation. In terms of DOFs this is can be illustrated as follows:
where the red DOFs are said to be dependent (i.e. they directly follow from the equalities listed above). The simplest this that we can do is construct a system with only the independent DOFs (in blue above) by directly assembling to the independent DOFs. To this end we employ the following DOF numbers:
where the yellow color of the lower left corner indicates that this node is used as reference. Firstly it is used to suppress rigid body deformation. Secondly we apply the macroscopic deformation as the initial condition.
Final equilibrium is then obtained by solving
(which has the dimensions of the number of independent DOFs), and then assembling to the entire system (including the dependent nodes). This is done in the first example, whereby the resulting system is partitioned to deal with the zero-displacement of the reference node.
Here we will enable mixed macroscopic boundary conditions by introduction extra DOFs for the macroscopic deformation gradient tensor, and its antagonist stress response. Since we work in two dimensions we introduce two virtual nodes, each with two DOFs:
We will now employ the following tying relation
To this end we first renumber the system to have all the dependent DOFs at the end
And then partition the system in independent and dependent DOFs
Finally, we obtain the following tying relations for the DOFs
We then employ the nodal dependency to obtain a system of the independent DOFs only:
which, after solving, we can reconstruct to the dependent DOFs using
Towards a production code
Although not (yet) pursued here, it would make sense to partition the system as follows:
Accompanied with the tying relations
And use the following condensation
with the following reconstruction: