This subroutine reads the Dirichlet boundary values from
the .dat file.
- stiff
- This subroutine calculates the element stiffness matrices
and writes them to the scratch file. It also calculates the right-hand
side load vector F.
- assmb
- This subroutine assembles the element stiffness matrices
to form the global stiffness matrix.
- solve
- This subroutine solves the global stiffness equation
(3.8).
Some more detail on each of these subroutines is now given. Note that the
source code is provided with comment lines explaining the sections, and
describing the variables used.
stiff
The global right-hand side vector F in (3.8) is prepared in the array aslod. For each element, the contribution to this is calculated by
(3.10), and the element stiffness matrix K^e by (3.9). The integrations
involved are performed numerically, by the one-point Gauss rule.
These local coordinates of the Gauss point are stored in (s,t), and
mapped to Cartesian coordinates (xg,yg) using the mapping (3.4-5).
The determinant of the Jacobian matrix J in (3.13) is found, and the
basis functions and Cartesian derivatives in (3.16-25) calculated at the
Gauss point. The entries in the element stiffness matrix estif are
formed according to (3.9), and written to the scratch file.
The contribution to the load vector by
(3.10) is added in also. For this, the right hand side vector f(x,y)
in (3.1) must be specified. This is done in the function rhsfn(x,y)
to be found in the file POISS.FOR after the result
subroutine. It is programmed as
f(x,y) = 1,
but you can change this
for your own problem.
assmb
The element stiffness matrices are read from the scratch file, and assembled
to form the global stiffness matrix K in the array gstif. As
K is a banded matrix, it is stored compactly, only the band itself being
stored. First the halfbandwidth nhbw is calculated;
this is the maximum, over all elements, of the maximum difference
between node numbers within the element, plus one. Then the lower diagonal
of the band only (since K is symmetric) is stored in gstif, in a
compact block with npoin rows and nhbw columns. The i,j'th
entry in K is stored at the location
K_ij -> gstif(i,j-i+nhbw+1) .
solve
This subroutine uses Choleski decomposition to solve (3.8). See Chapter 8 for
further details of the assembly and solution algorithms.
result
This section writes the results --- the nodal values of U --- into the
print file ( filename .prt). It also appends a one-line header
to the output file ( filename .out) and then for each node n
writes N and UN, the value of U at node n,
in the format I5, E10.3.
This file will be read by the post-processor.
Back to top
Program structure
After reading the mesh data from the input file, and echoing it to the
output and print files, the user is asked whether this is a plane stress, plane strain
or axisymmetric analysis.
The main program then performs the finite element calculation by calling the
following subroutines in turn:
- stiff
- This subroutine calculates the element stiffness matrices
and writes them to the scratch file.
- load
- This subroutine reads the loading data from the input file,
and calculates the equivalent nodal load vector f, stored in the array aslod.
- assmb
- This subroutine assembles the element stiffness matrices
to form the global stiffness matrix.
- solve
- This subroutine solves the global stiffness equation
(4.16) for the nodal displacements.
- result
- This subroutine uses the nodal displacement results to
calculate stresses at the Gauss points of each element. It then writes this
displacement and stress information to the print and output files.
Some more detail on each of these subroutines is now given. Much of the
coding is identical to the corresponding subroutines in the Poisson program,
and the description in Section 3.4 should be read first. Note that the
source code is provided with comment lines explaining the sections, and
describing the variables used.
stiff
Within the element loop, this subroutine first of all identifies which
material property
set that element is in, and then extracts its Young's modulus and
Poisson's ratio from PROPS.
To deal with the double integral in (4.24)
the matrix multiplications in the integrand are carried out in the do
loops 39,40, and 49, and the integral is approximated by the 2-by-2
Gaussian rule (4.25) in the do loop 75.
The element stiffness matrix is formed by using the following
subroutines:\
When the stiffness matrix ESTIF for each element has been formed, it is
written in a temporary scratch file.
load
This subroutine reads the loading data, and forms
the global load vector ASLOD. There are two types
of loading that we consider: point loads and surface tractions.
For point loads, data is read in from the statement
15 read(15,520) lcard,ld,np,fx,fy
where np is the global node number of the node where the
point load is applied, and fx and fy are the x and y components of the point load.
For surface tractions, there may be several different sets of surface
traction. Details of each set are read from the statements
25 read(15,530) lcard,iset,i1,i2,pn,pt
where iset is the surface traction set number,
and pn and pt
are the normal and tangential components of the traction. These tractions
are associated with loaded element edges by reading the data
30 read(15,540) lcard,nnas,l,n1,iset1,n3,iset3
where l is the element number, n1 and n2 are the
corner nodes defining the loaded edge of the element, and where the traction
load has the values in
surface traction sets iset1 and iset2 respectively.
assmb
The assembly subroutine operates as in the Poisson program, but with the
added complication that there are now two degrees of freedom (the x and
y displacements) at each node. The halfbandwith is thus
NHBW = 2*MDIFF + 1
where MDIFF is the maximum difference in node numbers within an element,
over all elements. The mapping of the entries in K into the compact storage
array gstif is also modified, for this reason.
solve
This subroutine solves the matrix equation K \vecu = \vec f
for \vec u, where:\
\qquad K is the global stiffness matrix GSTIF,
\qquad\vec u is the global displacement vector ASDIS, and
\qquad\vec f is the global load vector ASLOD.
The method of solution of this matrix equation is
provided by symmetric Choleski decomposition. See Chapter 8 for further details.
result
This subroutine appends the displacement and stress results
to the output file. In the case of axisymmetric analyses, there are four stress components
output; the last of these is the hoop stress sigma_theta.
Back to top
The PREFEL pre-processor allows you to use in your mesh one-dimensional
line elements containing two nodes or three nodes (at the two ends and the mid-point
of the line). In elasticity analyses, these may represent either bars (which
have axial stiffness) or beams (which have axial and flexural stiffness). The
advanced elasticity program ELADV.FOR interprets such elements as beams, although
the subroutines which evaluate the shape functions, etc. also contain the coding for
the simpler bar elements.
In FRAME.FOR, the theory of the standard 2-noded beam element, with Hermite cubic
interpolation functions for the displacements and rotation, is implemented. The
post-processor FELVUE also uses these interpolation functions in graphing beam elements.
The Euler beam theory for the transverse displacements is combined with the 1D linear
elasticity theory for axial displacements, producing an element with flexural and
axial rigidity (governed by the Moment of Inertia I and Young's modulus E).
The stiffness matrix and equivalent load vector (for both point loads and distributed loading)
are programmed explicitly in subroutines stiff and load.
Further, the beam element is related to the Cartesian xy-coordinate system by rotation
transformations, so that a multi-beam frame structure can be modelled.
There is no equivalent to the subroutine assmb. This is because in FRAME.FOR
an iterative solver has been used for solving the global matrix equation.
This is the conjugate gradient method with diagonal preconditioning.
It is described in Chapter 8 on global matrix storage and equation solution algorithms.
For this solver, there is no need to assemble the global stiffness matrix K; this matrix is
only used in forming matrix-vector multiplications Kv. This calculation in performed in
subroutine ematgv by reading each element matrix in turn from the scratch file where
they were written by stiff, multiplying each entry by the appropriate entry in v and
adding to the appropriate entry in the output vector.
Sample datafiles framexa1.dat, framexb1.dat are provided; these are documented in Chapter 9.
Three-noded beam elements
Engineers often desire to use beam elements in the same mesh as 2D elastic elements --
in problems involving soil-structure interaction, for example. If a two-noded beam element is
joined along an edge in a mesh of 8-noded quadrilateral elements -- to represent a beam buried
in soil, for example -- an unrealistic stress distribution results, because forces are only
transmitted between beam and soil at the beam nodes. A better, though not perfect, model can
be obtained if a 3-noded beam element is used. We now describe such an element, where the
(u,v,\theta) degrees of freedom exist at the two end nodes, and at the mid-node there is
only (u,v) displacement degrees of freedom. This is analogous to the composite element
used for soil consolidation in program CONSL.FOR (where the 3rd d.o.f. represents pore
pressure). Such an element was derived in the PhD thesis of Chao (see Chapter 10).
Five quartic interpolation functions can be derived which satisfy the continuity
requirements.
Explicit formulas for the stiffness matrix entries, and for converting distributed loads to
equivalent nodal loads can then be derived following the process of the previous section.
This element has not been programmed into FRAME.FOR, but it has been included in the
advanced elasticity program ELADV.FOR; both 2-noded and 3-noded beam elements are
included there, and the explicit form of the element stiffness matrices, after the rotation
transformation has been applied, are given in subroutine estifb. Subroutine load
shows the formulas for equivalent nodal loads (The 2-noded and 3-noded beam elements are
identified by NODSID=21 and NODSID=31 respectively).
There is a datafile eladbm1.dat which models a beam part-embedded in an elastic block,
with a normal surface traction applied on the element at its tip; see Chapter 9 for more details.
Back to top
The ELADV program is suitable for full-scale analyses, with up to 800 elements, 2000 nodes.
The full range of 2D elements may be used, namely linear and quadratic triangles and
quadrilaterals, the 6-noded quadrilateral, and the infinite side and corner elements.
Each element has its own shape functions and Gauss integration rule, and the coding for
these is included in subroutines sfr and gauss. An additional subroutine mapfun provides mapping functions needed in place of shape functions by the
infinite elements in some stages. The subroutines jacob and bmat are also
generalized. To identify the element type, the parameters nnode and nside
defining no. of nodes/sides, are passed into the subroutines, and in general the DO loops
which in ELAST ran from 1 to 8, now run from 1 to nnode, etc.
Subroutines sfr and gauss also contain coding for shape functions and Gauss
rules for the remaining elements available in PREFEL , namely linear and quadratic 1D
truss elements, and finite and infinite joints. Users wishing to develop main programs
using joint elements, should consult the specialist literature for the various models.
For further generality, the coding in ELADV does not assume two degrees of freedom at
each node, but reads the no. of d.o.f.s at node n from array ndof(npoin). To relate
each d.o.f. to the global d.o.f., it constructs an array nvar(npoin); the i'th
d.o.f. at node n will correspond to global component nvar(n)+i. Note that this is not
strictly needed in ELADV (though see next paragraph), but is useful for developing programs
for coupled analyses such as CONSL.FOR.
ELADV also contains coding for one-dimensional line elements. Linear and quadratic
(2-noded and 3-noded) bar or truss elements can be analysed in exactly the
same way as 2D elements, and the subroutines sfr, jacob contain the coding
for such elements (distinguished by NODSID = 21, 31 respectively). This coding
is not however used, since such element types are in fact identified in stiff
and load as cubic and quartic beam elements. The theory of these elements has
beend described in Chapter 5. Their stiffness matrces are formed explicitly, not by
numerical integration, and the coding for this is in subroutine estifb, called
from stiff. The equivalent nodal loads for distributed loading on these elements,
is also coded explicitly in load. Note that the end-nodes of these elements have
displacement and rotation degrees of freedom, so that ndof = 3.
For large-scale analyses, the Gaussian elimination solver solve is impractical
because it requires assembly of the global stiffness materix. Even with the band
storage used in assmb this can take up too much storage. ELADV therefore replaces
the assembly and solve subroutines with subroutine front, which uses the frontal
solution algorithm of Irons. This is a direct solver which essentially performs the
assembly and elimination processes simultaneously, only holding a so-called `front'
of variables currently active, in storage.
The amount of storage needed depends on the maximum frontwidth.
The parameters MFRON and MSTIF replace MHBW at the start of ELADV, as global storage
dimensions. The vector GSTIF of dimension MSTIF holds the part of the global matrix
within the current front (of max. dimension MFRON), and so MSTIF should be set to
MFRON*(MFRON+1)/2 for a symmetric matrix.
The subroutine writes information to scratch files on channels 19, 20 and 21 as it
proceeds through the forward decomposition, reading them back when performing
back-substitution. If a re-solution with a new right-hand-side vector is required, much
of the former process is unnecessary, so a parameter IRSOL is set to 2 if a re-solution
using an already-decomposed matrix is to be performed (e.g. with incremental loading);
IRSOL=0 for the first solution, as in ELADV. To avoid continually writing to the scratch
files, FRONT has a buffering feature. The size of the buffer is given by the parameter
MBUFA, set to 100 in ELADV; the user can adjust this for maximum efficiency. Note that
FRONT uses the array nvar for global d.o.f. positions, so is not restricted to meshes
with 2 d.o.f.s per node. This automatically copes with the beam elements, which have
3 degrees of freedom at their end-nodes.
FRONT reads two vectors of specified d.o.f.s (NFVAR lists the global d.o.f.s which are
fixed, and SPDIS lists the corresponding values) rather than the fixity codes in NKODE,
so this information is converted from NKODE in the main program after input.
In load, any specified nodal displacements are read, and this information is
added to NFVAr and SPDIS; the variable NSDIS tells the total number of such d.o.f.s.
A new parameter MFIXV sets the maximum number of fixed and specified d.o.f.s, dimensioning
NFVAR and SPDIS.
Body force loading is handled in load, if NBSET > 0, being converted to equivalent
nodal forces as with surface tractions.
As described in the pre-processor chapter, an in situ stress field can be prescribed.
This is useful in geomechanics applications, where the soil/rock mass is in a state
of stress before loading. The type of stress field is decided by ISTYP, read from
the first dataline in this section. The two parameters defining the stress field
are also read from this line, and are:
\begin{tabular}{l} ISTYP = Stress field type (1 or 2) \\ ISET = Set no. \\ REAL1 = horizontal stress \sigma_x (ISTYP=1) or soil density \gamma (ISTYP=2) \\ REAL2 = vertical stress \sigma_y (ISTYP=1) or lateral stress ratio K_0 (ISTYP=2) \\
\end{tabular}
Thus, for ISTYP=1 the stress at all Gauss points of soil elements (with material
type > 1) is (\sigma_x , \sigma_y) . For ISTYP=2, the stress state at a soil
element Gauss point with coordinates (x,y), y<0 is (-K_0 \gamma y , -\gamma y) .
Note that in ELADV, material type LMTYPE=1 is assumed to mean the element is a
structural element with no initial stresses; soil/rock elements have LMTYPE=2,3,....
Also, with an overburden stress field, the plane y=0 is taken as the ground level.
If no excavation boundary nodes are given in this section, ELADV assumes that the
equivalent nodal forces from this stress field should be applied in the load vector
ASLOD; you can use this to `turn on the gravity' to represent constructing an
embankment, for example. If an excavation boundary is defined, the equivalent nodal
forces are stored in TLOAD, and equal and opposite forces normal to the boundary
(representing unloading of the boundary) are added to ASLOD. Note that an initial
equilibrium stress and load state must be carefully stored if a nonlinear material
model is to be used.
Stresses at Gauss points are stored in an array strsg(3,4,melem); the i'th stress
component of the ij'th Gauss point of element l is stored in strsg(i,ij,l). In result, the stress increments are added to any initial stress state to give
the final stress state, for output. ELADV has, like ELAST, the facility for testing the minor
principal stress for yield, against a tensile strength input as the 4th material property.
This is described in Chapter 4, and the yield zones can be viewed in the post-processor.
Back to top
This program analyses plane thermoelasticity problems, in which we have the temperature
field \theta(x,y) as an additional, scalar degree of freedom.
The theory for this is summarised
in the Appendix, and results in the following augmented stiffness matrix equation:
{\samepage
/ T \ / \ / \
| K Q | | u | | f |
| | | | = | |
| 0 -H | | t | | r |
\ / \ / \ /
Here, K, u, f are the elastic global stiffness matrix, nodal displacement and
equivalent nodal load vectors, as described in Chapter 4; t is the temperature.
As in ELAST, both plane
stress and plane strain is allowed for, and decided by the user at run-time.
The matrix H results from the Poisson equation for temperature diffusion:
\ben
k_x \frac{\partial^2 \theta}{\partial x^2} + k_y \frac{\partial^2 \theta}{\partial y^2} = -r
\ee
which must be discretized over the mesh. The coefficients of thermal conductivity in the x
and y directions are k_x, k_y ; in THERM.FOR we assume isotropy: k_x=k_y=k . This is
exactly equivalent to the Poisson equation theory described in Chapter 3, and the form
of H is given by equation (3.9). We do not consider radiating boundaries in this program;
the temperature boundary conditions are either insulated (no normal flow, which is the
natural boundary condition and requires no special coding) or a Dirichlet boundary on
which the temperature is fixed. This is specified by fixing the temperature d.o.f.s for
these boundary nodes when preparing the mesh, and then inputting the Dirichlet boundary
plane data as described in Section 2.6.
The coupling of the two fields (displacement and temperature) occurs through the
constitutive equation
\ben
\vec{\sigma} = D \vec{\epsilon} - \beta \vec{m} \theta
\ee
where \vec{m} = (1,1,0)^T , and \beta is the coefficient of thermal expansion. This is
inserted in the equilibrium equation
\ben
B^T \vec{\sigma} + \vec{f} = 0
\ee
and introduces a coupling matrix Q . The form of Q is
\ben
Q = \int N^T M^T B dV
\ee
where
\[
M = {\left[\begin{array}{cc}
\beta & 0 \\
0 & \beta \\
0 & 0 \\
\end{array}\right]}
\]
In practice, the temperature degree of freedom is inserted after the displacement degrees
of freedom, for each node. We use the arrays ndof, nvar to relate the nodal and global
degrees of freedom, as described in ELADV above. THERM has coding for two types of element:
the linear 4-noded quadrilateral element ( NODSID=44) and the 8-noded serendipity
quadrilateral element ( NODSID=84). In both cases we have temperature degrees of
freedom at all nodes, so ndof(np)=3 for each node np in the mesh.
Subroutine stiff works element-by-element, and should be compared with the corresponding
subroutine in ELAST. It first notes the material properties:
- Young's modulus E
- Poisson ratio \nu
- Thickness T (default = 1.0)
- Thermal conductivity k
- Coefficient of thermal expansion \beta
and stores these as ym,pr,thick,conduc,expand respectively. \\
The elastic stiffness matrix is then formed in estif(nevab,nevab), just as was done
in ELAST (with the generalisation that the variable nnode (which may be 4 or 8) is used
in loops instead of assuming 8).
The other two matrices Q,H are then formed in qmatx, hmatx, according
to the formulae above. Once the process has been performed at each Gauss-point, and added
to perform the gaussian integration, the augmented element matrix is assembled from these
matrices, in xmatx. The rows and columns of this array are then shuffled, to place
the row/column corresponding to the temperature d.o.f. for node n , immediately after the
rows/columns for the displacement d.o.f.s for that node. It is xmatx which is
written to the scratch-file on channel 18.
Subroutine load is just as in ELAST, with the generalisation that surface tractions
over edges of linear elements use the linear 1D shape functions to calculate equivalent
nodal loads.
Subroutine bndry is lifted from POISS.FOR, and reads the Dirichlet boundary plane
data for the 3rd (temperature) nodal degree of freedom.
Subroutine assmb is very similar to that in ELAST, and uses the same mapping
of the band of K into the array gstif, namely:
\[
K_{i,j} \rightarrow \mbox{gstif(i,j-i+nhbw+1)}
\]
(where K is now the augmented element matrix, including the temperature rows/columns)
with one important difference: as the element stiffness matrices are non-symmetric, the whole
band -- above as well as below the diagonal -- must be stored. The array gstif
thus has nbw columns, where nbw = 2*nhbw+1.
Subroutine solve solves the linear system. As this is non-symmetric, we cannot use
the symmetric Choleski decomposition. Instead, we use the familiar Gaussian elimination
algorithm. The subroutines called in succession are:
- fixdof
- to add a large value to the diagonal coefficient for fixed degrees of
freedom. Temperature d.o.f.s on Dirichlet boundaries are dealt with as in POISS.FOR;
- reduce
- to perform the reduction of K to upper triangular form, with
simultaneous elementary row operations on the r.h.s. vector aslod;
- backsu
- to perform the back-substitution, the solution appearing in asdis;
- reactn
- to back-calculate the reactions at fixed displacement d.o.f.s.
The results are then output: nodal displacements and temperatures, and Gauss-point stresses
and strains.
The example datafiles for use with THERM, described in Chapter 9, model a cooling fin
joined to a hot boiler at its base. From these results it becomes plain that the element
formulations used here are prone to instability. This is particularly evident with the
linear quadrialteral element, used in the file thermln2.out, where the
famous `hourglass' instability occurs. The mesh of quadratic elements also exhibits
some instability of results, albeit much less. The solution to this problem, is to use an
element formulation in which the extra degree of freedom -- temperature, in this case --
is given a bilinear discretization, which the displacements are modelled using the 8-noded
serendipity formulation. The resulting 8-noded quadrialteral element would thus have
(u,v,\theta) degrees of freedom at its four corner nodes, and (u,v) degrees of
freedom at the midside nodes. The programming of such an element is demonstrated in the
soil consolidation program CONSL.FOR, described next.
Back to top
This program is designed to model soil-structure interaction, with an elastic structure
(elements with lmtyp(l)=1) in contact with a saturated soil/rock mass (elements with lmtyp(l)=2).
The program performs a time-dependent analysis, starting with an initial solution for the
undrained deformations immediately the load is applied, and then stepping forward in time.
The fluid flows through the soil mass, away from the loaded area, and further elastic deformation
occurs. After a sufficiently long time there is no further appreciable deformation/flow,
and the drained solution has been reached.
The CONSL program uses the eight-noded quadrilateral element from ELAST, but now there are
additional degrees of freedom for the pore-pressures at the corner nodes of LMTYP=2
elements. The material properties for soil elements are:
- Young's modulus E
- Poisson ratio \nu
- Horizontal permeability k_x
- Vertical permeability k_y
- Effective porosity \eta / K_f , where K_f is the fluid bulk modulus.
The program is restricted to plane strain (as plane stress makes little sense in this
application), so the material thickness defaults to 1.0.
Boundary conditions for displacements are handled using nodal fixities, as in ELAST. This
is also used for the pore pressures. Note that p is properly the excess pore
pressure, so we do not need to input an initial pore pressure distribution, and at
permeable boundaries we have p=0 ; the natural boundary condition is for an
impermeable boundary.
There are now four blocks within the augmented element stiffness matrices (see the
theory in the Appendix), as was the case with THERM.FOR. All the coding is generalised to allow for variable degrees
of freedom per node; the array nvar described in ELADV above is used. In particular,
the assmb and solve subroutines are generalised to use nvar.
The theory of elastic consolidation is summarized in the Appendix. The only difference between
the notation there and that implemented in CONSL, is that the time discretization parameter
\theta is used in CONSL with the sense
\ben
u(t) \doteq (1-\theta) u(t_0) + \theta u(t_1)
\ee
where t_1 = t_0 + \Delta t . Thus, theta should be replaced by 1-\theta in the theory
in the Appendix.
Suppose that we have the nodal displacements and pore pressures at time t_0 stored in the
vector u_0 , and we wish to calculate the displacements and pore pressures at time
t_1=t_0+\Delta t , in vector u_1 . The theory tells us that:
\ben
(\bar{X} + Y ) u_1 = \bar{\theta}f_n + f_{n+1} - (\bar{\theta}\bar{X}-Y)u_0
\ee
where \bar{\theta} = (1-\theta)/\theta and
/ T \ Form K in estif(16,16)
| K -Q | Form Q in qmatx(4,16)
X = | | Form H in hmatx(4,4)
| 0 -H | Assemble whole matrix in xmatx(20,20)
\ / Write xmatx to file
/ \
| 0 0 | Form Q in qmatx(4,16)
Y = | | Form S in smatx(4,4)
| -Q -S | Assemble whole matrix in ymatx(20,20)
\ / Write ymatx to file
\bar{X} is formed from X by multiplying the entries of the submatrix H by
\theta \Delta t .
The coupling matrix Q and diffusion matrix H correspond to those in THERM; the
porosity matrix S has a similar structure to H and is of minor importance. The
programming of subroutine stiff is a straightforward extension of that in
THERM; in this case we write the the augmented matrices xmatx,ymatx to file
for each element in turn. X , in xmatx, is turned into \bar{X} in the
assembly subroutine; this is done to avoid having to repeat stiff each time
the timestep \Delta t is changed. The global right-hand-side of the timestepping
equation above, is assembled in the subroutine newrhs.
In CONSL.FOR we have introduced skyline storage, more efficient than band storage --
see Chapter 8 for details. The subroutine mask constructs the pointer array
needed.
For the equation solution, note that the matrix \bar{X}+Y is symmetric, but is not
positive definite because the diagonal entries of H are negative. Thus the symmetric
Choleski decomposition cannot be used; instead we use the decomposition L.D.L^T
where L is lower triangular with 1's down the diagonal, and D is diagonal. The
decomposition and substitution subroutines are ldldec and ldlsub; see
Chapter 8 for full details of the coding.
As mentioned in the previous program THERM, the standard quadrilateral element
for this application is an 8-noded serendipity element in which the pore
pressure variable is interpolated linearly, using d.o.f.s which exist only at
the four corner nodes. A further complication, is that we wish to allow elements
representing a structure (in which there are no pore pressure d.o.f.s) to be used
in the mesh. Then corner nodes on an interface between a soil element and a
structure element will have (u,v,p) degrees of freedom, but the third d.o.f. will
be ignored when constructing the stiffness matrix for the structure element. This is
coped with by checking the material type of the element, in stiff; if lmtyp(l) = 1 it is a structural element, and we use ndof = 2.
A suitable timestepping r\'egime for use in CONSL is as follows. (The various
parameters defining it are read in from the final line of the data-file --
see section 2.10 above). A special timestep dtone is used for the first
solution (starting from u(0) = 0 ). With the dual-level element used in CONSL,
it is possible to take dtone = 0.0. Timestepping begins with the initial
timestep dtint; five steps are taken using this timestep. The timestep is then
multiplied by the factor ftime, and used for a further five timesteps. At
the end of each set of 5 steps, the user is prompted to say whether s/he wishes to
continue timestepping. If ftime is taken as 2.0, this regime will give
roughly equally-spaced points along the time axis on a logarithmic plot of
displacement against time, as is often used. Moreover, the stiffness matrix
is constant if the timestep \Delta t is unchanged, so it only needs to be
decomposed in the first of each set of 5 timesteps. (The user will observe that
after the first step, the subsequent steps are performed much more quickly.) There is
therefore no solve subroutine; the routines ldldec and ldlsub
are called directly from the main program. Timestepping
stops when the user chooses, or when the input maximum time ttend is reached.
Results are appended to the .out file after each set of five timesteps, and when
the program ends. These result sets can be viewed consecutively using FELVUE.
To enable the user to know how consolidation is progressing, the user is invited at run-time
to input the node numbers of five nodes at which displacements and pore pressures will be
written to screen at each step. The user should choose these nodes using the `Pick node
by mouse' option in the pre-processor, before running CONSL. Suitable nodes are
those lying on the surface under loaded areas, or down the centreline of a load. See
the example in Chapter 9.
Back to top
The theory of Mohr-Coulomb elasto-plasticity is provided in the Appendix. For
elasto-plastic material elements, five property values must be prescribed: in addition
to the Young's modulus and Poisson ratio, the third component is the triaxial stress
factor k and the fourth is the yield strength \sigma_c (denoted by triax
and sigc). (The relation of these
parameters to cohesion c and angle of friction \phi is given in the Appendix).
The PLAST program assumes associated flow, so that the resulting elstoplastic
stiffness matrix is symmetric and can be solved by the same routines as in ELAST,
namely band storage and Choleski decomposition. The subroutine mask, introduced
in CONSL.FOR, is however used to determine the halfbandwidth nhbw for the mesh.
See below for the storage and solution routines used in PLADV.FOR. The 8-noded serendipity
quadrialteral element is used, as in ELAST.
Should PLAST can be generalised to non-associated flow, a solver for non-symmetric
systems (such as the Gaussian elimination used in THERM.FOR, or the AFRONT subroutine in
VPLAS should be used.
The load must be applied incrementally, and within each increment a stress relaxation
iteration seeks equilibrium. The load increment loop terminates at statement 200 in
the main program. Within this, the iteration for equilibrium starts at statement 100.
Subroutine check evaluates the residual norm after each iteration, and if this
is less than gtoler (set in the main program at 10^{-4} ) the iteration is
taken as having converged. For greater flexibility, the user may wish to input gtoler
and other tolerances used, as extra parameters in the datafile, when preparing it
with the pre-processor (see \S2.9, block 2).
A new subroutine stress checks for yield, and
updates the stresses; a new subroutine residu calculates the residual or
out-of-balance forces. Yield at a Gauss point ij of element l is denoted by setting
mohrg(ij,l) to 1; if mohrg(ij,l)=0, the material is still behaving elastically.
In the program, material type lmtyp(l)=1 is reserved for purely elastic material, and
the plasticity parameters and yield criterion are only used for elements l
where lmtyp(l)=2,3,...
The Gauss point stresses are stored in strsg as in ELADV, and separate arrays keep
track of the accumulated load applied so far ( accld), and the total load to be applied
( totld), in addition to the current incremental load in aslod.
The yield criterion used is the Mohr-Coulomb criterion
\ben
F(\vec{\sigma}) \equiv \sigma_1 - k\sigma_3 - \sigma_c > 0.
\ee
The yield function F is evaluated in subroutine ycrit. The flow vector
\vec{a} = \partial F/ \partial \sigma are found in flowv, and the matrix of
second derivatives in flowm. These are stored in avect and dbds respectively.
If an alternative yield function was to be used, these subroutines would need to
be changed.
To avoid complicating the plasticity program code, a separate `main engine' PLADV has been
developed containing alternative methods of storage and equation solution. A parameter ISOLA is input at run-time to decide the solver used, as follows: \\
- ISOLA = 0 for direct solution using Choleski ( L.L^T ) decomposition and skyline storage;
- ISOLA = 1 for direct solution using L.D.L^T decomposition and skyline storage;
- ISOLA = 2 for iterative solution using conjugate gradients with diagonal preconditioning;
- ISOLA = 3 for Incomplete Choleski preconditioning with conjugate gradient solution.
These solvers, and the additional parameters requested for the iterative solvers, are
described in more detail in Chapter 8. PLADV uses skyline rather than band storage. There is
a trade-off between space efficiency and CPU time required; PLADV with ISOLA=0 runs more
slowly than PLAST.
Back to top
It will be found that the elasto-plasticity program PLAST is not very robust for
large-size problems. This is due to the sensitivity of the equilibrium iteration, and the
restriction to associated flow. A much more robust program
results from using viscoplasticity theory, where the stress-equilibrium iteration is replaced
by a time-stepping algorithm. It is also desirable to allow non-associated flow,
in which an angle of dilation is included. Such a program is given in VPLAS.FOR, which
also includes a sophisticated inner
iteration which attempts to reach local stress equilibrium at yielded Gauss points. In line
with the recommendation that VPLAS be used for solving practical plasticity problems,
a frontal solver has been included, which is very space-efficient. The symmetric solver front has already been described in ELADV.FOR above.
See the Appendix for Mohr-Coulomb elasto-viscoplasticity theory. The VPLAS program is
developed from PLAST, and the variables described in the previous section apply, with some
additions. Thus, the material property sets need, in addition to E and \nu , the
triaxial stress factor k and yield strength \sigma_c as in PLAST.
The fifth component is the flow rate \gamma , denoted gamma.
The VPLAS program has been generalised to allow for non-associated flow. For this,
the sixth component is the dilation factor \alpha (denoted by dilat in stiff),
related to the angle of dilation \psi . If \psi = \phi (so that k = \alpha )
the flow rule is associated, and the resulting stiffness
matrix is symmetric; the frontal solver FRONT is used. If 0 <= \psi < \phi , the
flow is non-associated, and the stiffness matrix becomes nonsymmetric. In this case,
a frontal solver adapted for nonsymmetric matrices is needed. This is provided in
the subroutine afront. Now the relation between the dimensioning parameters MFRON
and MSTIF is MSTIF=MFRON*MFRON, since the full matrix of active variables must be
stored. Note that in either front or afront, some stored data can be used
in re-solutions, even if the stiffness matrix has changed, since its structure is still
the same. This is exploited by setting IRSOL=1 in re-solutions. Further details of
the frontal solvers is given in Chapter 8.
The viscoplastic D matrix is calculated in the subroutine dmatvp. The H matrix
involved in the theory is calculated in subroutine hmat. Subroutine strain
tests for yield, and calculates the viscoplastic strain velocities (stored in the array vivel). Subroutine steady decides if these are small enough throughout the mesh
to say that equilibrium has been reached.
The example datafiles