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.
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.
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.
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.
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.
Again, a finer discretization would produce more meaningful stress contours.
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.
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.
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.
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 .