The xlwings Scipy spreadsheet has been updated with a new example of the xl_SolveF function, that uses the Scipy Optimize root function.  The new spreadsheet can be downloaded from:

xlScipy3.zip

The new example uses a python function ic_calc (included in the download file) to find the ultimate reaction forces and moment for a group of bolts with specified load eccentricity and angle, and yield displacement and force.  The calculation procedure is detailed at:
CE591eccentric_shear_F13.pdf

The shear force on each bolt is non-linear, depending on the bolt strength and shear displacement, and the displacement in each bolt depends on the position of the centre of rotation.  The Scipy root function adjusts the centre of rotation coordinates and the maximum bolt displacement so that the reaction forces and moments are equal to the applied forces and moment:

Results can be compared with a VBA spreadsheet doing the same calculation but using the Excel Solver, from Yakpol’s Spreadsheet Solutions for Structural Engineering.

Note that to generate the same result as my spreadsheet

• The XY coordinates must be adjusted to give the same perpendicular eccentricity from the centroid of the bolts to the line of action of the force.
• The force angle is the angle to the x-axis, whereas my function uses the angle to the Y axis.

The xl_SolveF function passes data to the function to solved as two arguments:

• A 1D vector containing the variables
• A single array containing all other data

An interface function is therefore used to extract the input parameters for the ic_calc function in the correct format, as shown below:

def ic_check(IC, vals):   
    xloc = vals[1][0]
    yloc  = vals[1][1]
    vals2 =vals[1][2]
    
    res1 = xlic_calc(IC,xloc, yloc, vals2)
    
    totX = res1[0]
    totY = res1[1]
    totM = res1[2]
            
    Pu =IC[2]
    Pux = res1[3]
    Puy = res1[4]
    Mu = res1[5]
    
    return [totX-Pux, totY-Puy, totM-Mu]

@xw.func
@xw.ret(transpose = True)
def xlic_calc(IC, xloc, yloc, vals):
    xloc = np.array(xloc)
    yloc = np.array(yloc)
    ecc =vals[0]
    theta = np.radians(vals[1])
    deltamax = vals[2]
    Rult = vals[3]
    num_bolts = len(xloc)
    
    ICx = IC[0]
    ICy = IC[1]
    Pu = IC[2]

    Pux = Pu * np.sin(theta)
    Puy = Pu * np.cos(theta)
    Mu = Pu * (ecc-ICx*np.cos(theta)-ICy*np.sin(theta))
    
    xIC = xloc - ICx
    yIC = yloc - ICy
    
    di = np.sqrt((xIC*xIC)+(yIC*yIC))
    dmax = max(di)
    
    deltai = di/dmax * deltamax
    ri = Rult * (1-np.e**(-10.0*deltai))**0.55
    fx = yIC * ri / di
    fy = xIC * ri / di
    moment = ri * di
    totX = sum(fx)
    totY = sum(fy)
    totM = sum(moment)
    
    return [totX, totY, totM, Pux, Puy, Mu]