# Sample datafiles for POISS.FOR: solving Poisson's equation over an H-shaped region

File pohex1.dat contains Part A mesh details for a mesh of three linear quadrilaterals. The coefficients a_x and a_y are set at 1.0 in the material properties section.

In file pohex2.dat the nodes on the left-hand vertical edge have been given Dirichlet boundary fixities; the mesh has then been refined, and doubled by reflection to the right and upwards, and finally triangulated, to create an H-shaped region of width 4.0 and heaight 3.0, containing 384 linear triangles, with Dirichlet boundaries on the outer vertical sides. Part B data has been added, fixing the variable U at 0.0 on the left-hand Dirichlet boundary, and at 60.0 on the right-hand boundary. The other edges act as insulated boundaries.

The program POISS.FOR solves the quasi-harmonic equation, using the coefficients a_x, a_y specified in the data file and with the right-hand-side function equal to 1.0. The result of running POISS.FOR with pohex2.dat as input, is provided in the output file pohex2.out, which can be viewed using FELVUE.

The file pohex3.dat is produced from pohex2.dat by adding a radiating boundary on the lower edge of the crossbar of the H shape, with parameters QBAR=0.0 and ALPHA=4.0. THe effect of this can be seen in the contour plot of U, and the flowlines, in pohex3.out. An example mesh involving fewer than 50 nodes and 50 elements, which can be run in the demonstration version of POISS.EXE, is provided in podemo.dat. This models only the lower half of the H-shaped region in pohex2.dat, with a coarser mesh. Because of the top boundary being an insulated boundary, it is the same problem which is being solved, as with the full H-shaped mesh.

# Sample datafiles for ELAST.FOR:

## Thick cylinder, compressed on outer rim

File elcyl1.dat defines a coarse mesh comprising just one 8-noded element, representing one quarter of a thick cylinder: inner radius 1.0 and outer radius 10.0. The coordinates of the midside nodes of the curved edges are (0.707,0.707) and (7.071,7.071); this is able to model a true quarter-circle to within 1 E=200,000 KN/m^2 and \nu = 0.3 . The nodes on the vertical and horizontal edges are constrained in the x and y directions respectively, to provide by symmetry a solution for a complete circular annulus.

In elcyl4.dat the mesh has been refined to 100 elements, and a compressive loading applied to the outer edges of the outermost row of elements. The top element (adjoining the vertical boundary) has an applied pressure of 35 KN/m^2 (traction set 2); the neighbouring element has a ramp' loading which drops from 35 KN/m^2 to 30 KN/m^2 (traction set 1), and then the next elements in have a uniform pressure of 30 KN/m^2 . The ramp loading up to 35 KN/m^2 is reproduced at the horizontal boundary. ELAST.EXE should be run with this datafile, using the plane strain option at run-time. The resulting displacements and stresses can be viewed in FELVUE in file elcyl4.out.

The data-file eldemo.dat is a coarse mesh (just 12 elements) for a thick cylinder, inner radius 3.25, outer radius 10.0, elastic parameters E=20,000 and Poisson Ratio = 0.3, subject to a uniform load of 10 KN/m^2 on the outer rim. This can be run with the demonstration version of ELAST.EXE, again in plane strain. If the resulting stress contours are plotted in FELVUE there will be seen a small non-uniformity around the cylinder due to the coarseness of the discretization. Remember that FELIPE does not apply any stress-smoothing, and contours are plotted element-by-element, so that inaccuracies caused by an insufficiently fine discretization are not disguised.
The file eldemaxi.dat is a model of the same problem using an axisymmetric mesh; if this is run with axisymmetric analysis (choice 3) identical displacements of the inner and outer circumference to those in eldemo will be obtained.

## Plane stress example: bracket under tension

File elbrac1.dat shows an initial mesh for analysing the problem of a plastic bracket, 15cm long and 6cm wide, with a hole of diameter 3cm drilled through it, centred 3cm from the right-hand end. The bracket is 1cm thick, and its analysis is a case of plane stress. The left-hand end is fixed to the wall. A metal pipe is inserted through the hole (so that the bracket holds it to the wall), and we wish to examine the stresses and displacements caused in the bracket if the pipe attempts to push away from the wall. Our mesh only models half the problem: there is reflection about the centreline y=0 . The refined mesh (using the Automatic refine option, then further refining around the pipe) is contained in elbrac2.dat, which shoulf be input into ELAST.EXE, run using the plane stress option.

We model only the left half of the pipe cross-section, since in practice the right half will not exert any influence on the bracket; it will become detached from the plastic, which we cannot cope with in our program. We apply a uniform load of 10^5 KN/m^2 (equivalent to 10KN per square centimetre, as we are working in cm), pushing leftwards on the centreline of the pipe. The material properties used are:
For the pipe (material set 2): E=10^5, \nu=0.4
For the bracket (material set 1): E=2,500, \nu=0.2
The material thickness is set at 1.0 in the 3rd component of each set. A tensile strength of 10 is also specified for the bracket, in the 4th component. The material type (LMTYP) facility is also used: the bracket and pipe elements are assigned material types 1 and 2 respectively.

The resulting deformations and stresses of the bracket are best seen in FELVUE if you choose to plot only those elements with material type 1. The stress concentration is seen to be around the rear of the pipe section (use the Zoom facility to examine these in more detail), and if the Yield Zones'' plotting option (choosing yield code 1) is used, you will see the Gauss points at which the tensile strength has been exceeded. This is largely due to the plastic flowing in behind the pipe, which is not realistic: a better mesh would include a wing' of pipe extending part-way around the rear of the half-pipe section.

Note that this is a hypothetical example, and the material property values do not bear resemblance to real materials.

# Sample datafiles for FRAME.FOR

For FRAME.FOR, two simple examples of plane frames are given. File framexa1.dat describes a two-beam frame which is pin-jointed at either end, with a uniform load distributed along the horizontal beam (Note the negative sign for this, according to the way in which the compression-positive sign convention applies to 1D elements -- see Chapter 5). The dimensions are actually in inches, and the load is 1000 lbs/ft. The material properties of the beams are:
E = 3\times10^7 , I = 100, A = 10
This example is actually worked example 4.14 in the Schaum series text (Buchanan 1994), and the solution for nodal displacements and rotations can be checked against this. The conjugate gradient solver converges in 5 iterations, using diagonal preconditioning. In viewing the results in FELVUE, use the Plot 1D elts'' option and answer y' to the question whether the 1D elements represent beams. Note that the View displts'' option only displays 2D elements, so a blank screen will result from this option. The Plot d.o.f.s'' option will also not apply for frame analyses.

The second example, in framexb1.dat, is also based on a problem in Buchanan(1994), and Imperial units (inches, p.s.i.) are used. It consists of a three-beam frame, pin-jointed at its right-hand end and free to roll in the horizontal plane at the other, loaded with a uniform pressure on its middle member, and with a point-load acting downwards at the mid-point of the left-hand beam. The PCG iteration converges (to a tolerance of 10^{-5} ) after 19 iterations. The results can be viewed in FELVUE in the same way as with the previous example.

The datafiles provided for ELAST can also be used with the Advanced-level program ELADV. The following additional datafiles demonstrate some of the extra capabilities of ELADV:

## Beam part-embedded in elastic block

Following on from the beam examples above, datafile eladbm1.dat demonstrates that 1D beam elements can be used in the same mesh as 2D plane elements in ELADV. The mesh shows a vertical beam with its base embedded in an elastic block which is anchored at its base. It is loaded with a distributed force applied normally on the top element. The beam elements used here are the three-noded quartic elements, which are more closely compatible with the serendipity quadrilaterals than the standard 2-noded cubic elements.

In FELVUE, the deformed beam, together with the deformed boundary of the block, can be seen using the View 1D elts'' option with an exaggeration factor of 1.0. The Displacements'' option displays the deformation of the block itself, without the 1D elements. A mesh with more refinement would produce better results for stress contours, although the plot of Gauss point stress crosses is reasonable; principal stresses drawn in blue are tensile.

## Footing on infinite soil mass

The file eladft3.dat shows a very coarse mesh for a concrete footing on a semi-infinite soil mass. It is included to demonstrate how different element types can be combined in the same mesh. As well as serendipity quadrilaterals, there are 6-noded triangles and infinite elements (including a doubly-infinite element at the mesh corner). The load is applied as a normal surface traction on the top surface of the concrete footing.

Again, a finer discretization would produce more meaningful stress contours.

## Cavern excavated in pre-stressed rock mass

The example illustrating some of the extra advanced elasticity facilities in ELADV.FOR is of a cavern excavated in a pre-stressed rock mass of infinite extent. Half of the cavern is modelled in file eladcv1.dat, consisting of a curved roof (modelled by a curved Serendipity element as in the thick-cylinder example) above a 2m deep-by-1m wide rectangular space. Beyond a layer of serendipity elements surrounding the cavern boundary lie mapped infinite elements (including a 3-noded corner (doubly-infinite) element) representing the infinite rock mass. Elastic parameters are E=5,000 KN/m^2 and \nu = 0.2 . Because infinite elements are used all around the mesh, there is no need to fix any nodes in the y -direction; only nodes on the centreline are constrained in the x direction, for symmetry. This is a plane strain analysis.

In eladcv2.dat the mesh has been refined, and the excavation boundary around the cavern edge defined. A uniform stress field of \sigma_x = 20, \sigma_y = 40 is stipulated for all the elements, which have been assigned material type 2 as required in ELADV. A further feature is that a small 6-noded triangular element has been added after refinement to cut off the corner'' at the bottom of the cavern. The add elements'' feature in modify mesh'' was used, and by leaving blank the coordinates of the new node (on the midside of the hypotenuse), these were linearly interpolated by the pre-processor. The node was then moved using the change coordinates'' option, to produce a rounded profile; this greatly improves the accuracy of the stress field around the corner singularity.

The deformations and stress field resulting from the unloading around the cavern boundary, can be seen by viewing eladcv2.out. A realistic swelling of the cavern roof and floor, where areas of greatest deviator stress occur, is noted. (The initial exaggeration factor is too large, and a factor of 20 should be tried, in conjuction with the Zoom facility, to see the deformation around the cavern wall). If running ELADV with this example, notice that the equation solution is performed by the frontal algorithm FRONT.

This datafile can also be run as an axisymmetric analysis, to model a cylindrical, dome-roofed cavern. The hoop stress contours can be plotted as sigma_z, in this case.

# PLAST.FOR, VPLAS.FOR, PLADV.FOR: thick-cylinder example revisited

The file plcyl4.dat, for use with PLAST.FOR, is derived from elcyl1.dat, but written at Advanced level with additional material properties: triaxial stress factor k=3.0 and yield strength \sigma_c = 20.0 defining an elasto-plastic Mohr-Coulomb rock (stored as material property components 3,4). The associated-flow plasticity algorithm often requires very small load increments in order to achieve convergence, when the applied loads are asymmetrical, so for this example a uniform compressive load of 25 KN/m^2 over the outer circumference is applied. This load is applied in 5 increments. Viewing the output in plcyl4.out, you can see the drop in major stress on the innermost two rows of Gauss points, which have yielded plastically; the yield zones'' plot (choosing yield code 1) will show this zone. Remember that an associated flow rule is used, which keeps the elasto-plastic stiffness matrix symmetric.

The file vpcyl4.dat is an extension of elcyl4.dat for the viscoplasticity program VPLAS.EXE. Here, in addition to the material parameters k and \sigma_c given above for PLAST, the flow parameter \gamma = 0.01 , and dilation l are defined in material property components 5,6; in this case l=1.0 , i.e. a non-associated zero-dilation flow rule. VPLAS.FOR uses the frontal algorithm subroutine AFRONT, adapted for nonsymmetric matrices, in this solution. The loading is similar to that in elcyl4.dat, with 25 KN/m^2 around the outer boundary, rising to 35 KN/m^2 at each end; the viscoplasticity algorithm involving timestepping is much more stable than the plasticity algorithm involving iteration, especially if a non-associated flow rule is used, and it converges happily even with the non-uniform loading. In addition to an incremental loading regime with 5 load increments, a timestepping regime is defined in the last line of the data set. The Crank-Nicholson time discretisation parameter \theta = 0.5 has been used. Compare the stresses, displacements and yield zones output in vpcyl4.out, with those from the elasticity and elasto-plasticity analyses.

The advanced plasticity program PLADV.FOR can also be run with the datafile plcyl4.dat. Results will be identical to those from PLAST.FOR, but can be obtained using a selection of different solvers. If the preconditioned conjugate gradient solver (ISOLA=2) with diagonal preconditioning is used, each solution takes about 180 PCG iterations (for a convergence tolerance of 10^{-6} ). If the same convergence tolerance is used with the Incomplete Choleski preconditioner (ISOLA=3), the performance depends on the amount of fill-in which occurs. this is controlled by the fill-in factor which is also specified at run-time. if a factor of 0.1 is chosen, none of the 47,052 holes' within the envelope of K are filled in, and convergence takes 56 iterations. If a factor of 10^{-4} is chosen, then about 2,500 holes' get filled in, and the iteration speeds up, needing 24 iterations to convergence. With a fill-in factor of 10^{-6} , about 50% of the holes' (some 25,000) are filled-in, and only 5 iterations are needed to converge.

# THERM.FOR: cooling fin

To illustrate the plane thermoelasticity application, we consider the problem of a rectangular cooling fin, 0.5m wide and 1m tall, anchored at its base to a hot vessel, temperature 100 degrees, while its other three sides are exposed to the air, temperature 20 degrees. The elastic constants are E=10^6, \nu=0.4 , and this is a plane stress analysis with a thickness of 0.5cm. The coefficients of thermal conductivity and of thermal expansion, entered as material properties 4 and 5, are 1.0 and 5.0 respectively (these were chosen arbitrarily).

The fin is modelled with meshes of 4-noded linear and 8-noded quadratic quadrilateral elements, in datafiles thermln2.dat and thermqd2.dat respectively. Examining the displacements in FELVUE, you will see how the fin has expanded under this thermal loading (choose an exaggeration factor of 1,000 to get a proper view), and can also view the distribution of tensile stresses generated in the metal. Proceeding to the d.o.f. plots, a plot of d.o.f. 3 will display the temperature field in the fin.

As mentioned in Chapter 7, the results using quadratic elements are good, but in the results from the mesh of linear elements the well-known hourglass' instability is very evident. To avoid this instability, a two-level discretization should be used, with quadratic discretization of displacement d.o.f.s but only a bilinear interpolation for the temperature using values at corner nodes. This type of element is implemented for the soil consolidation application CONSL.FOR.

# CONSL.FOR: footing on saturated soil

In this problem, a concrete block (one element, material type 1) rests on a rectangular mesh of elements (material type 2) representing a soil mass. Elastic parameters are E=10,000 KN/m^2, \nu = 0.4 for the concrete block, and E=100 KN/m^2, \nu = 0.2 for the soil. In addition, the soil material set has component 3 set to 0.001, representing the porosity of the interstitial pore fluid. The horizontal and vertical soil permeability coefficients are set at 0.004.

In the refined mesh in file conslex3.dat the modify mesh'' option has been used to write at Advanced level and to add an extra degree of freedom at corner nodes only, for elements with material type 2. Note that corner nodes on the interface between block and soil elements, have the third d.o.f., but the CONSL program ignores this when processing the stiffness matrix for the block element (which has material type 1). This extra d.o.f. -- the pore pressure -- has been fixed at zero on the free surface y=0 , apart from at nodes underneath the impermeable concrete block; the fixity code has been used for this. A uniform normal surface traction of 10 KN/m^2 is applied to the top of the concrete block.

The timestepping r\'egime involves making an initial solution with timestep dtone=0.0, the results of which are written to file. Then fives timesteps will be made with a timestep of dtint=1.0, after which the timestep will be doubled ( ftime=2.0) and a further five steps made, and so on. Timestepping will halt when ttend=100.0 is reached, unless the user has chosen to halt earlier. A fully implicit timestepping algorithm is used ( theta=1.0).

In order to judge how consolidation is progressing, the user can enter at run-time the numbers of five nodes at which displacement/pressure results will be written to screen at each timestep. The nodes to be used should be chosen by viewing the mesh of conslex3.dat in PREFEL and using the Pick node by mouse'' facility (You will need to Zoom first!). In the present mesh, suggested nodes are:
Node 151 (on the centreline, on top surface of block)
Node 145 (on the centreline, on interface between block and soil)
Node 185 (on surface at one corner of the block)
Node 139 (on centreline, soil depth of 0.1m)
Node 137 (on centreline, soil depth of 0.3m)

If CONSL.EXE is run with datafile conslex3.dat, it will be seen that after 15 timesteps, when t=35.0 , consolidation is virtually complete, and the timestepping may be ended. The program appends the results after each set of 5 timesteps, to the file conslex3.out. Running FELVUE with this file, first the undrained deformation is displayed (an exaggeration factor of 10 for the displacements is appropriate. Note that the soil elements close to the footing bulge unrealistically -- a finer refinement in this part of the mesh is needed). You can see the undrained pore pressure distribution in the soil, by proceeding to the d.o.f. plots, and plotting d.o.f. 3 for material type 2 (soil elements). Remember to answer `y' to the question, if this d.o.f. exists only at corner nodes. Again, it is necessary to zoom to see the excess pore pressures under the footing, in detail. (A consequence of the bulging is that a zone of negative pore pressure arises in these elements).

Proceeding, past the option of printing PostScript files, you are now invited to proceed to the next result set. This will be the solution at t=5.0 . View the displacements, stresses and pore pressures as before. Repeat the process for the results at t=15.0 and t=35.0 .