Corrected embedment factors for gazetas stiffness formulas

geofound/capacity.py

@@ -1126,26 +1126,49 @@ def calc_norm_moment_bs_via_butterfield_et_al_1994(n, v, n_max, b, qv_max, qm_ma
     m = 0.5 * (2 * c * s * v / t - s ** 2 * np.sqrt((4 * c ** 2 * v ** 2 / (s ** 2 * t ** 2) -
         4 * (-n ** 4 / n_max ** 2 + 2 * n ** 3 / n_max - n ** 2 + v ** 2 / t ** 2) / s ** 2)))
 
-    m = c * s * v / t - s * (np.sqrt((c ** 2 - 1) * v ** 2 / t ** 2 + n ** 2 * (n - n_max) ** 2 / n_max ** 2))
+    m = c * s * v / t - s * (np.sqrt((c ** 2 - 1) * v ** 2 / t ** 2 + n ** 2 * (n - n_max) ** 2 / n_max ** 2)) / t
+    # m = s * np.sqrt(c ** 2 * v ** 2 + n ** 4 * t ** 2 - 2 * n ** 3 * t ** 2 - v ** 2) / t + c * s * v / t
+    return m
 
-    return -m
+
+def calc_norm_shear_bs_via_butterfield_et_al_1994(n, r, qv_max, qm_max):
+    c = 0.22
+    t = qv_max
+    s = qm_max
+    v2 = (n - 1) * n * s * t / (np.sqrt(2 * c * r * s * t - r ** 2 * t ** 2 - s ** 2))
+    return v2
+
+def calc_norm_shear_bs(n, m):
+    c = 0.22
+    u = (c ** 2 - 1) * m ** 2 + (n - 1) ** 2 * n ** 2
 
 
 def run():
     import matplotlib.pyplot as plt
     n = 0.5
     n_max = 1
-    b = 3
-    qv_max = 0.13
-    qm_max = 0.48
-    v = np.linspace(0, qv_max, 10)
+    b = 1
+    qv_max = 0.52
+    qm_max = 0.35
+    v = np.linspace(0.1, 0.14, 20)
     # m = calc_norm_moment_bs_via_butterfield_et_al_1994(n, v, n_max, b, qv_max, qm_max)
-    # plt.plot(v / n_max, m / (b * n_max))
+    # for i in range(len(m)):
+    #     print(i, v[i], m[i])
+    # plt.plot(v / n_max, m / (n_max))
     v = 0
     n = np.linspace(0.1, 0.9, 10)
-    m = calc_norm_moment_bs_via_butterfield_et_al_1994(n, v, n_max, b, qv_max, qm_max)
-    plt.plot(n / n_max, m / (b * n_max))
+    n = 0.5
+    r = np.logspace(-1, 2, 5)
+    qv = calc_norm_shear_bs_via_butterfield_et_al_1994(n, r, qv_max, qm_max)
+    print(qv)
+    qm = r * qv
+    plt.plot(qv, qm)
+    # m = calc_norm_moment_bs_via_butterfield_et_al_1994(n, v, n_max, b, qv_max, qm_max)
+    # plt.plot(n / n_max, m / (b * n_max))
     plt.show()
 
+
+
+
 if __name__ == '__main__':
     run()

geofound/stiffness.py

@@ -82,7 +82,7 @@ def calc_rotational_via_pais_1988(sl, fd, ip_axis='width', axis=None, a0=0.0, **
     return k_f_0 * n_emb
 
 
-def calc_rotational_via_gazetas_1991(sl, fd, ip_axis='width', axis="length", a0=0.0, f_contact=1.0, **kwargs):
+def calc_rotational_via_gazetas_1991(sl, fd, ip_axis='width', axis=None, a0=0.0, f_contact=1.0, **kwargs):
     """
     Rotation stiffness of foundation from Gazetas (1991) and Mylonakis et al. (2006)
 
@@ -126,24 +126,29 @@ def calc_rotational_via_gazetas_1991(sl, fd, ip_axis='width', axis="length", a0=
         i_by = fd.i_ll
         i_bx = fd.i_ww
     if (ip_axis == 'width' and len_dominant) or (ip_axis == 'length' and not len_dominant):
-        x_axis = True
+        xx_axis = True  # weaker rotation
     else:
-        x_axis = False
+        xx_axis = False
 
     v = sl.poissons_ratio
     n_emb = 1.0
-    if x_axis:
+    if xx_axis:
         k_rx = 1 - 0.2 * a0
         k_f_0 = (sl.g_mod / (1 - v) * i_bx ** 0.75 * (l / b) ** 0.25 * (2.4 + 0.5 * (b / l))) * k_rx
         if fd.depth > 0.0:
             dw = min(fd.height, fd.depth) * f_contact
-            n_emb = 1 + 1.26 * (dw / b) ** 0.6 * (1.5 + (dw / fd.depth) ** 1.9 * (b / l) ** -0.6)
-    else:
-        k_ry = 1 - 0.3 * a0
+            n_emb = 1 + 1.26 * (dw / b) * (1. + (dw / b) * (dw / fd.depth) ** -0.2 * (b / l) ** 0.5)
+    else:  # yy_axis (rotation about y-axis)
+        if v < 0.45:
+            k_ry = 1 - 0.3 * a0
+        else:
+            k_ry = 1 - 0.25 * a0 * (l / b) ** 0.3
         k_f_0 = (sl.g_mod / (1 - v) * i_by ** 0.75 * (3 * (l / b) ** 0.15)) * k_ry
         if fd.depth > 0.0:
             dw = min(fd.height, fd.depth) * f_contact
-            n_emb = 1 + 0.92 * (dw / b) ** 0.6 * (1.0 + (dw / fd.depth) ** -0.2 * (b / l) ** 0.5)
+            # Note that the original Gazetas (1991) paper has (dw / L) ** 1.9 * (dw / D) ** -0.6
+            # whereas Mylonakis has this form:
+            n_emb = 1 + 0.92 * (dw / b) ** 0.6 * (1.5 + (dw / fd.depth) ** 1.9 * (b / l) ** -0.6)
     return k_f_0 * n_emb