【水利计算程序】python计算心墙土石坝有棱体排水（不透水地基）浸润线

# -*- coding:utf-8 -*-
# @Time     :2021-03-24 9:44

def Phreatic_line(x,m1=3,m2=2.5,m3=0.3,m4=0.3,m5=1.5,k1=0.2,k2=0.01,k3=0.5,B1=8,B2=26,H1=30,H2=4,d=6,L=50,hej=17.285):
    '''心墙土石坝有棱体排水（不透水地基）浸润线计算公式'''
    # 上游楔形体宽度m
    AL = k2 / k1 * ((H1 * m1 / (2 * m1 + 1) - m3 / (2 * m4 + 1)) + d)
    # 心墙顶部的计算宽度
    B11 = B1 + m3 / (2 * m3 + 1) * H1 + AL
    # 心墙底部计算宽度
    B22 = B2 - 2 * pow(m3, 2) / (2 * m3 + 1) * H1 + AL
    # 其值大于0.5时，按宽心墙面坝计算，小于0.5时按窄心墙坝计算
    # B22/H1=B22/H1
    # k/k1大于1000时ht=H2,K/K1小于1000时，ht稍大于H2,假设
    htj = H2
    # 单宽渗流量m3/d
    q = k2 * ((pow(H1 - htj, 2) / (
                (B22 - m4 * htj) + pow(pow(B22 - m4 * htj, 2) - pow(m4, 2) * pow(H1 - htj, 2), 0.5)) + (
                           H1 - htj) * htj / (B22 - 0.5 * m4 * htj)))
    # 棱体左下角上方水深
    h0 = 0.5 * (m5 * q / k3 + pow(pow(m5 * q / k3, 2) + 4 * pow(H2, 2), 0.5))
    # k/k1大于1000时ht=H2,K/K1小于1000时，ht稍大于H2
    ht = pow(pow(m3 * q / k3, 2) + 2 * L * q / k3 + pow(h0, 2) - m1 * q / k3, 0.5)
    # 心墙前的渗流深度m
    hw = pow(pow(H1, 2) - 2 * q / k1 * AL, 0.5)

    # 心墙下游坡渗出点处水深
    he = ((2 * pow(m4 + 0.5, 2) * (hej - ht) + m4 * ht) * (m4 + 0.5)) / (
                2 * pow(m4 + 0.5, 2) * hej + m4 * ht) * q / k2 + ht
    # 以心墙后he点的坝底为y坐标原点，方向向上
    y = pow(x * (pow(hw, 2) - pow(he, 2)) / (B2 - m3 * hw - m4 * he) + pow(he, 2), 0.5)
    return y

if __name__=="__main__":
    import matplotlib.pyplot as plt
    x = [0.5, 1, 3, 5, 7, 9, 10, 11]
    # x=[110.61,105.26,125.38,102.04,101.60,103.99,94,0]
    y=list(map(Phreatic_line,x))
    #print(y)
    #计算结果
    # [18.00535017357371, 18.695870580834278, 21.2345789782722, 23.500619395957514, 25.56659694998201, 27.477675433194936,
    #  28.38500539046262, 29.264217310927336]
    plt.plot(x,y)
    plt.show()