Fundacoes

fundacoes.py

import math
import numpy as np
import scipy.optimize as optimize
import matplotlib.pyplot as plt
import sys
from tabulate import tabulate


def _round(n):
    if n - math.floor(n) < 0.5:
        return math.floor(n)
    return math.ceil(n)


tabela_peso_especifico_argila = [
    # [menor, maior, valor]
    [False, 3,      13],
    [2,     6,      15],
    [5,     11,     17],
    [10,    20,     19],
    [19,    False,  21],
]

tabela_peso_espesifico_areia = [
    # [menor, maior, umida, seca, saturada]
    [False, 5,      18, 16, 19],
    [4,     9,      18, 16, 19],
    [8,     19,     19, 17, 20],
    [18,    41,     20, 18, 21],
    [40,    False,  20, 18, 21],
]

tabela_coesao = [
    # [menor, maior, menor, maior]
    [False, 2,      False,  10],
    [2,     4,      10,     25],
    [4,     8,      25,     50],
    [8,     15,     50,     100],
    [15,    30,     100,    200],
    [30,    False,  200,    False],
]


def angulo_de_atrito_teixeira(nspt):
    """
        segundo Teixeira (1996)
        citado por Cintra et al (2003)
    """
    return math.sqrt(20 * nspt) + 15

def peso_espesifico_argila(nspt):
    for t in tabela_peso_especifico_argila:
        if t[0] < nspt <t[1] or t[1] == False and t[0] < nspt:
            return t[2]

def peso_espesifico_areia(nspt, condicao=0):
    """
        condicao:
            0: umida
            1: seca
            2: saturada
    """
    for t in tabela_peso_espesifico_areia:
        if t[0] < nspt <t[1] or t[1] == False and t[0] < nspt:
            return t[2 + condicao]

def coesao(nspt):
    for t in tabela_coesao:
        if t[0] <= nspt < t[1] or t[1] == False and t[0] < nspt or \
        t[0] == 15 and t[0] <= nspt <= t[1]:
            if t[1] == False:
                return nspt * t[2] / t[0]
            elif t[0] == False:
                return nspt * t[3] / t[1]
            else:
                return t[3] - (t[1] - nspt) * (t[3] - t[2]) \
                / (t[1] - t[0])

def parametros_do_solo_spt(spt, condicao='umida'):
    """
        -> A FAZER: nivel agua

        spt é o resultado do perfil de solo,
        seguindo o seguinte formato:
        [
            [nspt, solo, h],
            ...
        ]
        onde:
            nspt:   -> valor medio de golpes ex: 6  ou
                    -> array de golpes ex: [5,6,5]

            solo:   -> 0 solo predominantemente arenoso ou
                    -> 1 solo predominantemente argiloso

            h:      -> profundidade da camada em metros

        exemplo:
            spt = [
                [6, 0, 4],
                [12, 0, 4],
                [17, 1, 5],
                [21, 1, 3]
            ]

        retorna Array com um subArray para cada camada no formato:
            [
                [phi, g, c, h],
                [...]
            ]
        onde:
            phi:    -> angulo de atrito por teixeira
            g:      -> peso especifico
            c:      -> coesao
            h:      -> altura da camada
    """
    parametros = []
    for camada in spt:
        if isinstance(camada[0], list):
            camada[0] = _round(np.mean(camada[0]))
        else:
            camada[0] = _round(camada[0])
        phi = angulo_de_atrito_teixeira(camada[0])
        c = coesao(camada[0])
        if camada[1] == 1:
            g = peso_espesifico_argila(camada[0])
        elif camada[1] == 0:
            g = peso_espesifico_areia(camada[0])
            c = c * 0.5
        h = camada[2]
        parametros.append([phi, g, c, h])
    return parametros

tabela_fatores_de_forma = [
    [1.1, 0.9, 1 ], #retangular
    [1.3, 0.8, 1 ], #quadrada
    [1.0, 1.0, 1 ], #corrida
    [1.3, 0.6, 1 ], #circular
]

# tabela_fatores_de_carga = np.array([
#     [0, 5,   10,  15,  20,  25,  30,  35,  40, ],        # angulos
#     [5.14,  5.5,  7.5,  10, 14.5, 20.5, 30, 45, 90,],    # nc
#     [0, 0, 0.05, 0.5, 2, 4, 18, 40, 85],                 # ng
#     [1, 1.1, 2.5, 4, 7, 10.5, 17, 32.5, 60],             # nq
# ])

tabela_fatores_de_carga = np.array([
    [0, 5,   10,  15,  20,  25,  30,  35,  40],              # angulos
    [5.7,  7.3,  9.6,  12.9, 17.7, 25.1, 37.2, 57.8, 95.7],  # nc
    [0, 0.5, 1.2, 2.5, 5, 9.7, 19.7, 42.4, 100.4],           # ng
    [1, 1.6, 2.7, 4.4, 7.4, 12.7, 22.5, 41.4, 81.3],         # nq
])

def terzaghi_fator_carga(phi, index):
    """
        np.interp(x, xp, fp)
        Retorna o interpolante linear unidimensional por
        partes para uma função com determinados pontos de
        dados discretos (xp, fp), avaliados em x.
    """
    return np.interp(phi, tabela_fatores_de_carga[0],
        tabela_fatores_de_carga[index])

def terzaghi(c, b, h, g, phi, q, g2=None, phi2=None, forma=0, full=False):
    """
        c = coesao camada inferior
        b = menor dimensao
        q = tensao efetiva
        g = peso especifico do solo lateral
        g2 = peso especifico do solo abaixo da cota de assentamento
        phi = angulo de atrito do solo lateral
        phi2 = angulo de atrito do solo abaixo da cota de assentamento
        forma = [0,1,2,3] -> [retangular,quadrada,corrida,circular]
    """
    # se não foi dado phi2 ou g2, considera o mesmo solo
    if phi2 == None:
        phi2 = phi
    if g2 == None:
        g2 = g

    sc, sg, sq = tabela_fatores_de_forma[forma]

    nc = terzaghi_fator_carga(phi, 1)
    ng = terzaghi_fator_carga(phi, 2)
    nq = terzaghi_fator_carga(phi2, 3)

    if full:
        coesao = c*nc*sc
        atrito = b/2*g2*ng*sg
        sobrecarga = q*nq*sq
        return {
            "nc": nc,
            "ng": ng,
            "nq": nq,
            "sc": sc,
            "sg": sg,
            "sq": sq,
            "coesao": coesao,
            "atrito": atrito,
            "sobrecarga": sobrecarga,
            "total": coesao + atrito + sobrecarga
        }

    return c*nc*sc + b/2*g2*ng*sg + q*nq*sq

# def fundacao_rasa_terzaghi(h, b, parametros, forma=0, margem=0.0):
#     """
#         h: profundidade da fundcao
#         b: menor dimensao
#         parametros: parametros do solo, obtidos em: parametros_do_solo_spt()
#         margem: distancia minima para conciderar a camada inferior,
#                 impede a otimizacao de considerar [x].99 diferente de
#                 [x+1].00 na profundidade, para adotar fatores de carga
#                 que forcem melhores parametros (de forma errada)
#
#         Distingue entre quais são os parametros a serem
#         considerados para o solo inferiro e lateral à fundação
#
#         retorna:
#             terzaghi(...) -> chamada com os parametros corretos
#     """
#     profundidade = 0
#     ok = False
#     for p in parametros:
#         profundidade += p[3]
#         if profundidade - p[3] < h :
#             phi = p[0]
#             g = p[1]
#             _c = p[2]
#         if ok == False and profundidade >= h + margem:
#             phi2 = p[0]
#             g2 = p[1]
#             c = p[2]
#             ok = True
#     if ok == False:
#         c = _c
#         phi2 = phi
#         g2 = g
#
#     return terzaghi(c, b, h, g, phi, g2, phi2, forma)


def _delta(v, carga, _q, phi, phi2, g, g2, c, z, forma=1, fs=3, _h=None, _b=None,
    l=None, inverte=False):
    """
        retorna a diferenca entra t_adm e t_ap
    """
    if inverte == True:
        b = v[0]
        h = _h
    else:
        h = v[0]
        b = _b
    if forma == 1: # quadrada
        a = b * b
    elif forma == 3: # circular
        a = math.pi * (b/2)**2
    else:
        a = b * l
    q = (h - z) * g + _q
    t_ap = carga / a
    t_adm = terzaghi(c, b, h, g, phi, q, g2, phi2, forma) / fs
    return t_adm - t_ap

def fundacao_rasa_terzaghi_otimizada(carga, spt, b=None, h=None, l=None,
    forma=1, fs=3, b_min=0.25, b_max=2, h_max=4):
    """
        otimiza para achar a menor profundidade ou a menor dimensao,
        baseado nas entradas de dado.

        se entrar com b, vai achar a menor h
        se entrar com h, vai achat a menor b

        se não encontrar uma otimizacao, retorna None
        se a altura for negativa, retorna 0, significa
        que para dado b, poderia estar na superficie
    """

    _params = parametros_do_solo_spt(spt)
    params = []

    profundidade = 0

    for i, p in enumerate(_params):
        if i > 0:
            _q = params[i-1]["_q"] + _params[i-1][1]*_params[i-1][3]
        else:
            _q = 0
        params.append({
            "phi": p[0],
            "g": p[1],
            "c": p[2],
            "phi2": p[0],
            "g2": p[1],
            "h": (profundidade, p[3] + profundidade),
            "_q": _q,
            "z": p[3]
        })
        profundidade += p[3]

    res = None

    z = 0

    for p in params:
        if h ==  None:
            args = (carga, p["_q"], p["phi"], p["phi2"], p["g"], p["g2"],
                p["c"], z, forma, fs, False, b, l)
            root = optimize.root(_delta, [1], args)
            if res == None and root.success == True:
                if p["h"][0] <= root.x <= p["h"][1] and root.x <= h_max:
                    res = root.x[0]
                    break
                elif root.x[0] <= 0 :
                    root = 0
                    break
        elif p["h"][0] <= h < p["h"][1]:
            args = (carga, p["_q"], p["phi"], p["phi2"], p["g"], p["g2"],
                p["c"], z, forma, fs, h, False, l, True)
            root = optimize.root(_delta, [1], args)
            if res == None and root.success == True:
                res = root.x[0]
                break
        z += p['z']

    return res


def comportamento_fundacao_rasa_terzaghi(carga, spt, show=False, b_min=0.25, b_max=4,
    formas=[1,3], h_min=0, h_max=4, **kwargs):
    """
        Faz o grafico das dimensoes otimizadas de uma sapata,
        dado perfil de solo e carga
    """
    fig, ax = plt.subplots()

    # h = np.arange(0, h_max+0.01, 0.01)
    # bs = []
    # for forma in formas:
    #     b = []
    #     for i in h:
    #         b.append(fundacao_rasa_terzaghi_otimizada(carga, spt, h=i,
    #         **kwargs, forma=forma, h_max=h_max))
    #     bs.append(b)
    #
    # fig, ax = plt.subplots()
    # ax.set(xlabel='Dimensão mínima (m)', ylabel='Profundidade mínima (m)',
    #        title='Dimensões otimizadas, para dado spt e carga de ' + str(carga) + 'kPa')
    # ax.grid()

    b = np.arange(b_min, b_max+0.01, 0.01)
    hs = []
    for forma in formas:
        h = []
        for i in b:
            h.append(fundacao_rasa_terzaghi_otimizada(carga, spt, b=i,
            **kwargs, forma=forma, h_max=h_max))
        hs.append(h)

    fig, ax = plt.subplots()
    ax.set(xlabel='Dimensão mínima (m)', ylabel='Profundidade mínima (m)',
           title=f'Dimensões otimizadas, para dado spt e carga de {carga} kPa')
    ax.grid()

    for i, h in enumerate(hs):
        label = None
        if formas[i] == 0:
            label = 'retangular'
        elif formas[i] == 1:
            label = 'quadrada'
        elif formas[i] == 2:
            label = 'corrida'
        elif formas[i] == 3:
            label = 'circular'
        plt.plot(b, h, label=label)

    plt.legend()

    fig.savefig("comportamento_fundacao_rasa_terzaghi.png")

    if show == True:
        plt.show()

def vista_superior(b, p=0.3, show=False):
    """
        A fazer: retangular e quadrada
    """
    p = p/2
    fig, ax = plt.subplots()
    ax.set(title='Vista Superior')

    verts = [
        [(0, 0),(b, 0), (b, b), (0, b), (0,0)],
        [(-p+b/2, -p+b/2),(p+b/2, -p+b/2),(p+b/2,p+b/2),(-p+b/2, p+b/2),(-p+b/2, -p+b/2)],
        [(-p+b/2, -p+b/2),(0, 0)],
        [(p+b/2, p+b/2),(b, b)],
        [(-p+b/2, p+b/2),(0, b)],
        [(p+b/2, -p+b/2),(b, 0)],
    ]
    for v in verts:
        xs, ys = zip(*v)
        ax.plot(xs, ys, '-', color='black')

    verts = [(0, -p),(b, -p)]
    xs, ys = zip(*verts)
    ax.plot(xs, ys, '|-', color='black')
    ax.text(b/2, -p*1.5, f'{b} m')

    verts = [(-p, 0),(-p, b)]
    xs, ys = zip(*verts)
    ax.plot(xs, ys, '_-', color='black')
    ax.text(-p*2, b/2, f'{b} m', rotation=-90)

    plt.axis('off')
    plt.gca().set_aspect('equal', adjustable='box')

    fig.savefig(f"vista_superior_{b}.png")

    if show == True:
        plt.show()

def vista_frontal(b, h, p=0.3, h2=0.5, show=False):
    """
        A fazer: retangular e quadrada
    """
    p = p/2
    fig, ax = plt.subplots()
    ax.set(title='Vista Frontal')

    verts = [
        [(-b, 0), (b, 0)],
        [(-p, 0),(-p, -h+h2)],
        [(p, 0),(p, -h+h2)],
        [(-p, -h+h2),(-b/2, -h+h2), (-b/2, -h), (b/2, -h), (b/2, -h+h2), (p, -h+h2)],
    ]
    for v in verts:
        xs, ys = zip(*v)
        ax.plot(xs, ys, '-', color='black')

    verts = [(-b/2, -h-p),(b/2, -h-p)]
    xs, ys = zip(*verts)
    ax.plot(xs, ys, '|-', color='black')
    ax.text(0, -h-p*2, str(b)+' m')

    verts = [(-b/2-p, 0),(-b/2-p, -h)]
    xs, ys = zip(*verts)
    ax.plot(xs, ys, '_-', color='black')
    ax.text(-b/2-p*2, -h/2, str(h)+' m', rotation=-90)

    plt.axis('off')
    plt.gca().set_aspect('equal', adjustable='box')

    fig.savefig(f"vista_frontal_{b}.png")

    if show == True:
        plt.show()

def test():
    def peso_espesifico_argila_test():
        tests = [
            (1,13),
            (4,15),
            (5,15),
            (6,17),
            (17,19),
            (20, 21)
        ]
        for t in tests:
            assert peso_espesifico_argila(t[0]) == t[1]

    def peso_espesifico_areia_test():
        tests = [
            (1,0,18),
            (1,1,16),
            (1,2,19),
            (10,0,19),
            (40,0,20),
            (41,0,20),
            (41,2,21)
        ]
        for t in tests:
            assert peso_espesifico_areia(t[0], t[1]) == t[2]

    def coesao_test():
        tests = [
            (6,37),
            (12,78),
            (17,113),
            (21,140)
        ]
        for t in tests:
            assert math.floor(coesao(t[0])) == t[1]

    def angulo_de_atrito_teixeira_test():
        tests = [(6,25), (12,30), (17,33), (21,35)]
        for t in tests:
            assert math.floor(angulo_de_atrito_teixeira(t[0])) == t[1]

    def parametros_do_solo_spt_test():
        spt = [
            [6, 0, 4],
            [12, 0, 4],
            [17, 1, 5],
            [21, 1, 3],
        ]
        esperado = [
            [26, 18, 18, 4],
            [31, 19, 39, 4],
            [34, 19, 113, 5],
            [36, 21, 140, 3],
        ]
        resultado = parametros_do_solo_spt(spt)
        for index, camada in enumerate(resultado):
            assert math.ceil(camada[0]) == esperado[index][0]
            assert camada[1] == esperado[index][1]
            assert math.floor(camada[2]) == esperado[index][2]
            assert camada[3] == esperado[index][3]

    def terzaghi_fator_carga_test():
        tests = [
            (10, 1, 9.6),
            (10, 2, 1.2),
            (10, 3, 2.7),
            (12, 1, 10.92),
        ]
        for t in tests:
            assert terzaghi_fator_carga(t[0], t[1]) == t[2]


    def terzaghi_test():
        kwargs = {'c': 1, 'b': 1, 'h': 2, 'g': 20, 'phi': 20}
        esperado = 360
        assert math.floor(terzaghi(**kwargs)) == esperado


    def fundacao_rasa_terzaghi_test():
        """
            spt = [
                [11, 0, 3],  -> g=19 , phi=29.83 , c=35.71 , h=3
                [21, 0, 10], -> g=20 , phi=35.49 , c=70.00 ,  h=10
            ]
            h = 3
            b = 1
            para profundidade de 3m, como esta na interface entre
            duas camadas, sera diferenciada os parametros laterais
            e inferiroes da seguinte forma:
                parcela coesao: c=70 (inferior)
                parcela atrito: g=20 (inferior)
                parcela sobrecarga: g=19, h=3 (lateral)
            os fatores de carga:
                calculados da camada lateral -> nc = 37, ng = 19
                calculado da camada inferior -> nq = 45
            forma retangular => sc=1.1 , sg=0.9, sq=1

            logo capacidade de carga =
                70*37*1.1 = 2849
                1/2*20*19*0.9 = 171
                19*3*45*1 = 2565
                = 5585 kpa
                 +- 10 kpa de erro devido olho nu no abaco
                 e erredondamentos nos angulos e fatores
        """
        spt = [
            [11, 0, 3],
            [21, 0, 10],
        ]
        h = 3
        b = 1
        parametros = parametros_do_solo_spt(spt)
        esperado = 5585
        f = fundacao_rasa_terzaghi(h, b, parametros, forma=0, margem=1)
        # computa erro do calculo feito a mao, devido arredondamentos
        # e olho nu no abaco de fatores de carga, erro +- 10 kpa
        assert esperado - 10 < f < esperado + 10

    peso_espesifico_argila_test()
    peso_espesifico_areia_test()
    coesao_test()
    angulo_de_atrito_teixeira_test()
    parametros_do_solo_spt_test()
    terzaghi_fator_carga_test()
    terzaghi_test()
    # fundacao_rasa_terzaghi_test()
    print('Todos os testes completos a bem sucedidos')

spt = [
    #[nspt, 0-1, profundidade]
    [4,1,2],
    [9,0,3],
    [16,0,4],
    [22,0,6]
]

if __name__ == "__main__":
    # spt = [[4,1,2],[9,0,3],[16,0,4],[22,0,6]]
    # carga = 755
    try:
        comando = sys.argv[1]
    except:
        raise SyntaxError("Entre com um comando.")

    args = []
    if len(sys.argv) > 2:
        for i, a in enumerate(sys.argv[2:]):
            if "=" not in a:
                args.append(eval(a))
        kwargs = dict(arg.split('=') for arg in sys.argv[2+len(args):])
        for i in kwargs:
            kwargs[i]=eval(kwargs[i])

    if "-o" == comando:
        print(fundacao_rasa_terzaghi_otimizada(*args, **kwargs))
    elif "-c" == comando:
        comportamento_fundacao_rasa_terzaghi(*args, **kwargs)
    elif "-vs" == comando:
        vista_superior(*args, **kwargs)
    elif "-vf" == comando:
        vista_frontal(*args, **kwargs)
    elif "-ps" == comando:
        params = parametros_do_solo_spt(*args, **kwargs)
        headers = ['camada','angulo (°)','peso esp. (km/m³)','coesão (kPa)','altura (m)']
        linhas = []
        for p in params:
            p[0] = round(p[0],2)
            p[2] = round(p[2],2)
            linhas.append([*p])
        print('\nParametros do solo:')
        print(tabulate(linhas, headers=headers, showindex="always", floatfmt=".2f"))
    elif "-pt" == comando:
        # c=25 b=1.5 h=0.95 g=15 phi=23.94 q=14.25 forma=1
        t = terzaghi(*args, **kwargs, full=True)
        linhas = []
        for item in t:
            linhas.append([item, t[item]])
        print('\nTerzaghi:')
        print(tabulate(linhas, numalign="right", floatfmt=".2f"))