geofound/capacity.py
def capacity_sp_meyerhof_and_hanna_1978(sp, fd, verbose=0):
    """
    Calculates the two-layered foundation capacity according Meyerhof and Hanna (1978)

    :param sp: Soil profile object
    :param fd: Foundation object
    :param wtl: water table level
    :param verbose: verbosity
    :return: ultimate bearing stress
    """
    assert isinstance(sp, sm.SoilProfile)

    sl_0 = sp.layer(1)
    sl_1 = sp.layer(2)
    h0 = sp.layer_depth(2)
    gwl = sp.gwl

    sl_0.nq_factor_0 = (
    (np.tan(np.pi / 4 + np.deg2rad(sl_0.phi / 2))) ** 2 * np.exp(np.pi * np.tan(np.deg2rad(sl_0.phi))))
    if sl_0.phi == 0:
        sl_0.nc_factor_0 = 5.14
    else:
        sl_0.nc_factor_0 = (sl_0.nq_factor_0 - 1) / np.tan(np.deg2rad(sl_0.phi))
    sl_0.ng_factor_0 = (sl_0.nq_factor_0 - 1) * np.tan(1.4 * np.deg2rad(sl_0.phi))

    sl_1.nq_factor_1 = (
        (np.tan(np.pi / 4 + np.deg2rad(sl_1.phi / 2))) ** 2 * np.exp(np.pi * np.tan(np.deg2rad(sl_1.phi))))
    if sl_1.phi == 0:
        sl_1.nc_factor_1 = 5.14
    else:
        sl_1.nc_factor_1 = (sl_1.nq_factor_1 - 1) / np.tan(np.deg2rad(sl_1.phi))
    sl_1.ng_factor_1 = (sl_1.nq_factor_1 - 1) * np.tan(1.4 * np.deg2rad(sl_1.phi))

    if verbose:
        log("Nc: ", sl_1.nc_factor_1)
        log("Nq: ", sl_1.nq_factor_1)
        log("Ng: ", sl_1.ng_factor_1)

    sl_0.kp_0 = (np.tan(np.pi / 4 + np.deg2rad(sl_0.phi / 2))) ** 2
    sl_1.kp_1 = (np.tan(np.pi / 4 + np.deg2rad(sl_1.phi / 2))) ** 2

    # shape factors
    if sl_0.phi >= 10:
        sl_0.s_c_0 = 1 + 0.2 * sl_0.kp_0 * (fd.width / fd.length)
        sl_0.s_q_0 = 1.0 + 0.1 * sl_0.kp_0 * (fd.width / fd.length)
    else:
        sl_0.s_c_0 = 1 + 0.2 * (fd.width / fd.length)
        sl_0.s_q_0 = 1.0
    sl_0.s_g_0 = sl_0.s_q_0

    if sl_1.phi >= 10:
        sl_1.s_c_1 = 1 + 0.2 * sl_1.kp_1 * (fd.width / fd.length)
        sl_1.s_q_1 = 1.0 + 0.1 * sl_1.kp_1 * (fd.width / fd.length)
    else:
        sl_1.s_c_1 = 1 + 0.2 * (fd.width / fd.length)
        sl_1.s_q_1 = 1.0
    sl_1.s_g_1 = sl_1.s_q_1

    # Note: this method explicitly accounts for the foundation depth, so there are no depth factors
    # TODO: inclination factors, see doi.org/10.1139/t78-060

    # Capacity
    a = 1  # assumed to  be one but can range between 1.1 and 1.27 for square footings according to Das (1999) Ch 4
    s = 1

    r = 1 + (fd.width / fd.length)

    # put the same things before that condition
    # effective weight not in the soil object

    if gwl == 0:  # case 1: GWL at surface
        q_at_interface = sl_0.unit_bouy_weight * h0
        unit_eff_weight_0_at_fd_depth = sl_0.unit_bouy_weight
        unit_eff_weight_0_at_interface = sl_0.unit_bouy_weight
        unit_eff_weight_1_below_foundation = sl_1.unit_bouy_weight

    elif 0 < gwl <= fd.depth:  # Case 2: GWL at between foundation depth and surface
        q_at_interface = (sl_0.unit_dry_weight * gwl) + (sl_0.unit_bouy_weight * (h0 - gwl))
        q_d = (sl_0.unit_dry_weight * gwl) + (sl_0.unit_bouy_weight * (fd.depth - gwl))
        unit_eff_weight_0_at_fd_depth = q_d / fd.depth
        unit_eff_weight_0_at_interface = sl_0.unit_bouy_weight
        unit_eff_weight_1_below_foundation = sl_1.unit_bouy_weight

    elif fd.depth < gwl <= fd.width + fd.depth:
        if gwl < h0:  # Case 3: GWL at between foundation depth and foundation depth plus width, and GWL < layer 1 depth

            average_unit_bouy_weight = sl_0.unit_bouy_weight + (
            ((gwl - fd.depth) / fd.width) * (sl_0.unit_dry_weight - sl_0.unit_bouy_weight))

            q_at_interface = (sl_0.unit_dry_weight * gwl) + (sl_0.unit_bouy_weight * (h0 - gwl))
            unit_eff_weight_0_at_fd_depth = sl_0.unit_dry_weight
            unit_eff_weight_0_at_interface = average_unit_bouy_weight
            unit_eff_weight_1_below_foundation = sl_1.unit_bouy_weight

        else:  # Case 4: GWL at between foundation depth and foundation depth plus width, and GWL > layer 1 depth
            average_unit_bouy_weight = sl_1.unit_bouy_weight + (
            ((gwl - h0) / fd.width) * (sl_1.unit_dry_weight - sl_1.unit_bouy_weight))

            q_at_interface = sl_0.unit_dry_weight * h0
            unit_eff_weight_0_at_fd_depth = sl_0.unit_dry_weight
            unit_eff_weight_0_at_interface = sl_0.unit_dry_weight
            unit_eff_weight_1_below_foundation = average_unit_bouy_weight

    elif gwl > fd.depth + fd.width:  # Case 5: GWL beyond foundation depth plus width
        q_at_interface = sl_0.unit_dry_weight * h0
        unit_eff_weight_0_at_fd_depth = sl_0.unit_dry_weight
        unit_eff_weight_0_at_interface = sl_0.unit_dry_weight
        unit_eff_weight_1_below_foundation = sl_1.unit_dry_weight
    else:
        raise ValueError("Could not interpret inputs")  # never reached

    # maximum value (qu <= qt)
    q_ult6 = q_at_interface - unit_eff_weight_0_at_fd_depth * fd.depth
    q_0 = (sl_0.cohesion * sl_0.nc_factor_0) + (0.5 * unit_eff_weight_0_at_interface * fd.width * sl_0.ng_factor_0)
    q_b2 = (q_at_interface * sl_1.nq_factor_1 * sl_1.s_q_1)
    q_1 = (sl_1.cohesion * sl_1.nc_factor_1) + (0.5 * unit_eff_weight_1_below_foundation * fd.width * sl_1.ng_factor_1)
    q_b3 = (unit_eff_weight_1_below_foundation * fd.width * sl_1.ng_factor_1 * sl_1.s_g_1 / 2)
    q_ult5 = r * (unit_eff_weight_0_at_interface * ((h0 - fd.depth) ** 2)) * (1 + (2 * fd.depth / (h0 - fd.depth))) * (
        np.tan(np.deg2rad(sl_0.phi)) / fd.width) * s
    q_t2 = (unit_eff_weight_0_at_fd_depth * fd.depth * sl_0.nq_factor_0 * sl_0.s_q_0)
    q_t3 = (unit_eff_weight_0_at_interface * fd.width * sl_0.ng_factor_0 * sl_0.s_g_0 / 2)

    # qb
    q_b1 = (sl_1.cohesion * sl_1.nc_factor_1 * sl_1.s_c_1)
    q_b = q_b1 + q_b2 + q_b3

    q1_q0 = q_1 / q_0

    # calculate the ca factor
    # if sl_0.cohesion == 0:
    #     c1_c0 = 0
    # else:
    #     c1_c0 = sl_1.cohesion / sl_0.cohesion
    x = np.array([0.000, 0.082, 0.206, 0.298, 0.404, 0.509, 0.598, 0.685, 0.772])
    y = np.array([0.627, 0.700, 0.794, 0.855, 0.912, 0.948, 0.968, 0.983, 0.997])

    # raise Warning("ca should be interpolated using q1/q2 not cohesion, see Figure 4 in MH1978")
    ca_c0 = np.interp(q1_q0, x, y)
    ca = ca_c0 * sl_0.cohesion

    # ks
    x_0 = np.array([0, 20.08, 22.42, 25.08, 27.58, 30.08, 32.58, 34.92, 37.83, 40.00, 42.67, 45.00, 47.00, 49.75])
    y_0 = np.array([0.93, 0.93, 0.93, 0.93, 1.01, 1.17, 1.32, 1.56, 1.87, 2.26, 2.72, 3.35, 3.81, 4.82])
    x_2 = np.array([0, 20.08, 22.50, 25.08, 27.58, 30.08, 32.50, 35.00, 37.67, 40.17, 42.67, 45.00, 47.50, 50.00])
    y_2 = np.array([1.55, 1.55, 1.71, 1.86, 2.10, 2.33, 2.72, 3.11, 3.81, 4.43, 5.28, 6.14, 7.46, 9.24])
    x_4 = np.array([0, 20.00, 22.51, 25.10, 27.69, 30.11, 32.45, 35.04, 37.88, 40.14, 42.65, 45.07, 47.33, 50.08])
    y_4 = np.array([2.49, 2.49, 2.64, 2.87, 3.34, 3.81, 4.43, 5.20, 6.29, 7.38, 9.01, 11.11, 14.29, 19.34])
    x_10 = np.array([0, 20.00, 22.50, 25.08, 28.00, 30.00, 32.50, 34.92, 37.50, 40.17, 42.42, 45.00, 47.17, 50.08])
    y_10 = np.array([3.27, 3.27, 3.74, 4.44, 5.37, 6.07, 7.16, 8.33, 10.04, 12.30, 15.95, 21.17, 27.47, 40.00])
    x_int = sl_0.phi

    if sl_0.phi < 1:
        fd.ks = 0
    else:

        if q1_q0 == 0:
            fd.ks = np.interp(x_int, x_0, y_0)

        elif q1_q0 == 0.2:
            fd.ks = np.interp(x_int, x_2, y_2)

        elif q1_q0 == 0.4:
            fd.ks = np.interp(x_int, x_4, y_4)

        elif q1_q0 == 1.0:
            fd.ks = np.interp(x_int, x_10, y_10)

        elif 0 < q1_q0 < 0.2:
            ks_1 = np.interp(x_int, x_0, y_0)
            ks_2 = np.interp(x_int, x_2, y_2)
            fd.ks = (((ks_2 - ks_1) * q1_q0) / 0.2) + ks_1

        elif 0.2 < q1_q0 < 0.4:
            ks_1 = np.interp(x_int, x_2, y_2)
            ks_2 = np.interp(x_int, x_4, y_4)
            fd.ks = (((ks_2 - ks_1) * (q1_q0 - 0.2)) / 0.2) + ks_1

        elif 0.4 < q1_q0 < 1.0:
            ks_1 = np.interp(x_int, x_4, y_4)
            ks_2 = np.interp(x_int, x_10, y_10)
            fd.ks = (((ks_2 - ks_1) * (q1_q0 - 0.4)) / 0.6) + ks_1
        else:
            raise DesignError(
                "Cannot compute 'ks', bearing ratio out-of-range (q1_q0 = %.3f) required: 0-1." % q1_q0)

    # qu
    q_ult4 = (r * (2 * ca * (h0 - fd.depth) / fd.width) * a)
    q_ult5_ks = q_ult5 * fd.ks
    q_ult = q_b + q_ult4 + q_ult5_ks - q_ult6

    q_t1 = (sl_0.cohesion * sl_0.nc_factor_0 * sl_0.s_c_0)
    q_t = q_t1 + q_t2 + q_t3

    if q_ult > q_t:
        if h0 > fd.width/2:
            fd.q_ult = q_t

        else:
            vert_eff_stress_interface = sp.vertical_effective_stress(h0)
            vert_eff_stress_lowest = sp.vertical_effective_stress(fd.width+fd.depth)
            average_eff_stress = (vert_eff_stress_interface + vert_eff_stress_lowest) / 2

            c_2_eff = sl_1.cohesion + average_eff_stress * np.tan(np.radians(sl_1.phi))

            if sl_0.cohesion > c_2_eff:
                fd.q_ult = q_t

            else:
                # vd = {}
                # vd[1] =[1, 1, 1, 1, 1]
                # vd[0.667] = [1, 1.033, 1.064, 1.088, 1.109]
                # vd[0.5] = [1, 1.056, 1.107, 1.152, 1.193]
                # vd[0.333] = [1, 1.088, 1.167, 1.241, 1.311]
                # vd[0.25] = [1, 1.107, 1.208, 1.302, 1.389]
                # vd[0.2] = [1, 1.121, 1.235, 1.342, 1.444]
                # vd[0.1] = [1, 1.154, 1.302, 1.446, 1.584]

                h_over_b = (h0 - fd.depth)/fd.width
                c1_over_c2 =sl_0.cohesion/c_2_eff

                c_1_over_c_2 = [0.1, 0.2, 0.25, 0.333, 0.5, 0.667, 1.]
                m_1 = [1.584, 1.444, 1.389, 1.311, 1.193, 1.109, 1.]
                m_125 = [1.446, 1.342, 1.302, 1.241, 1.152, 1.088, 1.]
                m_167 = [1.302, 1.235, 1.208, 1.167, 1.107, 1.064, 1.]
                m_25 = [1.154, 1.121, 1.107, 1.088, 1.056, 1.033, 1.]
                m_5 = [1, 1, 1, 1, 1, 1, 1]

                if h_over_b == 0.1:
                    m = np.interp(c1_over_c2, c_1_over_c_2, m_1)
                elif h_over_b == 0.125:
                    m = np.interp(c1_over_c2, c_1_over_c_2, m_125)
                elif h_over_b == 0.167:
                    m = np.interp(c1_over_c2, c_1_over_c_2, m_167)
                elif h_over_b == 0.250:
                    m = np.interp(c1_over_c2, c_1_over_c_2, m_25)
                elif h_over_b >= 0.5:
                    m = np.interp(c1_over_c2, c_1_over_c_2, m_5)
                elif 0.1 < h_over_b < 0.125:
                    m_a = np.interp(c1_over_c2, c_1_over_c_2, m_1)
                    m_b = np.interp(c1_over_c2, c_1_over_c_2, m_125)
                    m = np.interp(h_over_b, [0.1,0.125], [m_a,m_b])
                elif 0.125 < h_over_b < 0.167:
                    m_a = np.interp(c1_over_c2, c_1_over_c_2, m_125)
                    m_b = np.interp(c1_over_c2, c_1_over_c_2, m_167)
                    m = np.interp(h_over_b, [0.125, 0.167], [m_a, m_b])
                elif 0.167 < h_over_b < 0.25:
                    m_a = np.interp(c1_over_c2, c_1_over_c_2, m_167)
                    m_b = np.interp(c1_over_c2, c_1_over_c_2, m_25)
                    m = np.interp(h_over_b, [0.167, 0.250], [m_a, m_b])
                elif 0.25 < h_over_b < 0.5:
                    m_a = np.interp(c1_over_c2, c_1_over_c_2, m_25)
                    m_b = np.interp(c1_over_c2, c_1_over_c_2, m_5)
                    m = np.interp(h_over_b, [0.250, 0.500], [m_a, m_b])

                fd.q_ult = (sl_0.cohesion * m * sl_0.nc_factor_0) + (unit_eff_weight_0_at_fd_depth * fd.depth)

    else:
        fd.q_ult = q_ult

    return fd.q_ult
def capacity_sp_meyerhof_and_hanna_1978(sp, fd, verbose=0):
    """
    Calculates the two-layered foundation capacity according Meyerhof and Hanna (1978)

    :param sp: Soil profile object
    :param fd: Foundation object
    :param wtl: water table level
    :param verbose: verbosity
    :return: ultimate bearing stress
    """
    assert isinstance(sp, sm.SoilProfile)

    sl_0 = sp.layer(1)
    sl_1 = sp.layer(2)
    h0 = sp.layer_depth(2)
    gwl = sp.gwl

    sl_0.nq_factor_0 = (
    (np.tan(np.pi / 4 + np.deg2rad(sl_0.phi / 2))) ** 2 * np.exp(np.pi * np.tan(np.deg2rad(sl_0.phi))))
    if sl_0.phi == 0:
        sl_0.nc_factor_0 = 5.14
    else:
        sl_0.nc_factor_0 = (sl_0.nq_factor_0 - 1) / np.tan(np.deg2rad(sl_0.phi))
    sl_0.ng_factor_0 = (sl_0.nq_factor_0 - 1) * np.tan(1.4 * np.deg2rad(sl_0.phi))

    sl_1.nq_factor_1 = (
        (np.tan(np.pi / 4 + np.deg2rad(sl_1.phi / 2))) ** 2 * np.exp(np.pi * np.tan(np.deg2rad(sl_1.phi))))
    if sl_1.phi == 0:
        sl_1.nc_factor_1 = 5.14
    else:
        sl_1.nc_factor_1 = (sl_1.nq_factor_1 - 1) / np.tan(np.deg2rad(sl_1.phi))
    sl_1.ng_factor_1 = (sl_1.nq_factor_1 - 1) * np.tan(1.4 * np.deg2rad(sl_1.phi))

    if verbose:
        log("Nc: ", sl_1.nc_factor_1)
        log("Nq: ", sl_1.nq_factor_1)
        log("Ng: ", sl_1.ng_factor_1)

    sl_0.kp_0 = (np.tan(np.pi / 4 + np.deg2rad(sl_0.phi / 2))) ** 2
    sl_1.kp_1 = (np.tan(np.pi / 4 + np.deg2rad(sl_1.phi / 2))) ** 2

    # shape factors
    if sl_0.phi >= 10:
        sl_0.s_c_0 = 1 + 0.2 * sl_0.kp_0 * (fd.width / fd.length)
        sl_0.s_q_0 = 1.0 + 0.1 * sl_0.kp_0 * (fd.width / fd.length)
    else:
        sl_0.s_c_0 = 1 + 0.2 * (fd.width / fd.length)
        sl_0.s_q_0 = 1.0
    sl_0.s_g_0 = sl_0.s_q_0

    if sl_1.phi >= 10:
        sl_1.s_c_1 = 1 + 0.2 * sl_1.kp_1 * (fd.width / fd.length)
        sl_1.s_q_1 = 1.0 + 0.1 * sl_1.kp_1 * (fd.width / fd.length)
    else:
        sl_1.s_c_1 = 1 + 0.2 * (fd.width / fd.length)
        sl_1.s_q_1 = 1.0
    sl_1.s_g_1 = sl_1.s_q_1

    # Note: this method explicitly accounts for the foundation depth, so there are no depth factors
    # TODO: inclination factors, see doi.org/10.1139/t78-060

    # Capacity
    a = 1  # assumed to  be one but can range between 1.1 and 1.27 for square footings according to Das (1999) Ch 4
    s = 1

    r = 1 + (fd.width / fd.length)

    # put the same things before that condition
    # effective weight not in the soil object

    if gwl == 0:  # case 1: GWL at surface
        q_at_interface = sl_0.unit_bouy_weight * h0
        unit_eff_weight_0_at_fd_depth = sl_0.unit_bouy_weight
        unit_eff_weight_0_at_interface = sl_0.unit_bouy_weight
        unit_eff_weight_1_below_foundation = sl_1.unit_bouy_weight

    elif 0 < gwl <= fd.depth:  # Case 2: GWL at between foundation depth and surface
        q_at_interface = (sl_0.unit_dry_weight * gwl) + (sl_0.unit_bouy_weight * (h0 - gwl))
        q_d = (sl_0.unit_dry_weight * gwl) + (sl_0.unit_bouy_weight * (fd.depth - gwl))
        unit_eff_weight_0_at_fd_depth = q_d / fd.depth
        unit_eff_weight_0_at_interface = sl_0.unit_bouy_weight
        unit_eff_weight_1_below_foundation = sl_1.unit_bouy_weight

    elif fd.depth < gwl <= fd.width + fd.depth:
        if gwl < h0:  # Case 3: GWL at between foundation depth and foundation depth plus width, and GWL < layer 1 depth

            average_unit_bouy_weight = sl_0.unit_bouy_weight + (
            ((gwl - fd.depth) / fd.width) * (sl_0.unit_dry_weight - sl_0.unit_bouy_weight))

            q_at_interface = (sl_0.unit_dry_weight * gwl) + (sl_0.unit_bouy_weight * (h0 - gwl))
            unit_eff_weight_0_at_fd_depth = sl_0.unit_dry_weight
            unit_eff_weight_0_at_interface = average_unit_bouy_weight
            unit_eff_weight_1_below_foundation = sl_1.unit_bouy_weight

        else:  # Case 4: GWL at between foundation depth and foundation depth plus width, and GWL > layer 1 depth
            average_unit_bouy_weight = sl_1.unit_bouy_weight + (
            ((gwl - h0) / fd.width) * (sl_1.unit_dry_weight - sl_1.unit_bouy_weight))

            q_at_interface = sl_0.unit_dry_weight * h0
            unit_eff_weight_0_at_fd_depth = sl_0.unit_dry_weight
            unit_eff_weight_0_at_interface = sl_0.unit_dry_weight
            unit_eff_weight_1_below_foundation = average_unit_bouy_weight

    elif gwl > fd.depth + fd.width:  # Case 5: GWL beyond foundation depth plus width
        q_at_interface = sl_0.unit_dry_weight * h0
        unit_eff_weight_0_at_fd_depth = sl_0.unit_dry_weight
        unit_eff_weight_0_at_interface = sl_0.unit_dry_weight
        unit_eff_weight_1_below_foundation = sl_1.unit_dry_weight
    else:
        raise ValueError("Could not interpret inputs")  # never reached

    # maximum value (qu <= qt)
    q_ult6 = q_at_interface - unit_eff_weight_0_at_fd_depth * fd.depth
    q_0 = (sl_0.cohesion * sl_0.nc_factor_0) + (0.5 * unit_eff_weight_0_at_interface * fd.width * sl_0.ng_factor_0)
    q_b2 = (q_at_interface * sl_1.nq_factor_1 * sl_1.s_q_1)
    q_1 = (sl_1.cohesion * sl_1.nc_factor_1) + (0.5 * unit_eff_weight_1_below_foundation * fd.width * sl_1.ng_factor_1)
    q_b3 = (unit_eff_weight_1_below_foundation * fd.width * sl_1.ng_factor_1 * sl_1.s_g_1 / 2)
    q_ult5 = r * (unit_eff_weight_0_at_interface * ((h0 - fd.depth) ** 2)) * (1 + (2 * fd.depth / (h0 - fd.depth))) * (
        np.tan(np.deg2rad(sl_0.phi)) / fd.width) * s
    q_t2 = (unit_eff_weight_0_at_fd_depth * fd.depth * sl_0.nq_factor_0 * sl_0.s_q_0)
    q_t3 = (unit_eff_weight_0_at_interface * fd.width * sl_0.ng_factor_0 * sl_0.s_g_0 / 2)

    # qb
    q_b1 = (sl_1.cohesion * sl_1.nc_factor_1 * sl_1.s_c_1)
    q_b = q_b1 + q_b2 + q_b3

    q1_q0 = q_1 / q_0

    # calculate the ca factor
    # if sl_0.cohesion == 0:
    #     c1_c0 = 0
    # else:
    #     c1_c0 = sl_1.cohesion / sl_0.cohesion
    x = np.array([0.000, 0.082, 0.206, 0.298, 0.404, 0.509, 0.598, 0.685, 0.772])
    y = np.array([0.627, 0.700, 0.794, 0.855, 0.912, 0.948, 0.968, 0.983, 0.997])

    # raise Warning("ca should be interpolated using q1/q2 not cohesion, see Figure 4 in MH1978")
    ca_c0 = np.interp(q1_q0, x, y)
    ca = ca_c0 * sl_0.cohesion

    # ks
    x_0 = np.array([0, 20.08, 22.42, 25.08, 27.58, 30.08, 32.58, 34.92, 37.83, 40.00, 42.67, 45.00, 47.00, 49.75])
    y_0 = np.array([0.93, 0.93, 0.93, 0.93, 1.01, 1.17, 1.32, 1.56, 1.87, 2.26, 2.72, 3.35, 3.81, 4.82])
    x_2 = np.array([0, 20.08, 22.50, 25.08, 27.58, 30.08, 32.50, 35.00, 37.67, 40.17, 42.67, 45.00, 47.50, 50.00])
    y_2 = np.array([1.55, 1.55, 1.71, 1.86, 2.10, 2.33, 2.72, 3.11, 3.81, 4.43, 5.28, 6.14, 7.46, 9.24])
    x_4 = np.array([0, 20.00, 22.51, 25.10, 27.69, 30.11, 32.45, 35.04, 37.88, 40.14, 42.65, 45.07, 47.33, 50.08])
    y_4 = np.array([2.49, 2.49, 2.64, 2.87, 3.34, 3.81, 4.43, 5.20, 6.29, 7.38, 9.01, 11.11, 14.29, 19.34])
    x_10 = np.array([0, 20.00, 22.50, 25.08, 28.00, 30.00, 32.50, 34.92, 37.50, 40.17, 42.42, 45.00, 47.17, 50.08])
    y_10 = np.array([3.27, 3.27, 3.74, 4.44, 5.37, 6.07, 7.16, 8.33, 10.04, 12.30, 15.95, 21.17, 27.47, 40.00])
    x_int = sl_0.phi

    if sl_0.phi < 1:
        fd.ks = 0
    else:

        if q1_q0 == 0:
            fd.ks = np.interp(x_int, x_0, y_0)

        elif q1_q0 == 0.2:
            fd.ks = np.interp(x_int, x_2, y_2)

        elif q1_q0 == 0.4:
            fd.ks = np.interp(x_int, x_4, y_4)

        elif q1_q0 == 1.0:
            fd.ks = np.interp(x_int, x_10, y_10)

        elif 0 < q1_q0 < 0.2:
            ks_1 = np.interp(x_int, x_0, y_0)
            ks_2 = np.interp(x_int, x_2, y_2)
            fd.ks = (((ks_2 - ks_1) * q1_q0) / 0.2) + ks_1

        elif 0.2 < q1_q0 < 0.4:
            ks_1 = np.interp(x_int, x_2, y_2)
            ks_2 = np.interp(x_int, x_4, y_4)
            fd.ks = (((ks_2 - ks_1) * (q1_q0 - 0.2)) / 0.2) + ks_1

        elif 0.4 < q1_q0 < 1.0:
            ks_1 = np.interp(x_int, x_4, y_4)
            ks_2 = np.interp(x_int, x_10, y_10)
            fd.ks = (((ks_2 - ks_1) * (q1_q0 - 0.4)) / 0.6) + ks_1
        else:
            raise DesignError(
                "Cannot compute 'ks', bearing ratio out-of-range (q1_q0 = %.3f) required: 0-1." % q1_q0)

    # qu
    q_ult4 = (r * (2 * ca * (h0 - fd.depth) / fd.width) * a)
    q_ult5_ks = q_ult5 * fd.ks
    q_ult = q_b + q_ult4 + q_ult5_ks - q_ult6

    q_t1 = (sl_0.cohesion * sl_0.nc_factor_0 * sl_0.s_c_0)
    q_t = q_t1 + q_t2 + q_t3

    if q_ult > q_t:
        if h0 > fd.width/2:
            fd.q_ult = q_t

        else:
            vert_eff_stress_interface = sp.vertical_effective_stress(h0)
            vert_eff_stress_lowest = sp.vertical_effective_stress(fd.width+fd.depth)
            average_eff_stress = (vert_eff_stress_interface + vert_eff_stress_lowest) / 2

            c_2_eff = sl_1.cohesion + average_eff_stress * np.tan(np.radians(sl_1.phi))

            if sl_0.cohesion > c_2_eff:
                fd.q_ult = q_t

            else:
                # vd = {}
                # vd[1] =[1, 1, 1, 1, 1]
                # vd[0.667] = [1, 1.033, 1.064, 1.088, 1.109]
                # vd[0.5] = [1, 1.056, 1.107, 1.152, 1.193]
                # vd[0.333] = [1, 1.088, 1.167, 1.241, 1.311]
                # vd[0.25] = [1, 1.107, 1.208, 1.302, 1.389]
                # vd[0.2] = [1, 1.121, 1.235, 1.342, 1.444]
                # vd[0.1] = [1, 1.154, 1.302, 1.446, 1.584]

                h_over_b = (h0 - fd.depth)/fd.width
                c1_over_c2 =sl_0.cohesion/c_2_eff

                c_1_over_c_2 = [0.1, 0.2, 0.25, 0.333, 0.5, 0.667, 1.]
                m_1 = [1.584, 1.444, 1.389, 1.311, 1.193, 1.109, 1.]
                m_125 = [1.446, 1.342, 1.302, 1.241, 1.152, 1.088, 1.]
                m_167 = [1.302, 1.235, 1.208, 1.167, 1.107, 1.064, 1.]
                m_25 = [1.154, 1.121, 1.107, 1.088, 1.056, 1.033, 1.]
                m_5 = [1, 1, 1, 1, 1, 1, 1]

                if h_over_b == 0.1:
                    m = np.interp(c1_over_c2, c_1_over_c_2, m_1)
                elif h_over_b == 0.125:
                    m = np.interp(c1_over_c2, c_1_over_c_2, m_125)
                elif h_over_b == 0.167:
                    m = np.interp(c1_over_c2, c_1_over_c_2, m_167)
                elif h_over_b == 0.250:
                    m = np.interp(c1_over_c2, c_1_over_c_2, m_25)
                elif h_over_b >= 0.5:
                    m = np.interp(c1_over_c2, c_1_over_c_2, m_5)
                elif 0.1 < h_over_b < 0.125:
                    m_a = np.interp(c1_over_c2, c_1_over_c_2, m_1)
                    m_b = np.interp(c1_over_c2, c_1_over_c_2, m_125)
                    m = np.interp(h_over_b, [0.1,0.125], [m_a,m_b])
                elif 0.125 < h_over_b < 0.167:
                    m_a = np.interp(c1_over_c2, c_1_over_c_2, m_125)
                    m_b = np.interp(c1_over_c2, c_1_over_c_2, m_167)
                    m = np.interp(h_over_b, [0.125, 0.167], [m_a, m_b])
                elif 0.167 < h_over_b < 0.25:
                    m_a = np.interp(c1_over_c2, c_1_over_c_2, m_167)
                    m_b = np.interp(c1_over_c2, c_1_over_c_2, m_25)
                    m = np.interp(h_over_b, [0.167, 0.250], [m_a, m_b])
                elif 0.25 < h_over_b < 0.5:
                    m_a = np.interp(c1_over_c2, c_1_over_c_2, m_25)
                    m_b = np.interp(c1_over_c2, c_1_over_c_2, m_5)
                    m = np.interp(h_over_b, [0.250, 0.500], [m_a, m_b])

                fd.q_ult = (sl_0.cohesion * m * sl_0.nc_factor_0) + (unit_eff_weight_0_at_fd_depth * fd.depth)

    else:
        fd.q_ult = q_ult

    return fd.q_ult
def capacity_sp_meyerhof_and_hanna_1978 ( sp , fd , verbose = 0 ) : assert isinstance ( sp , sm . SoilProfile ) sl_0 = sp . layer ( 1 ) sl_1 = sp . layer ( 2 ) h0 = sp . layer_depth ( 2 ) gwl = sp . gwl sl_0 . nq_factor_0 = ( ( np . tan ( np . pi / 4 + np . deg2rad ( sl_0 . phi / 2 ) ) ) ** 2 * np . exp ( np . pi * np . tan ( np . deg2rad ( sl_0 . phi ) ) ) ) if sl_0 . phi == 0 : sl_0 . nc_factor_0 = 5.14 else : sl_0 . nc_factor_0 = ( sl_0 . nq_factor_0 - 1 ) / np . tan ( np . deg2rad ( sl_0 . phi ) ) sl_0 . ng_factor_0 = ( sl_0 . nq_factor_0 - 1 ) * np . tan ( 1.4 * np . deg2rad ( sl_0 . phi ) ) sl_1 . nq_factor_1 = ( ( np . tan ( np . pi / 4 + np . deg2rad ( sl_1 . phi / 2 ) ) ) ** 2 * np . exp ( np . pi * np . tan ( np . deg2rad ( sl_1 . phi ) ) ) ) if sl_1 . phi == 0 : sl_1 . nc_factor_1 = 5.14 else : sl_1 . nc_factor_1 = ( sl_1 . nq_factor_1 - 1 ) / np . tan ( np . deg2rad ( sl_1 . phi ) ) sl_1 . ng_factor_1 = ( sl_1 . nq_factor_1 - 1 ) * np . tan ( 1.4 * np . deg2rad ( sl_1 . phi ) ) if verbose : log ( "Nc: " , sl_1 . nc_factor_1 ) log ( "Nq: " , sl_1 . nq_factor_1 ) log ( "Ng: " , sl_1 . ng_factor_1 ) sl_0 . kp_0 = ( np . tan ( np . pi / 4 + np . deg2rad ( sl_0 . phi / 2 ) ) ) ** 2 sl_1 . kp_1 = ( np . tan ( np . pi / 4 + np . deg2rad ( sl_1 . phi / 2 ) ) ) ** 2 # shape factors if sl_0 . phi >= 10 : sl_0 . s_c_0 = 1 + 0.2 * sl_0 . kp_0 * ( fd . width / fd . length ) sl_0 . s_q_0 = 1.0 + 0.1 * sl_0 . kp_0 * ( fd . width / fd . length ) else : sl_0 . s_c_0 = 1 + 0.2 * ( fd . width / fd . length ) sl_0 . s_q_0 = 1.0 sl_0 . s_g_0 = sl_0 . s_q_0 if sl_1 . phi >= 10 : sl_1 . s_c_1 = 1 + 0.2 * sl_1 . kp_1 * ( fd . width / fd . length ) sl_1 . s_q_1 = 1.0 + 0.1 * sl_1 . kp_1 * ( fd . width / fd . length ) else : sl_1 . s_c_1 = 1 + 0.2 * ( fd . width / fd . length ) sl_1 . s_q_1 = 1.0 sl_1 . s_g_1 = sl_1 . s_q_1 # Note: this method explicitly accounts for the foundation depth, so there are no depth factors # TODO: inclination factors, see doi.org/10.1139/t78-060 # Capacity a = 1 # assumed to  be one but can range between 1.1 and 1.27 for square footings according to Das (1999) Ch 4 s = 1 r = 1 + ( fd . width / fd . length ) # put the same things before that condition # effective weight not in the soil object if gwl == 0 : # case 1: GWL at surface q_at_interface = sl_0 . unit_bouy_weight * h0 unit_eff_weight_0_at_fd_depth = sl_0 . unit_bouy_weight unit_eff_weight_0_at_interface = sl_0 . unit_bouy_weight unit_eff_weight_1_below_foundation = sl_1 . unit_bouy_weight elif 0 < gwl <= fd . depth : # Case 2: GWL at between foundation depth and surface q_at_interface = ( sl_0 . unit_dry_weight * gwl ) + ( sl_0 . unit_bouy_weight * ( h0 - gwl ) ) q_d = ( sl_0 . unit_dry_weight * gwl ) + ( sl_0 . unit_bouy_weight * ( fd . depth - gwl ) ) unit_eff_weight_0_at_fd_depth = q_d / fd . depth unit_eff_weight_0_at_interface = sl_0 . unit_bouy_weight unit_eff_weight_1_below_foundation = sl_1 . unit_bouy_weight elif fd . depth < gwl <= fd . width + fd . depth : if gwl < h0 : # Case 3: GWL at between foundation depth and foundation depth plus width, and GWL < layer 1 depth average_unit_bouy_weight = sl_0 . unit_bouy_weight + ( ( ( gwl - fd . depth ) / fd . width ) * ( sl_0 . unit_dry_weight - sl_0 . unit_bouy_weight ) ) q_at_interface = ( sl_0 . unit_dry_weight * gwl ) + ( sl_0 . unit_bouy_weight * ( h0 - gwl ) ) unit_eff_weight_0_at_fd_depth = sl_0 . unit_dry_weight unit_eff_weight_0_at_interface = average_unit_bouy_weight unit_eff_weight_1_below_foundation = sl_1 . unit_bouy_weight else : # Case 4: GWL at between foundation depth and foundation depth plus width, and GWL > layer 1 depth average_unit_bouy_weight = sl_1 . unit_bouy_weight + ( ( ( gwl - h0 ) / fd . width ) * ( sl_1 . unit_dry_weight - sl_1 . unit_bouy_weight ) ) q_at_interface = sl_0 . unit_dry_weight * h0 unit_eff_weight_0_at_fd_depth = sl_0 . unit_dry_weight unit_eff_weight_0_at_interface = sl_0 . unit_dry_weight unit_eff_weight_1_below_foundation = average_unit_bouy_weight elif gwl > fd . depth + fd . width : # Case 5: GWL beyond foundation depth plus width q_at_interface = sl_0 . unit_dry_weight * h0 unit_eff_weight_0_at_fd_depth = sl_0 . unit_dry_weight unit_eff_weight_0_at_interface = sl_0 . unit_dry_weight unit_eff_weight_1_below_foundation = sl_1 . unit_dry_weight else : raise ValueError ( "Could not interpret inputs" ) # never reached # maximum value (qu <= qt) q_ult6 = q_at_interface - unit_eff_weight_0_at_fd_depth * fd . depth q_0 = ( sl_0 . cohesion * sl_0 . nc_factor_0 ) + ( 0.5 * unit_eff_weight_0_at_interface * fd . width * sl_0 . ng_factor_0 ) q_b2 = ( q_at_interface * sl_1 . nq_factor_1 * sl_1 . s_q_1 ) q_1 = ( sl_1 . cohesion * sl_1 . nc_factor_1 ) + ( 0.5 * unit_eff_weight_1_below_foundation * fd . width * sl_1 . ng_factor_1 ) q_b3 = ( unit_eff_weight_1_below_foundation * fd . width * sl_1 . ng_factor_1 * sl_1 . s_g_1 / 2 ) q_ult5 = r * ( unit_eff_weight_0_at_interface * ( ( h0 - fd . depth ) ** 2 ) ) * ( 1 + ( 2 * fd . depth / ( h0 - fd . depth ) ) ) * ( np . tan ( np . deg2rad ( sl_0 . phi ) ) / fd . width ) * s q_t2 = ( unit_eff_weight_0_at_fd_depth * fd . depth * sl_0 . nq_factor_0 * sl_0 . s_q_0 ) q_t3 = ( unit_eff_weight_0_at_interface * fd . width * sl_0 . ng_factor_0 * sl_0 . s_g_0 / 2 ) # qb q_b1 = ( sl_1 . cohesion * sl_1 . nc_factor_1 * sl_1 . s_c_1 ) q_b = q_b1 + q_b2 + q_b3 q1_q0 = q_1 / q_0 # calculate the ca factor # if sl_0.cohesion == 0: #     c1_c0 = 0 # else: #     c1_c0 = sl_1.cohesion / sl_0.cohesion x = np . array ( [ 0.000 , 0.082 , 0.206 , 0.298 , 0.404 , 0.509 , 0.598 , 0.685 , 0.772 ] ) y = np . array ( [ 0.627 , 0.700 , 0.794 , 0.855 , 0.912 , 0.948 , 0.968 , 0.983 , 0.997 ] ) # raise Warning("ca should be interpolated using q1/q2 not cohesion, see Figure 4 in MH1978") ca_c0 = np . interp ( q1_q0 , x , y ) ca = ca_c0 * sl_0 . cohesion # ks x_0 = np . array ( [ 0 , 20.08 , 22.42 , 25.08 , 27.58 , 30.08 , 32.58 , 34.92 , 37.83 , 40.00 , 42.67 , 45.00 , 47.00 , 49.75 ] ) y_0 = np . array ( [ 0.93 , 0.93 , 0.93 , 0.93 , 1.01 , 1.17 , 1.32 , 1.56 , 1.87 , 2.26 , 2.72 , 3.35 , 3.81 , 4.82 ] ) x_2 = np . array ( [ 0 , 20.08 , 22.50 , 25.08 , 27.58 , 30.08 , 32.50 , 35.00 , 37.67 , 40.17 , 42.67 , 45.00 , 47.50 , 50.00 ] ) y_2 = np . array ( [ 1.55 , 1.55 , 1.71 , 1.86 , 2.10 , 2.33 , 2.72 , 3.11 , 3.81 , 4.43 , 5.28 , 6.14 , 7.46 , 9.24 ] ) x_4 = np . array ( [ 0 , 20.00 , 22.51 , 25.10 , 27.69 , 30.11 , 32.45 , 35.04 , 37.88 , 40.14 , 42.65 , 45.07 , 47.33 , 50.08 ] ) y_4 = np . array ( [ 2.49 , 2.49 , 2.64 , 2.87 , 3.34 , 3.81 , 4.43 , 5.20 , 6.29 , 7.38 , 9.01 , 11.11 , 14.29 , 19.34 ] ) x_10 = np . array ( [ 0 , 20.00 , 22.50 , 25.08 , 28.00 , 30.00 , 32.50 , 34.92 , 37.50 , 40.17 , 42.42 , 45.00 , 47.17 , 50.08 ] ) y_10 = np . array ( [ 3.27 , 3.27 , 3.74 , 4.44 , 5.37 , 6.07 , 7.16 , 8.33 , 10.04 , 12.30 , 15.95 , 21.17 , 27.47 , 40.00 ] ) x_int = sl_0 . phi if sl_0 . phi < 1 : fd . ks = 0 else : if q1_q0 == 0 : fd . ks = np . interp ( x_int , x_0 , y_0 ) elif q1_q0 == 0.2 : fd . ks = np . interp ( x_int , x_2 , y_2 ) elif q1_q0 == 0.4 : fd . ks = np . interp ( x_int , x_4 , y_4 ) elif q1_q0 == 1.0 : fd . ks = np . interp ( x_int , x_10 , y_10 ) elif 0 < q1_q0 < 0.2 : ks_1 = np . interp ( x_int , x_0 , y_0 ) ks_2 = np . interp ( x_int , x_2 , y_2 ) fd . ks = ( ( ( ks_2 - ks_1 ) * q1_q0 ) / 0.2 ) + ks_1 elif 0.2 < q1_q0 < 0.4 : ks_1 = np . interp ( x_int , x_2 , y_2 ) ks_2 = np . interp ( x_int , x_4 , y_4 ) fd . ks = ( ( ( ks_2 - ks_1 ) * ( q1_q0 - 0.2 ) ) / 0.2 ) + ks_1 elif 0.4 < q1_q0 < 1.0 : ks_1 = np . interp ( x_int , x_4 , y_4 ) ks_2 = np . interp ( x_int , x_10 , y_10 ) fd . ks = ( ( ( ks_2 - ks_1 ) * ( q1_q0 - 0.4 ) ) / 0.6 ) + ks_1 else : raise DesignError ( "Cannot compute 'ks', bearing ratio out-of-range (q1_q0 = %.3f) required: 0-1." % q1_q0 ) # qu q_ult4 = ( r * ( 2 * ca * ( h0 - fd . depth ) / fd . width ) * a ) q_ult5_ks = q_ult5 * fd . ks q_ult = q_b + q_ult4 + q_ult5_ks - q_ult6 q_t1 = ( sl_0 . cohesion * sl_0 . nc_factor_0 * sl_0 . s_c_0 ) q_t = q_t1 + q_t2 + q_t3 if q_ult > q_t : if h0 > fd . width / 2 : fd . q_ult = q_t else : vert_eff_stress_interface = sp . vertical_effective_stress ( h0 ) vert_eff_stress_lowest = sp . vertical_effective_stress ( fd . width + fd . depth ) average_eff_stress = ( vert_eff_stress_interface + vert_eff_stress_lowest ) / 2 c_2_eff = sl_1 . cohesion + average_eff_stress * np . tan ( np . radians ( sl_1 . phi ) ) if sl_0 . cohesion > c_2_eff : fd . q_ult = q_t else : # vd = {} # vd[1] =[1, 1, 1, 1, 1] # vd[0.667] = [1, 1.033, 1.064, 1.088, 1.109] # vd[0.5] = [1, 1.056, 1.107, 1.152, 1.193] # vd[0.333] = [1, 1.088, 1.167, 1.241, 1.311] # vd[0.25] = [1, 1.107, 1.208, 1.302, 1.389] # vd[0.2] = [1, 1.121, 1.235, 1.342, 1.444] # vd[0.1] = [1, 1.154, 1.302, 1.446, 1.584] h_over_b = ( h0 - fd . depth ) / fd . width c1_over_c2 = sl_0 . cohesion / c_2_eff c_1_over_c_2 = [ 0.1 , 0.2 , 0.25 , 0.333 , 0.5 , 0.667 , 1. ] m_1 = [ 1.584 , 1.444 , 1.389 , 1.311 , 1.193 , 1.109 , 1. ] m_125 = [ 1.446 , 1.342 , 1.302 , 1.241 , 1.152 , 1.088 , 1. ] m_167 = [ 1.302 , 1.235 , 1.208 , 1.167 , 1.107 , 1.064 , 1. ] m_25 = [ 1.154 , 1.121 , 1.107 , 1.088 , 1.056 , 1.033 , 1. ] m_5 = [ 1 , 1 , 1 , 1 , 1 , 1 , 1 ] if h_over_b == 0.1 : m = np . interp ( c1_over_c2 , c_1_over_c_2 , m_1 ) elif h_over_b == 0.125 : m = np . interp ( c1_over_c2 , c_1_over_c_2 , m_125 ) elif h_over_b == 0.167 : m = np . interp ( c1_over_c2 , c_1_over_c_2 , m_167 ) elif h_over_b == 0.250 : m = np . interp ( c1_over_c2 , c_1_over_c_2 , m_25 ) elif h_over_b >= 0.5 : m = np . interp ( c1_over_c2 , c_1_over_c_2 , m_5 ) elif 0.1 < h_over_b < 0.125 : m_a = np . interp ( c1_over_c2 , c_1_over_c_2 , m_1 ) m_b = np . interp ( c1_over_c2 , c_1_over_c_2 , m_125 ) m = np . interp ( h_over_b , [ 0.1 , 0.125 ] , [ m_a , m_b ] ) elif 0.125 < h_over_b < 0.167 : m_a = np . interp ( c1_over_c2 , c_1_over_c_2 , m_125 ) m_b = np . interp ( c1_over_c2 , c_1_over_c_2 , m_167 ) m = np . interp ( h_over_b , [ 0.125 , 0.167 ] , [ m_a , m_b ] ) elif 0.167 < h_over_b < 0.25 : m_a = np . interp ( c1_over_c2 , c_1_over_c_2 , m_167 ) m_b = np . interp ( c1_over_c2 , c_1_over_c_2 , m_25 ) m = np . interp ( h_over_b , [ 0.167 , 0.250 ] , [ m_a , m_b ] ) elif 0.25 < h_over_b < 0.5 : m_a = np . interp ( c1_over_c2 , c_1_over_c_2 , m_25 ) m_b = np . interp ( c1_over_c2 , c_1_over_c_2 , m_5 ) m = np . interp ( h_over_b , [ 0.250 , 0.500 ] , [ m_a , m_b ] ) fd . q_ult = ( sl_0 . cohesion * m * sl_0 . nc_factor_0 ) + ( unit_eff_weight_0_at_fd_depth * fd . depth ) else : fd . q_ult = q_ult return fd . q_ult
Calculates the two-layered foundation capacity according Meyerhof and Hanna (1978)

    :param sp: Soil profile object
    :param fd: Foundation object
    :param wtl: water table level
    :param verbose: verbosity
    :return: ultimate bearing stress
Calculates the two - layered foundation capacity according Meyerhof and Hanna ( 1978 )