All the above programs, with one exception, need to be linked to the file FELSB.FOR containing input/output and other common subroutines. The exception is ELADV.FOR, which is completely self-standing. These programs should all be compiled in double precision (use the /dreal option in the Salford compiler). Remeber to compile FELSB.FOR with this option also -- unpredictable errors at run-time will be generated if a `main engine' compiled in double precision calls a subroutine in FELSB.FOR which is compiled in single precision, as the double-precision arguments of the call will not be mapped properly onto the subroutine's arguments.
The Appendix to the full manual summarises the finite element theory used in each program; a full description is beyond the scope of the manual, and the user is referred to the standard texts such as Zienkiewicz and Taylor (see Chapter 10). In this chapter some programming and data input aspects of the `main engines' will be discussed. Note that the different algorithms for global matrix storage and solution, are described in the following Chapter and will not be covered below.
When POISS.EXE is run, the program will prompt you to enter the name of the .dat data file -- which should be entered without the .dat extension -- and will open this on unit 15 and read it using the subroutine input. It will write the input data in readable form, to a file with a .prt filename extension, opened on unit 16, using the subroutine print. Also, the data is echoed to an output file with .out filename extension, opened on unit 17. During the finite element calculation, a scratch file is opened on unit 18 to store the element stiffness matrices. After performing the finite element calculation, the results are appended to the print file, and the output file, by the subroutine result. The .prt file can be printed, and the .out file is used by the post-processor FELVUE.
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.
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:\
It also forms their local derivatives with respect to \xi and \eta, and returns these in the array DERIV.
For point loads, data is read in from the statement
15 read(15,520) lcard,ld,np,fx,fywhere 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,ptwhere 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,iset3where 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.
NHBW = 2*MDIFF + 1where 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.
The method of solution of this matrix equation is provided by symmetric Choleski decomposition. See Chapter 8 for further details.
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.
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.
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.
/ 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:
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:
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.
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:
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.
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: \\
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.