Retaining wall analysis – FEA and hand calculations

I have used the PlateMC spreadsheet described in the previous post to compare the results of a finite element analysis of a retaining wall with a similar analysis using the commercial package Strand7, and a simple hand calculation using Coulomb soil pressure theory.

The spreadsheet may be downloaded from:

The finite element model used is shown in the screen-shot below, showing the deformed shape under final load, with 100 times magnification of displacements:


The main features of the analysis were:

  • The retaining wall was a 10 metre high reinforced concrete structure, with an 8 metre wide base.
  • Different soil properties were assigned to the foundations, the fill elements immediately behind the wall, the remaining elements over the wall heel, and the remaining fill.
  • The fill was placed in 2 metre high layers, with a 20 kPa compaction load applied and then removed from each layer.
  • The final load case included a 30 kPa surcharge load, over the full width of the fill top surface.
  • The spreadsheet and Strand7 analyses used the same soil properties, except that the first spreadsheet run included a dilatency angle for the fill elements.

The hand calculation used a fixed active pressure coefficient, Ka, based on the standard Coulomb equation, using the properties of the fill over the heel:


For this structure the wall back face angle to the vertical (omega), and the top fill surface angle to the horizontal (beta) were both zero.  The fill friction angle (phi) and the soil/wall interface friction (delta) are listed in the screen shots below, together with the fill elastic modulus values used in the finite element analyses.

The results shown below are all for the completed fill plus 30 kPa surcharge load.  The first analysis used the same elastic modulus value for all the fill (50 MPa), and equal friction values of 35 degrees for the structural fill and the friction layer:


The Plate MC and Strand7 results were similar over most of the wall, but the Strand7 results increased more quickly towards the base, with the moment at the base being about 8% higher.  The hand calculation was about 20% higher than the Strand7 results.

Results of setting the soil dilatency angle to zero in the PlateMC analysis are shown below:


Bending moments at the base were slightly increased, but the results are still about 6% lower than the Strand7 results.

Reducing the elastic modulus of the general fill to 20 MPa slightly increased the PlateMC bending moment, and reduced the Strand7 moment, such that the results at the base were very close:PlateMC-RW1-3

Reducing the friction angle of the friction layer (and the hand calculation delta angle) to 17.5 degrees increased both finite element analyses by about 20%, with the PlateMC results now being a little higher.  The hand calculation moment was reduced a little (about 1%), reducing it to less than the FEA results.


Reducing the elastic modulus of the friction layer to 20 MPa had only a small effect, increasing the Strand7 moment slightly and reducing it for PlateMC:


Reducing the wall friction to zero in the hand calculation increased the bending moment by about 8%, to a value just over that found in the finite element analyses with a 17.5 degree friction layer:PlateMC-RW1-5a

Reducing the friction angle of the friction layer to 10 degrees increased the maximum bending moment in the FEA results by about 20%:PlateMC-RW1-6

Removing the compaction loads from the finite element analyses had only a small effect on the results, slightly reducing maximum bending moments for both programs:PlateMC-RW1-7

In summary:

  • Results from the spreadsheet, PlateMC_staged, and Strand7 were reasonably consistent, especially with reduced friction angles for the friction layer.
  • Maximum bending moments from the hand calculation were significantly higher than the finite element results when the full friction angle was used.
  • Reducing the friction layer friction had much less effect on the hand calculation results than the FEA results.
  • Reasonably consistent results were found (in this case) when the FEA results with friction zone angle of 17.5 degrees were compared with the hand calculation results with zero interface friction.
Posted in Excel, Finite Element Analysis, Fortran, Frame Analysis, Geotechnical Engineering, Link to dll, Link to Python, Newton, NumPy and SciPy, Strand7, UDFs, VBA | Tagged , , , , , , , , , , , , , | Leave a comment

Plane Strain FEA – Non linear staged analysis

Continuing the recent series of posts featuring Fortran based FEA routines, the spreadsheet from the previous post has been revised:

  • The spreadsheet allows for a staged analysis, with additional elements placed in layers for any number of construction stages.
  • Temporary loads may be applied to the top surface of each layer, and then removed.
  • The mesh details are entered as node coordinates, and plate node numbers, rather than the automated mesh generation methods used previously.
  • The PlateStress4 user defined function (UDF) has been added, allowing beam forces and moments to be extracted from the Gauss Point stresses of plate elements.

The new file, including full open-source code, may be downloaded from:

As before, the spreadsheet requires Python (including Numpy) and xlwings, all of which are included in the free Anaconda Python package.

Example input and results are shown below, featuring analysis of a 10 metre high cantilever retaining wall, with fill placed in 5 stages:


Any number of materials may be defined.  Properties are:

  • E: Young’s Modulus
  • v: Poisson’s Ratio
  • c: Cohesion
  • phi: Friction angle
  • psi: Dilatency angle
  • gamma: density (force units)

All plates are 8-noded plane strain elements, with nodes defined in the clockwise direction, starting at the bottom left corner.


The construction sequence is defined with the last plate number in each applied layer (column W).  Both plate and node numbers must form a continuous sequence in each layer.  Compaction loads (column Z) are applied to then removed from the top surface of the plates listed in columns X and Y.  The final layer may also have a surcharge load applied to the specified plates as the final load case.

Restrained nodes are indicated with 0 for restrained freedoms and 1 for unrestrained.

Any number of pairs of nodes may be pinned, forcing equal X and Y deflections.


Typical results for the example structure are shown above.  The next post will provide more details of results of different soil properties, and also compare results with those from a commercial FEA package (Strand7), and standard retaining wall analysis procedures.


Numerical summary results are shown above, including results after placing each layer, and after application and removal of the compaction load.  Total execution time was just over 5 seconds.


The spreadsheet is set up to calculate the axial force, bending moment, and shear force over the height of the wall for any specified load case.


Calculation of the wall actions is carried out on the ConcRes sheet.  The procedure to generate these results is:

  • List the element numbers where results are required, with 4 rows for each element.
  • Calculate the row and column numbers from the StressRes range for the required Gauss points and load case.
  • Use the Index function to return the required stress values.
  • Use the PlateStress4 function to calculate the axial force, bending moment, and shear force for each of the selected plates (see previous post)
Posted in Beam Bending, Concrete, Excel, Finite Element Analysis, Fortran, Geotechnical Engineering, Link to dll, Link to Python, Newton, NumPy and SciPy, UDFs, VBA | Tagged , , , , , , , , , , , | 1 Comment

PlateStress Update

The functions from the PlateStress spreadsheet (details here) will be incorporated in the FEA spreadsheets currently being presented.

In the process of using these functions I made a couple of improvements:

  • The required sequence of the Gauss Points numbers was clarified:
  • The plate angle input was modified, so that for a straight beam a single value would be accepted:

Download from

Posted in Beam Bending, Excel, Finite Element Analysis, Frame Analysis, Newton, Strand7, UDFs, VBA | Tagged , , , , , , | 1 Comment

Dave Swarbrick

Dave Swarbrick shared with Mark Twain the privilege of saying that reports of his death had been greatly exaggerated.  Sadly, that is no longer true; he died on 3rd June this year, a fact not widely reported in Australia (although he did get a decent obituary in The Sydney Morning Herald).

Here follows a small sample of his work, including two live BBC recordings with the guitarist Richard Thompson, one from 1969 with Fairport Convention, and finally  from the Goodbye Television Centre Concert in 2013:


Posted in Bach | Tagged , , | Leave a comment

FEA slope analysis; effect of mesh size

The last version of the FEA slope analysis spreadsheet (presented here) allows the slope shape and element density to be easily varied.  This post looks at the effect of varying the density from very coarse to fine, with constant soil properties.

The element density was varied from a minimum of 24, up to a maximum of 5600, as shown in the screenshots below:

24 Elements:

5600 Elements:PlaneStrain5-3

Calculated maximum deflections with increasing reduction of the soil strength properties (friction angle and cohesion) are shown below:
It can be seen that:

  • The coarsest mesh (24 element) shows almost none of the non-linear behaviour found in the other analyses.
  • The 80 element mesh is significantly better, and the 350 element mesh (shown below) gives almost the same results as the two finer meshes.
  • The 1400 and 5600 element meshes gave almost the same results
  • The finest mesh gave the best visual representation of the circular arc of the slip failure surface.

Note that the analysis does not model non-linear geometric effects, and the deformations at the base of the slope are not realistic; nonetheless geometric effects (which are magnified by a factor of 100 in the diagrams) have only a small effect on the overall behaviour.

350 Element mesh:PlaneStrain5-1

Posted in Excel, Finite Element Analysis, Fortran, Geotechnical Engineering, Link to dll, Link to Python, Newton, VBA | Tagged , , , , , , , , , , | Leave a comment

Linking Excel to Fortran with ctypes; update

In the previous post arrays were created in Python and passed to a Fortran dll as ctype pointers.  The data written to these arrays in the Fortran routines was then converted back to Python arrays with statements such as: nodedisp = npc.as_array(nodedisp).  In fact this conversion is not necessary.  The data passed to the Fortran code is a pointer to the location of the arrays created in Python.  The Fortran results are therefore written to these arrays, which may be read directly from Python using the original array names.  More importantly, I found the results of using the npc.as_array() method to be inconsistent; sometimes returning the expected results, but sometimes empty or corrupted arrays.  Reading the results of the Fortran code using the original Python arrays is therefore both simpler and more reliable.

I have updated the previous post to reflect the revised process, and also updated the download file with the revised code.  The new file, including full open-source code, may be downloaded from:

Posted in Excel, Finite Element Analysis, Fortran, Link to dll, Link to Python, Newton, NumPy and SciPy, VBA | Tagged , , , , , , , , , , | 3 Comments

Linking Excel to Fortran, via xlwings, Python and ctypes

Continuing the recent series on using Fortran code published in Programming the Finite Element Method (5th ed. John Wiley & Sons, I.M. Smith, D.V. Griffiths and L. Margetts (2014)), this post looks at linking Excel to Fortran code compiled as a dll (dynamic link library).

It is quite possible to link Excel VBA directly to Fortran based dlls, but because Python provides a convenient link to numerous scientific and maths routines in Numpy, Scipy, and other libraries, I have chosen to create the link through Python, using xlwings.  The spreadsheet, dll files, and Python and Fortran code can be downloaded from:

To run the code requires Python and xlwings to be installed.  The Anaconda Python package includes xlwings by default, and also allows The MinGW package to be installed, which includes Gfortran.  GFortran is not required to run the application, but it is needed if you want to make any changes to the Fortran code and recompile.  The dll file may be copied to the same directory as the spreadsheet and Python file (, or to any directory on the system search path.

As for compiling an exe file and linking to Excel, the process is not difficult, but it is not well documented, and any deviation from the exact procedure is likely to result in errors accompanied with a collection of exceedingly unhelpful error messages.

In the previous version of this spreadsheet the data for the compiled program was written to the hard disk as text files, which were opened and read by the Fortran program, which then wrote back the results to the hard disk, to be opened and read by Excel.  This is a slow and cumbersome procedure, but it does remove all problems of differing data formats.  In contrast, when the data is transferred directly from Excel to the Fortran based dll it is essential all data and objects are in a format that is recognisable to both programs.  The procedure I have used to ensure data compatibility was:

  • In Excel VBA collect the input data as a series of named worksheet ranges, converted to variant arrays.
  • Transfer to Python with an xlwings UDF (user defined function).  Arrays are transferred as Python tuples.
  • Convert the tuples to Numpy arrays of the appropriate data type, and if necessary convert values to the correct data type.  For this application all values were either double64 or int (equivalent to VBA double and long).
  • Convert the arrays and values to c_type format, which is compatible with both C and Fortran.
  • Call the required function from the dll, passing all arguments by reference.
  • Return the result arrays and values to VBA using the original Python variables.
  • Paste the results to named ranges in the spreadsheet.
  • Update the deformed mesh plot in VBA.

The VBA code is straightforward.  The Sub P64py reads data from the ranges “Plate_def”, “Tolerance”, “srf”, and “props”, adjusting the range sizes where necessary, then calls the function xl_P64, which converts empty cells in the Plate_Def range to zeros, then calls the Python function py64 in the module

The main features of the Python code are:

  • Import numpy, os and xlwings.
  • Import the required data types and methods from the ctypes library (byref, cdll, c_int, c_double), and import numpy.ctypeslib to handle Numpy arrays.
  • Load the dll file, using cdll.LoadLibrary (see Python code below).  The code is set up to look for the dll first in the workbook directory, or if that fails look anywhere on the search path.
  • The data from the “Plate_def” worksheet range is passed from VBA as a variant array with all numerical values passed as doubles.  In both Numpy and Fortran arrays must have a single data type, and the Fortran code requires some of the input data as doubles and some as integers.  For this reason the dataa array must be split into two separate arrays, containing doubles (meshdims) and integers (meshnums).
  • Arrays for the output data must be sized and created in Python, and passed as ctype pointers to Fortran, together with the size of any variable dimensions.
  • Having converted and created the required arrays in Numpy format, ctype pointers are created with statements such as:
    coordp = npc.as_ctypes(coorda)
    Note that it is important that the pointer has a different name to the Numpy array.
  • Values are converted to ctypes using c_int or c_double, e.g:
    nsrf = c_int(nsrfi)
  • The p64 function is called from p64dll, with all arguments passed “byref” (see code below)
  • The data passed to the Fortran code consists of pointers to the Python variables, hence the arrays and values returned from p64 may be accessed using the original Python variables names.

The full Python code is shown below:

import numpy  as np
from ctypes import byref, cdll, c_int, c_double, c_char_p
import numpy.ctypeslib as npc
import os
import time
import xlwings as xw

#Load p64s.dll
    wkbpath = os.path.dirname(os.path.abspath(__file__))
    p64dll = cdll.LoadLibrary(wkbpath + '/p64s.dll')
    p64dll = cdll.LoadLibrary('./p64s.dll')
def py64(dataa, tola, srfa, propa):
# Convert dataa to Numpy array, and split into double and integer arrays
    dataa = np.array(dataa)
    meshdims = dataa[0:2,:]
    meshnums = dataa[2:4,0:2].astype(int)

#Calculate size of output arrays    
    np_types = int(dataa[4,0])
    nx1 = meshnums[0,0]
    nx2 = meshnums[0,1]
    ny1 = meshnums[1,0]
    ny2 = meshnums[1,1]
#Create Numpy arrays and convert to ctypes pointers
#Convert integers to c_int
    coorda = np.zeros((nn,2))
    coordp = npc.as_ctypes(coorda)
    g_num = np.zeros((nels,8),
    g_nump = npc.as_ctypes(g_num)
    srfa = np.array(srfa)[:,0]
    nsrfi = srfa.shape[0]
    nsrf = c_int(nsrfi)
    nodedisa = np.zeros((nn,nsrfi*3))
    nodedisp = npc.as_ctypes(nodedisa)
    nnp = c_int(nn)
    meshdimp = npc.as_ctypes(meshdims)
    meshnump = npc.as_ctypes(meshnums)
    tolp = npc.as_ctypes(np.array(tola)[0])
    srfp = npc.as_ctypes(srfa)
    propa = np.array(propa)
    propp = npc.as_ctypes(propa)
    resa = np.zeros((5,nsrfi))
    resp = npc.as_ctypes(resa)
    np_types = c_int(np_types)
    nelsp = c_int(nels)
 #Call p64 with all ctype arguments passed by reference   
    iErr = p64dll.p64(byref(resp), byref(meshdimp), byref(meshnump), byref(tolp), byref(srfp), byref(nsrf), 
    byref(propp), byref(np_types), byref(coordp), byref(nnp), byref(nodedisp), byref(g_nump), byref(nelsp))
 #Return results using original Python variables
    return [resa, coorda, nodedisa, g_num, nn,  nels] 

The Fortran code requires an interface to read input data and return results. Note that dynamic arrays must have their size defined by integer arguments passed from the calling program. Also note that Fortran arrays are defined in the opposite direction to Python and VBA arrays (i.e. columns, rows):

module main
implicit none
integer, parameter:: dp=SELECTED_REAL_KIND(15)


...  Functions called by p64

END module main

integer function p64(resa, meshdims, meshnums, tola, srfa, nsrf, prop, &
   np_types, g_coord, nn, NodeDisp,g_num, nels)
! Program 6.4 Plane strain slope stability analysis of an elastic-plastic
!             (Mohr-Coulomb) material using 8-node rectangular
!             quadrilaterals. Viscoplastic strain method.
  use main
  implicit none 
 integer, intent(in):: meshnums(2,2), nsrf, np_types, nn, nels
 integer, intent(out):: g_num(8,nels)
 real(dp),INTENT(IN)::meshdims(3,2), tola(3), srfa(nsrf), prop(7,np_types)
 real(dp),INTENT(OUT)::resa(nsrf,5), g_coord(2,nn), NodeDisp(nsrf*3,nn)


The Fortran code is compiled to a dll as shown below:

  • gfortran -shared -fno-underscoring -o p64.dll P64dll.f03

As for the compilation of an exe file, if the program is to run on a computer without MinGW it must be compiled with the -static option:

  • gfortran -shared -static -fno-underscoring -o p64s.dll P64dll.f03


The resulting output is shown below:


The VBA plotting routine has been re-written to allow changes in the number of elements to be plotted.



See the next post for further examples and discussion of results with varying element size.

Posted in Excel, Finite Element Analysis, Fortran, Geotechnical Engineering, Link to dll, Link to Python, Newton, NumPy and SciPy, VBA | Tagged , , , , , , , , , , | 2 Comments