Coupled Problems - Programming the Finite Element Method - page 419

store kc stores element kc matrices
val

nodal loads weighting factors

related to the volumetric strain

vol

used to compute coupling matrix

volf

The analysis of the behaviour of porous elastic solids (Biot 1941) under load is in
many ways analogous to the coupled flow analysis described in the previous program. The
displacements of the soil skeleton take over the role of the velocities u and v , and the
excess porewater pressure (now called u w ) the role of the fluid pressure p (nodal freedoms
are in the order u , v , u w ).
The differential equations to be solved are (2.136) and (2.137). Due to the coupling of
fluid and solid phases there arises the complication that the applied “total” stresses, { σ } ,

are divided between a portion carried by the soil skeleton, called effective stress σ σ σ ,and

a portion carried by the pore water, called in soil mechanics the pore pressure and denoted
in Chapter 2 by

.
After discretisation in space by finite elements, the coupled equations are given by
(2.139). These can be seen to be partly algebraic equations and partly first order differential
equations in time. In the incremental load method used here, discretisation in time by the
θ -method leads to equations (3.115), which are in principle no different to (3.94) for
uncoupled problems. If using an assembly approach, solutions will involve setting up the
coupled global “stiffness” matrix on the left side of these equations ( kv ), followed by an
equation solution for every time step to obtain the incremental solutions followed by an
update of the variables from (3.116). For constant element properties and time step t ,the
left hand side matrix kv needs to be factorised only once, the remainder of the solution
involving matrix-by-vector multiplication on the right hand side, followed by forward and
back substitution.
Closer examination of equation (3.115) will reveal that some of the diagonal terms
of the left hand side matrix will be negative, thus the usual Cholesky solution strat-
egy will fail due to the need to take square roots. The subroutines sparin gauss and
spabac gauss which operate on the skyline kv vector use Gaussian [ L ][ D ][ L ] T factori-
sation are therefore introduced for the first time (see Section 3.8). Other new subroutines
include, load function , which reads the input load-time function and linearly interpo-
lates the function to produce a load-time function at the calculation time-step resolution
held in al ,and formke which forms the lhs element matrix ke from equation (3.115).
The structure chart describing the algorithm is shown in Figure 9.7.
The problem chosen is of a plane strain “oedometer” specimen as shown by the mesh
and input data given in Figure 9.8. The base and sides of the mesh are impermeable “no-
flow” boundaries, and “smooth” roller boundary conditions are imposed on the sides. The
top of the specimen is drained and subjected to a “ramp” loading of the form shown in
Figure 9.9.
The first line of data reads the number of elements in each direction of the rectangular
mesh ( nxe and nye ), and the number of property types ( np types ). The properties are
read next in the order k x /γ w , k y /γ w , E and ν .Since np types=1 in this homogeneous
example, the etype data is not needed. Next comes the rectangular mesh coordinate data
x coords and y coords , followed by the time stepping and output data. In this case,
the data calls for nstep=300 calculation steps, with a time step of dtim=0.01 .The

{ u w }