Foundation Bearing Pressures using Jupyter and PyXLL (with videos) Mar 2023

Foundation Bearing Pressures using Jupyter and PyXLL (with videos) Mar 2023

Check out the videos below
May 2 2023 Update

Since early March, off and on (mostly on ๐Ÿ˜Š), I have been working on improvements to this application. The iterative method is exposed to a surprising number of corner cases; I offer a short summary below.

Limitations and corner cases

  • Polygon: defines the extent of the base of the foundation
  • Load point (LP): defines the point of application of the vertical load
  • Center of gravity (CG): defines the CG of a polygon
  • Neutral axis (NA): defines the line of zero stress for the reduced polygon
    Simple example of a rectangular footing showing reduced bearing area
  • Reduced polygon: defines the effective portion of the base of the foundation
  • Quadrant: mathematics term

The user enters:

  • the coordinates of a polygon, in clockwise or counterclockwise order
  • the load point

The method of analysis is iterative.

The iterative method requires that the polygon be located in the first quadrant of the plane. The program takes care of that internally; the user can locate the polygon on any XY coordinate system.

The CG of the input polygon is calculated.

Within the polygon, the iterative procedure will tend to converge if the LP is such that

    \[\begin{aligned} &{X}_{LP} < {X}_{CG} \text{$\;\;$and} \\ &{Y}_{LP} < {Y}_{CG} \end{aligned}\]

that is, if the LP is lower left of the CG.

Internally, therefore, the program must rearrange the coordinate system to place the load point in the as-needed position. Turns out there are 8 possible transformations for an XY coordinate system (right click image to see in separate tab).


After the iteration is complete, the program returns:

  • the load point
  • the input polygon
  • the input polygon CG
  • the input polygon area
  • the neutral axis of the reduced polygon
  • the reduced polygon
  • the reduced polygon CG
  • the reduced polygon area
  • the reduced area ratio

Hover over chart elements to see analysis results.
Click legend elements to toggle on an off.
Static Charts (right click chart to see in separate tab)

Foundation Bearing Pressures
I built a Jupyter notebook to calculate foundation bearing pressures.

My Jupyter notebook can be used stand-alone or in conjunction with PyXLL.

Here are a few of the features.

  • The footing is defined using a simple list of 2D coordinates, a closed polygon
  • The single concentrated loading point is defined using simple 2d coordinates, a single point
  • I built a Point class to allow easy manipulation of 2D coordinates
  • I built a Points class to allow easy manipulation of lists of 2D coordinates
  • I built a Section class with methods for determining and manipulating section properties such as determining the resultant polygon after a given line (neutral axis) makes multiple intersections in a given polygon (section perimeter).
  • I built a Solver function and a Plot function to simplify model analysis (one line for each function)
  • I used Plotly for my plot routine
    • The as-built foundation is plotted: coordinates, area and CG are shown
    • The effective foundation is plotted: coordinates, area and CG are shown
    • The neutral axis is shown
    • The load point is shown
    • The resultant stress are shown at each vertex
    • PyXLL displays a static image of the Plotly chart
    • The static image is displayed below the function call
    • An optional LIVE Plotly HTML chart is built (see function call below)
    • The optional LIVE Plotly HTML chart may be viewed using an Excel hyperlink
      • Allows verification of all geometry
      • Allows verification of loading
      • For extra clarity, all plot items can be toggled
    • The Excel function call is a one liner:
      = def fs_plot(P_load, l_pts, load_point_x, load_point_y,
                    plot_width=640, plot_height=640, 
      			  pad_ratio=1.375, live_plotQ=False)
      = fs_plot(P_load, Lpts, Load_Point_X, Load_Point_Y, 
                plot_width, plot_height, 
      		  pad_ratio, live_plotQ)
      • P_load Load_Point_X, and Load_Point_Y are single cell named ranges
        (the load magnitude)
      • Load_Point_X, and Load_Point_Y are single cell named ranges
        (the load X, Y coordinates)
      • Lpts is a 2D named range
        (these are the X, Y coordinates of the polygon)
      • plot_width is an optional integer input
      • plot_height is an optional integer input
      • pad_ratio is an optional float input
        (allows padding around the plotted polygon)
      • live_plotQ is an optional boolean input

        Click for example live Chart! — hover over load point, vertices, CGs.

        The plot elements can be toggled off and on (click on the legend). The plot can be panned and zoomed.


        1) My calculations are based on an article by Eli Czerniak.

        Computer Foundation Design
        “How to Calculate Footing Soil Bearing by Computer”
        Eli Czerniak, Fluor Corporation
        Gulf Publishing Company, Houston, Texas

        This was a great article 60 years ago, and remains so well into the 21st century.

        2) The excellent article by Dr. Drang presenting the use of line integrals to calculate section properties. In my Point class, I used the basic functions as given by Dr. Drang with very slight modifications. I added a couple functions for my specific calculations.

        For example, here is my trial_fxy method; this method uses line integrals to determine the ‘CG’ of the uniformly varying stress on the effective area of the foundation. With quite a bit of help from Sympy ๐Ÿ˜Š, I built this function to check my results. The neutral axis and the associated dimensions a and b are shown in the reference diagram below.

        Using area integrals:

            \[f(x,y) = 1-\frac{x}{a}-\frac{y}{b}\]

            \[\displaystyle s = \mathop{\int\int} f(x,y) \,dx \,dy\]

            \[\displaystyle sx = \mathop{\int\int} x f(x,y) \,dx \,dy\]

            \[\displaystyle sy = \mathop{\int\int} y f(x,y) \,dx \,dy\]

            \[\displaystyle StressCG_x = \frac{sx}{s}\]

            \[\displaystyle StressCG_y = \frac{sy}{s}\]

        Therefore, using line integrals:

            \[\displaystyle s = \mathop{\oint} [-y^2/(2 b) + y (a - x)/a] \,dx\]

            \[\displaystyle sx = \mathop{\oint} [-x y^2/(2 b) + y(a x - x^2)/a] \,dx\]

            \[\displaystyle sy = \mathop{\oint} [-y^3/(3 b) + y^2 (a - x)/(2 a)] \,dx\]

            \[\displaystyle StressCG_x = \frac{sx}{s}\]

            \[\displaystyle StressCG_y = \frac{sy}{s}\]

        The Python code:

        class Section():
            def trial_fxy(pts, a, b):
                '''Input: polygon as list of points, pts
                          Neutral axis x-intercept, a
                          Neutral axis y-intercept, b
                          points are type Point
                # polygon must be a 'closed' set of coordinates   
                # therefore, the last point must equal the first point
                if pts[0] != pts[-1]:
                    pts = pts + pts[:1]
                x = [ c.x for c in pts ]
                y = [ c.y for c in pts ]
                sx = sy = sx1 = sy1 = sk = 0
                for i in range(len(pts) - 1):
                    x0 = x[i]
                    x1 = x[i+1]
                    y0 = y[i]
                    y1 = y[i+1]
                    s = ((-2*a*b*x0*y0 + 2*a*b*x1*y0 + 
                          a*x0*y0**2 - a*x1*y0**2 + 2*b*x0**2*y0 -
                          2*b*x0*x1*y0)/(2*a*b) + 
                          (a*x0*y0**2 - 
                          2*a*x0*y0*y1 + a*x0*y1**2 -
                          a*x1*y0**2 + 2*a*x1*y0*y1 - a*x1*y1**2 + 
                          2*b*x0**2*y0 - 2*b*x0**2*y1 -
                          4*b*x0*x1*y0 + 4*b*x0*x1*y1 + 2*b*x1**2*y0 -
                          2*b*x1**2*y1)/(6*a*b) +
                          (a*b*x0*y0 - a*b*x0*y1 - a*b*x1*y0 + 
                          a*b*x1*y1 - a*x0*y0**2 + a*x0*y0*y1 +
                          a*x1*y0**2 - a*x1*y0*y1 - 2*b*x0**2*y0 + 
                          b*x0**2*y1 + 3*b*x0*x1*y0 -
                          b*x0*x1*y1 - b*x1**2*y0)/(2*a*b))
                    ssxx = ((-2*a*b*x0**2*y0 + 2*a*b*x0*x1*y0 + 
                             a*x0**2*y0**2 - a*x0*x1*y0**2 +
                             2*b*x0**3*y0 - 2*b*x0**2*x1*y0)/(2*a*b) +
                             (4*a*b*x0**2*y0 - 2*a*b*x0**2*y1 -
                             6*a*b*x0*x1*y0 + 2*a*b*x0*x1*y1 + 
                             2*a*b*x1**2*y0 - 3*a*x0**2*y0**2 +
                             2*a*x0**2*y0*y1 + 4*a*x0*x1*y0**2 - 
                             2*a*x0*x1*y0*y1 - a*x1**2*y0**2 -
                             6*b*x0**3*y0 + 2*b*x0**3*y1 + 
                             10*b*x0**2*x1*y0 - 2*b*x0**2*x1*y1 -
                             4*b*x0*x1**2*y0)/(4*a*b) + 
                             (-a*x0**2*y0**2 + 2*a*x0**2*y0*y1 - 
                             a*x0**2*y1**2 + 2*a*x0*x1*y0**2 - 
                             4*a*x0*x1*y0*y1 + 2*a*x0*x1*y1**2 - 
                             a*x1**2*y0**2 + 2*a*x1**2*y0*y1 - 
                             a*x1**2*y1**2 - 2*b*x0**3*y0 + 
                             2*b*x0**3*y1 + 6*b*x0**2*x1*y0 -
                             6*b*x0**2*x1*y1 - 6*b*x0*x1**2*y0 + 
                             6*b*x0*x1**2*y1 + 2*b*x1**3*y0 -
                             2*b*x1**3*y1)/(8*a*b) + 
                             (-2*a*b*x0**2*y0 +
                             2*a*b*x0**2*y1 + 4*a*b*x0*x1*y0 -
                             4*a*b*x0*x1*y1 - 2*a*b*x1**2*y0 + 
                             2*a*b*x1**2*y1 + 3*a*x0**2*y0**2 - 
                             4*a*x0**2*y0*y1 + a*x0**2*y1**2 - 
                             5*a*x0*x1*y0**2 + 6*a*x0*x1*y0*y1 - 
                             a*x0*x1*y1**2 + 2*a*x1**2*y0**2 -
                             2*a*x1**2*y0*y1 + 6*b*x0**3*y0 - 
                             4*b*x0**3*y1 - 14*b*x0**2*x1*y0 + 
                             8*b*x0**2*x1*y1 + 10*b*x0*x1**2*y0 - 
                             4*b*x0*x1**2*y1 - 2*b*x1**3*y0)/(6*a*b))
                    ssyy = ((-3*a*b*x0*y0**2 + 3*a*b*x1*y0**2 + 
                             2*a*x0*y0**3 - 2*a*x1*y0**3 +
                             3*b*x0**2*y0**2 - 
                             3*b*x0*x1*y0**2)/(6*a*b) + 
                             (2*a*b*x0*y0**2 - 2*a*b*x0*y0*y1 -
                             2*a*b*x1*y0**2 + 2*a*b*x1*y0*y1 - 
                             2*a*x0*y0**3 + 2*a*x0*y0**2*y1 + 
                             2*a*x1*y0**3 - 2*a*x1*y0**2*y1 - 
                             3*b*x0**2*y0**2 + 2*b*x0**2*y0*y1 + 
                             4*b*x0*x1*y0**2 - 2*b*x0*x1*y0*y1 -
                             b*x1**2*y0**2)/(4*a*b) + 
                             (-2*a*x0*y0**3 +
                             6*a*x0*y0**2*y1 - 6*a*x0*y0*y1**2 +
                             2*a*x0*y1**3 + 2*a*x1*y0**3 - 
                             6*a*x1*y0**2*y1 + 6*a*x1*y0*y1**2 - 
                             2*a*x1*y1**3 - 3*b*x0**2*y0**2 + 
                             6*b*x0**2*y0*y1 - 3*b*x0**2*y1**2 + 
                             6*b*x0*x1*y0**2 - 12*b*x0*x1*y0*y1 + 
                             6*b*x0*x1*y1**2 - 3*b*x1**2*y0**2 + 
                             6*b*x1**2*y0*y1 - 
                             3*b*x1**2*y1**2)/(24*a*b) + 
                             (-a*b*x0*y0**2 + 2*a*b*x0*y0*y1 - 
                             a*b*x0*y1**2 + a*b*x1*y0**2 - 
                             2*a*b*x1*y0*y1 + a*b*x1*y1**2 + 
                             2*a*x0*y0**3 - 4*a*x0*y0**2*y1 +
                             2*a*x0*y0*y1**2 - 2*a*x1*y0**3 + 
                             4*a*x1*y0**2*y1 - 2*a*x1*y0*y1**2 + 
                             3*b*x0**2*y0**2 - 4*b*x0**2*y0*y1 + 
                             b*x0**2*y1**2 - 5*b*x0*x1*y0**2 + 
                             6*b*x0*x1*y0*y1 - b*x0*x1*y1**2 +
                             2*b*x1**2*y0**2 - 2*b*x1**2*y0*y1)/(6*a*b))
                    sx += ssxx
                    sy += ssyy
                    sk += s
                return sx/sk, sy/sk

        What is PyXLL-Jupyter?

        Integration for Jupyter notebooks and Microsoft Excel

        PyXLL to Jupyter Interface (with video) Dec 2020

        What is PyXLL?

        Write Excel Add-Ins in Python

Example LIVE Rectangular Analysis Output
Example LIVE Circular Analysis Output
For some reason, this chart is ‘zoomed out’ in iOS (both Safari and Chrome); I haven’t tried MacOS.
In Windows, the chart appears as expected.
circle 0pt70.html

Leave a Reply

Your email address will not be published. Required fields are marked *