tensoes_talude.py

import numpy as np
from collections import namedtuple
import matplotlib.pyplot as plt


# Dados do problema:
gama_a = 20
gama_n = 23
h = 12
sup_base = 10
left_angle = 1/1
right_angle = 2/1
phi = 35
points = [
    (0,5),
    (12,5),
    (17,5),
    (22,5),
    (46,5),
]


# Parte 0: inicializacao e definicao de utilidades
Case = namedtuple('Case', 'case variant')
Tensions = namedtuple('Tensions', 'delta v h')

def get_x(angle, h):
    return angle * h

def get_alfas(b1, b2, z):
    alfa2 = np.arctan(b1/z)
    return np.arctan((b1 + b2)/z) - alfa2, alfa2

def get_delta(gama_a, h, z, b1, b2):
    alfa1, alfa2 =  get_alfas(b1, b2, z)
    return gama_a * h / np.pi * ((b1 + b2) / b2 * (alfa1 + alfa2) - b1 * alfa2/ b2)

def get_total_v(gama_n, z, delta):
    return gama_n * z + delta

def get_k0(phi):
    return 1 - np.sin(np.deg2rad(phi))

def get_total_h(total_v, phi):
    return total_v * get_k0(phi)

def evaluate_superposition(px, left, sup_base, right):
    """
    -> retorna:
        um caso de superposicao de acordo com a
        coordenada x e as dimencoes do talude.

    -> caveats:
        o talude deve ser definido considerando inicio de sua base px == 0
        posicoes anteriores ao talude devem ser negativas

    -> casos:
        0: ponto em baixo da extremidade mais a esquerda ou direita do talude
        1: ponto em baixo da interface triangulo/retangulo
        2: ponto em qualquer lugar em baixo do retangulo
        3: ponto nao em baixo do talude
        4: ponto abaixo da superficie inclinada do talude
    """
    if px == 0:
        return Case(0, 0)
    elif px == left + sup_base + right:
        return Case(0, 1)

    elif px == left:
        return Case(1, 0)
    elif px == left + sup_base:
        return Case(1, 1)

    elif px > left and px < left + sup_base:
        return Case(2, 0)

    elif px < 0:
        return Case(3, 0)
    elif px > left + sup_base + right:
        return Case(3, 1)

    elif px > 0 and px < left:
        return Case(4, 0)
    elif px > left + sup_base and px < left + sup_base + right:
        return Case(4, 1)

def get_tensions(x, z, left, sup_base, right, gama_a, gama_n, phi):
    """
    -> retorna:
        Tensions(delta=..., v=..., h=...)
        delta = aumento de tensão vertical
        v = tensão vertical total
        h = tensão horizontal total

    -> escopo:
        pontos fora e a baixo de qualquer parte do talude

    -> caveats:
        ver o caveats da funcao 'evaluate_superposition'

    -> limitacoes:
        nao foi testado para pontos muito distantes do talude,
        ou seja, casos em qua na verdade a tensao seria zero,
        pode resultar em valores inesperados
    """
    sp = evaluate_superposition(x, left, sup_base, right)
    if sp.case == 0:
        t2_b1 = 0
        if sp.variant == 0:
            t1_b1 = left + sup_base
            t1_b2 = right
            t2_b2 = left
        else:
            t1_b1 = right + sup_base
            t1_b2 = left
            t2_b2 = right
    elif sp.case == 1:
        t1_b1 = sup_base
        t2_b1 = 0
        if sp.variant == 0:
            t1_b2 = right
            t2_b2 = left
        else:
            t1_b2 = left
            t2_b2 = right
    elif sp.case == 2:
        t1_b1 = x - left
        t1_b2 = left
        t2_b1 = left + sup_base - x
        t2_b2 = right
    elif sp.case == 3:
        if sp.variant == 0:
            t1_b1 = np.absolute(x) + left + sup_base
            t1_b2 = right
            t2_b1 = np.absolute(x)
            t2_b2 = left
        else:
            t1_b1 = x - left
            t1_b2 = left
            t2_b1 = x - left - sup_base - right
            t2_b2 = right
    elif sp.case == 4:
        t2_b1 = 0
        if sp.variant == 0:
            t1_b1 = left - x + sup_base
            t1_b2 = right
            t2_b2 = x
            t2_h = t2_b2 / left * h
            t3_b2 = left - x
            t3_h = t3_b2 / left * h
        else:
            t1_b1 = x - left
            t1_b2 = left
            t2_b2 = left + sup_base + right - x
            t2_h = t2_b2 / right * h
            t3_b2 = x - left - sup_base
            t3_h = t3_b2 / right * h
    if sp.case != 4:
        t1_delta = get_delta(gama_a, h, z, t1_b1, t1_b2)
        t2_delta = get_delta(gama_a, h, z, t2_b1, t2_b2)
        if sp.case == 0 or sp.case == 3:
            delta = t1_delta - t2_delta
        else:
            delta = t1_delta + t2_delta
    else:
        t1_delta = get_delta(gama_a, h, z, t1_b1, t1_b2)
        t2_delta = get_delta(gama_a, t2_h, z, t2_b1, t2_b2)
        t3_delta = get_delta(gama_a, t3_h, z, 0, t3_b2)
        delta = t1_delta + t2_delta - t3_delta
    total_v = get_total_v(gama_n, z, delta)
    total_h = get_total_h(total_v, phi)
    return Tensions(delta, total_v, total_h)

left = get_x(left_angle, h)
right = get_x(right_angle, h)


# Graficos:
base = left + right + sup_base
Z = np.arange(1, base + 1, base / 20)
X = np.arange(-base, base * 2, 0.01)

fig, axs = plt.subplots(4, 1, figsize=(8.27, 11.69))
fig.subplots_adjust(hspace=0.8)
for ax in axs:
    ax.invert_yaxis()
    ax.set(xlabel='distancia (m)')

# grafico apenas dos pontos pedidos
x = []
z = []
file = []
for p in points:
    ts = get_tensions(p[0], p[1], left, sup_base, right, gama_a, gama_n, phi)
    z.append(ts.v)
    x.append(p[0])
    file.append((p[0], ts.v))
axs[0].plot(x, z)
file = np.asarray(file)
np.savetxt(
    'tensao_total_vertical_pontos_pedidos.csv',
    file,
    delimiter=';',
    fmt="%.2f",
    header="x;v"
)

# graficos generalizados
for z in Z:
    Y_delta = []
    Y_v = []
    Y_h = []
    for x in X :
        res = get_tensions(x, z, left, sup_base, right, gama_a, gama_n, phi)
        Y_delta.append(res.delta)
        Y_v.append(res.v)
        Y_h.append(res.h)
    axs[1].plot(X, Y_delta)
    axs[2].plot(X, Y_v)
    axs[3].plot(X, Y_h)

axs[0].set(ylabel='v (kpa)', title='Tensao vertical total')
axs[1].set(ylabel='delta (kpa)', title='Variacao de tensao vertical')
axs[2].set(ylabel='v (kpa)', title='Tensao vertical total')
axs[3].set(ylabel='h (kpa)', title='Tensao horizontal total')

plt.savefig('distribuicao_tensoes_talude.pdf')