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

Check out the videos below
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,
      <p style="margin-left: 93px;">plot_width=640, plot_height=640, pad_ratio=1.375, live_plotQ=False)</p>
      = fs_plot(<span style="color:red">P_load</span>, <span style="color:red">Lpts</span>, <span style="color:red">Load_Point_X</span>, <span style="color:red">Load_Point_Y</span>, 
      <p style="margin-left: 66px;"><span style="color:red">plot_width</span>, <span style="color:red">plot_height</span>, <span style="color:red">pad_ratio</span>, <span style="color:red">live_plotQ</span>)</p>
      <span style="color:red">P_load</span>, <span style="color:red">Load_Point_X</span>, and <span style="color:red">Load_Point_Y</span> are single cell named ranges
      <p style="margin-left: 25px;">(the load magnitude)</p>
      <span style="color:red">Load_Point_X</span>, and <span style="color:red">Load_Point_Y</span> are single cell named ranges
      <p style="margin-left: 25px;">(the load X, Y coordinates)</p>
      <span style="color:red">Lpts</span> is a 2D named range
      <p style="margin-left: 25px;">(these are the X, Y coordinates of the polygon)</p>
      <span style="color:red">plot_width</span> is an optional integer input
      <span style="color:red">plot_height</span> is an optional integer input
      <span style="color:red">pad_ratio</span> is an optional float input
      <p style="margin-left: 25px;">(allows padding around the plotted polygon)</p>
      <span style="color:red">live_plotQ</span> 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”
El 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:

<span style="color:blue">class</span> Section():

    <span style="color:blue">def</span> <span style="color:magenta">trial_fxy</span>(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

        <span style="color:blue">if</span> pts[0] != pts[-1]:
            pts = pts + pts[:1]

        x = [ c.x <span style="color:blue">for</span> c <span style="color:blue">in</span> pts ]
        y = [ c.y <span style="color:blue">for</span> c <span style="color:blue">in</span> pts ]

        sx = sy = sx1 = sy1 = sk = 0

        <span style="color:blue">for</span> i <span style="color:blue">in</span> 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
        <span style="color:blue">return</span> 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
circle 0pt70.html

Leave a Reply

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