mechan_rgr2/rgr.py

#!/usr/bin/env python3

# Copyright © 2019 Anton Tsyganenko <anton-tsyganenko@yandex.ru>
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

# Программа принимает на вход вариант в виде списков действующих сил,
# распределенных нагрузок и сосредоточенных моментов (вариант задается
# чуть ниже) и автоматически:
#
# - находит реакции опор;
#
# - из условия прочности определяет размеры прямоугольного сечения;
#
# - выбирает по сортаменту двутавр, удовлетворяющий условиям прочности и
# жесткости;
#
# - проверяет условие прочности по нормальным и касательным напряжениям
# для сечений в виде прямоугольника и двутавра;
#
# - считает коэффициент экономии;
#
# - строит эпюры распределения нормальных и касательных напряжений по
# высоте сечения для сечений в виде двутавра и прямоугольника;
#
# - строит эпюры поперечных сил, изгибающих моментов, углов поворота
# сечений и прогибов для сечения в виде двутавра
#
# - проверяет выполнение условия жесткости;
#
# - строит таблицу поперечных сил, изгибающих моментов, углов поворота
# сечений и прогибов для сечения в виде двутавра, выделяя в таблице
# строчки, где максимально значение Q, M или y.

forces, qs, Ms = \
    ([tuple(float(i) for i in f.split()) if f != "" else (0,0,0)\
    for f in input(msg+": ").split(";")]\
    for msg in ("Силы", "Распределенные нагрузки", "Моменты"))


Xa = float(input("Координата первой опоры/заделки: "))
Xb_text = input("Координата второй опоры (если есть): ")

Xmax = float(input("Длина балки: "))

import matplotlib
import matplotlib.pyplot as plt
import numpy as np
import math
import os

os.chdir(".") # Куда сохранять результат

# Введем параметры материалов, допустимие значения величин и и т.д. В
# программе в переменных сохраняются численные значения соответствующих
# величин в основных единицах СИ.

E = 2e11
sigma_dop = 160e6
tau_dop = 100e6
hb_rel = 2 # отношение высоты прямоугольного сечения к ширине


# Введем сортамент, по которому программа будет выбирать двутавры. Если
# ввести не все существующие варианты, программа будет выбирать из
# введенных.


sortament = [
#    №,      h,      b,       d,       t,       A,   Jx (Ix),      Wx,      Sx
( "10", 100e-3,  55e-3,  4.5e-3,  7.2e-3, 12.0e-4,    198e-8, 39.7e-6, 23.0e-6),
( "12", 120e-3,  64e-3,  4.8e-3,  7.3e-3, 14.7e-4,    350e-8, 58.4e-6, 33.7e-6),
( "14", 140e-3,  73e-3,  4.9e-3,  7.5e-3, 17.4e-4,    572e-8, 81.7e-6, 46.8e-6),
( "16", 160e-3,  81e-3,  5.0e-3,  7.8e-3, 20.2e-4,    873e-8,  109e-6, 62.3e-6),
( "18", 180e-3,  90e-3,  5.1e-3,  8.1e-3, 23.4e-4,   1290e-8,  143e-6, 81.4e-6),
("18a", 180e-3, 100e-3,  5.1e-3,  8.3e-3, 25.4e-4,   1430e-8,  159e-6, 89.8e-6),
( "20", 200e-3, 100e-3,  5.2e-3,  8.4e-3, 25.4e-4,   1840e-8,  184e-6,  104e-6),
("20a", 200e-3, 110e-3,  5.2e-3,  8.6e-3, 28.9e-4,   2030e-8,  203e-6,  114e-6),
( "22", 220e-3, 110e-3,  5.4e-3,  8.7e-3, 30.6e-4,   2550e-8,  232e-6,  131e-6),
("22a", 220e-3, 120e-3,  5.4e-3,  8.9e-3, 32.8e-4,   2790e-8,  254e-6,  143e-6),
( "24", 240e-3, 115e-3,  5.6e-3,  9.5e-3, 34.8e-4,   3460e-8,  289e-6,  163e-6),
("24а", 240e-3, 125e-3,  5.6e-3,  9.8e-3, 37.5e-4,   3800e-8,  317e-6,  178e-6),
( "27", 270e-3, 125e-3,  6.0e-3,  9.8e-3, 40.2e-4,   5010e-8,  371e-6,  210e-6),
("27а", 270e-3, 135e-3,  6.0e-3, 10.2e-3, 43.2e-4,   5500e-8,  407e-6,  229e-6),
( "30", 300e-3, 135e-3,  6.5e-3, 10.2e-3, 46.5e-4,   7080e-8,  472e-6,  268e-6),
("30а", 300e-3, 145e-3,  6.5e-3, 10.7e-3, 49.9e-4,   7780e-8,  518e-6,  292e-6),
( "33", 330e-3, 140e-3,  7.0e-3, 11.2e-3, 53.8e-4,   9840e-8,  597e-6,  339e-6),
( "36", 360e-3, 145e-3,  7.5e-3, 12.3e-3, 61.9e-4,  13380e-8,  743e-6,  423e-6),
( "40", 400e-3, 155e-3,  8.0e-3, 13.0e-3, 71.4e-4,  18930e-8,  947e-6,  540e-6),
( "45", 450e-3, 160e-3,  8.6e-3, 14.2e-3, 83.0e-4,  27450e-8, 1220e-6,  699e-6),
( "50", 500e-3, 170e-3,  9.5e-3, 15.2e-3, 97.8e-4,  39290e-8, 1570e-6,  905e-6),
( "55", 550e-3, 180e-3, 10.3e-3, 16.5e-3,  114e-4,  55150e-8, 2000e-6, 1150e-6),
( "60", 600e-3, 190e-3, 11.1e-3, 17.8e-3,  132e-4,  75450e-8, 2510e-6, 1450e-6),
( "65", 650e-3, 200e-3, 12.0e-3, 19.2e-3,  153e-4, 101400e-8, 3120e-6, 1800e-6),
( "70", 700e-3, 210e-3, 13.0e-3, 20.8e-3,  176e-4, 134600e-8, 3840e-6, 2230e-6),
("70a", 700e-3, 210e-3, 15.0e-3, 24.0e-3,  202e-4, 152700e-8, 4360e-6, 2550e-6),
("70б", 700e-3, 210e-3, 17.5e-3, 28.2e-3,  234e-4, 175370e-8, 5010e-6, 2940e-6),
]

if Xb_text:
    Xb = float(Xb_text)
    vartype = "double"
    Ydop = Xmax/300
else:
    vartype = "hard"
    Ydop = Xmax/200


print("Допустимый прогиб: {:g} м".format(Ydop))


# Определим функции, вычисляющие поперечные силы и изгибающие моменты в
# заданной точке

def q(x):
    result = 0
    for force in forces:
        if force[1]<x:
            result += force[0]

    for q in qs:
        if q[1]<x:
            result += q[0]*(x-q[1])
        if q[2]<x:
            result -= q[0]*(x-q[2])

    return result

def m(x):
    result = 0

    for q in qs:
        if q[1]<x:
            result += q[0]*(x-q[1])**2/2
        if q[2]<x:
            result -= q[0]*(x-q[2])**2/2

    for force in forces:
        if force[1]<x:
            result += force[0]*(x-force[1])

    for moment in Ms:
        if moment[1]<x:
            result += moment[0]

    return result


# Найдем реакции опор исходя из равенства нулю суммы всех сил и моментов
# сил, и выбрав точку с x=Xa.

SF = sum(f[0] for f in forces)\
   + sum(q[0]*(q[2]-q[1]) for q in qs)

SM = - sum(m[0] for m in Ms)\
     + sum(f[0]*(f[1]-Xa) for f in forces)\
     + sum(q[0]*(q[2]-q[1])*((q[2]+q[1])/2-Xa) for q in qs)

if vartype == "double":
    Ra = -(SM - Xb*SF)/(Xa - Xb)
    Rb = -(Xa*SF - SM)/(Xa - Xb)
    forces.append((Rb, Xb))
    forces.append((Ra, Xa))
    print("Ra = {:g} Н; Rb = {:g} Н".format(Ra, Rb))
elif vartype == "hard":
    Ra = -SF
    Ma = SM
    forces.append((Ra, Xa))
    Ms.append((Ma, Xa))
    print("Ra = {:g} Н; Ma = {:g} Н*м".format(Ra, Ma))


# Для расчета углов поворота сечений и прогибов, сначала определим
# функцию `EIy_rel`, вычисляющую условный прогиб, принимая начальную
# координату и угол поворота сечения в начале за 0, затем вычислим угол
# поворота сечения в начале (с точностью до множителей, зависящих только
# от материала и параметров сечения: E и I), после чего определим
# функцию, вычисляющую EIθ (`EIt`) и EIy


def EIy_rel(x):
    result = 0

    for q in qs:
        if q[1]<x:
            result += q[0]*(x-q[1])**4/24
        if q[2]<x:
            result -= q[0]*(x-q[2])**4/24

    for force in forces:
        if force[1]<x:
            result += force[0]*(x-force[1])**3/6

    for moment in Ms:
        if moment[1]<x:
            result += moment[0]*(x-moment[1])**2/2

    return result

def EIt_rel(x):
    result = 0

    for q in qs:
        if q[1]<x:
            result += q[0]*(x-q[1])**3/6
        if q[2]<x:
            result -= q[0]*(x-q[2])**3/6

    for force in forces:
        if force[1]<x:
            result += force[0]*(x-force[1])**2/2

    for moment in Ms:
        if moment[1]<x:
            result += moment[0]*(x-moment[1])

    return result

if vartype == "double":
    EIt0 = (EIy_rel(Xb) - EIy_rel(Xa)) / (Xa - Xb)
elif vartype == "hard":
    EIt0 = -EIt_rel(Xa)

EIy0 = -EIy_rel(Xa)-EIt0*Xa

print("EIθ_0 = {:g} Н*м^2; EIy_0 = {:g} Н*м^3".format(EIt0, EIy0))

def EIt(x):
    return EIt_rel(x) + EIt0

def EIy(x):
    return EIy_rel(x) + x*EIt0 + EIy0


# Определим список координат, при которых будем вычислять все параметры
# для построения эпюр (и других целей): поперечную силу, изгибающий
# момент, угол поворота сечения и прогиб. Будем вычислять их по длине
# балки каждый сантиметр, и через 1 нм после границ участков, чтобы на
# эпюре скачки были вертикальные.

X = list(np.arange(0, Xmax+0.01, 0.01))
X += [q[1]+1e-9 for q in qs]
X += [q[2]+1e-9 for q in qs]
X += [f[1]+1e-9 for f in forces]
X += [m[1]+1e-9 for m in Ms]
X = list(set(X)) # убрать повторяющиеся координаты
X.sort()


# Посчитаем поперечные силы и изгибающие моменты, после чего определим
# опасные сечения и требуемые моменты инерции и сопротивления, чтобы
# выполнялись условия жесткости и прочности.

Q = []
M = []
EIY=[]
for x in X:
    Q.append(q(x))
    M.append(m(x))
    EIY.append(EIy(x))

Qmax = max(Q, key=abs)
Mmax = max(M, key=abs)

I_need = abs(max(EIY, key=abs))/(Ydop*E)
W_need = abs(Mmax)/sigma_dop
print("I_тр = {:g} м^4; W_тр = {:g} м^3".format(I_need, W_need))

# Определим ширину прямоугольного сечения, для которого будут
# выполняться условия прочности по нормальным и касательным напряжениям
# и условие жесткости
b_need_sigma = (6/hb_rel**2*W_need)**(1/3)
b_need_tau = math.sqrt(2/hb_rel * abs(Qmax)/tau_dop)
b_need_I = (I_need*12/hb_rel**3)**(1/4)
b_need = max(b_need_sigma, b_need_tau, b_need_I)
# Округлим до целых см в большую сторону.
b_rect = math.ceil(b_need*100)/100
h_rect = b_rect*hb_rel # высота сечения

print("b_тр = {:g} м; b = {:g} м; h = {:g} м"\
      .format(b_need, b_rect, h_rect))


# Для дальнейших вычислений будет удобнее задействовать
# объектно-ориентированное программирование. Определим классы сечений в
# виде прямоугольника и двутавра (упрощенного). Объекты этих классов
# смогут выдавать свою ширину при заданной y, статический момент
# относительно нейтральной оси части площади сечения, лежащей по одну
# сторону от y, требующиеся для вычисления касательного напряжения по
# формуле Журавского, и изображение, которое потребуется для эпюр
# распределения τ и σ.

class Profile:
    def __init__(self, b, h):
        self.b = b
        self.h = h
        self.Wx = b*h**2/6
        self.Ix = b*h**3/12
        self.A = b*h
        self.name = "Прямоугольник {}x{} м".format(b, h)

    def S(self, y):
        return self.b/2*(self.h**2/4 - y**2)

    def wid(self, y):
        return self.b

    def image(self):
        x = self.b/2
        y = self.h/2
        return ([x, -x, -x, x], [y, y, -y, -y])


class Dvutavr(Profile):
    def __init__(self, l):
        self.h = l[1]
        self.b = l[2]
        self.d = l[3]
        self.t = l[4]
        self.A = l[5]
        self.Ix = l[6]
        self.Wx = l[7]
        self.Sx = l[8]
        self.name = "Двутавр №{}".format(l[0])

    def wid(self, y):
        # за пределами сечения возвращаем ненулевую ширину, чтобы не
        # было неопределенности 0/0
        if abs(y)>=self.h/2-self.t:
            return self.b
        else:
            return self.d

    def S(self, y):
        ycr = self.h/2-self.t
        if abs(y) >= ycr:
            return self.wid(y)/2*(self.h**2/4 - y**2)
        else:
            return self.Sx-(self.Sx-self.S(ycr))*(y/ycr)**2

    def image(self):
        y1 = self.h/2
        y2 = self.h/2-self.t
        x1 = self.b/2
        x2 = self.d/2
        return ([x1, -x1, -x1, -x2, -x2, -x1, -x1, x1, x1, x2, x2, x1],
                [y1, y1, y2, y2, -y2, -y2, -y1, -y1, -y2, -y2, y2, y2])


# По сортаменту найдем двутавр, для которого будут выполняться условия
# прочности и жесткости. Создадим объекты сечений, определим опасные
# сечения (точнее, значения поперечных сил и изгибающих моментов в них)
# и посчитаем нормальные и касательные напряжения в опасных сечениях.
# Попутно проверим условие прочности.

items = []
for item in sortament:
    if item[7] >= W_need and item[6] >= I_need:
        items.append(Dvutavr(item))
        print("Выбран двутавр №{}".format(item[0]))
        break
else:
    print("Не удалось найти подходящий двутавр по сортаменту")

items.append(Profile(b_rect, h_rect))

p_plots = []
for item in items:
    tau = []
    sigma = []
    image = list(item.image())
    ym = item.h/2
    # для эпюр напряжений используем 243 точки.
    yp = np.linspace(-ym, ym, 243)

    for y in yp:
        tau.append(Qmax*item.S(y)/(item.Ix*item.wid(y)))
        sigma.append(-y*Mmax/item.Ix)

    p_plots.append((yp, image, sigma, tau))

    tau_max, sigma_max = max(tau, key=abs), max(sigma, key=abs)
    ok = abs(tau_max)<=tau_dop and abs(sigma_max)<=sigma_dop

    print("Условие прочности {}выполняется для сечения: {}"\
        .format("" if ok else "НЕ ", item.name))

    print("τ_max = {:.4g} Па; σ_max = {:.4g} Па"\
        .format(tau_max, sigma_max))


# Построим эпюры распределения нормальных и касательных напряжений по
# сечению. По горизонтальным осям масштаб одинаковый. Обратите внимание
# на множители, применяемые к значениям по осям напряжений.


fig2 = plt.figure(2, figsize=(10,10))
axs = fig2.subplots(2, 3, sharey="all", sharex="col")

for i, dat in enumerate(p_plots):
    ix, iy = dat[1][0], dat[1][1] # Изображение сечения
    axs[i, 0].fill(ix, iy)
    axs[i, 0].set_aspect(1)
    axs[i, 0].set_ylabel("y / м")
    axs[i, 0].set_xlabel("x / м")

    axs[i, 1].plot(dat[2], dat[0]) # нормальные напряжения
    axs[i, 1].axvline(color="black")
    axs[i, 1].set_xlabel("σ / Па")

    axs[i, 2].axvline(color="black") # касательные
    axs[i, 2].plot(dat[3], dat[0])
    axs[i, 2].set_xlabel("τ / Па")

    for j in range(3):
        axs[i, j].grid()

fig2.savefig("profile.png")


# Сосчитаем углы поворота сечений и прогиб. Проверим условие жесткости


I = items[0].Ix

T = []
Y = []
for x in X:
    Y.append(EIy(x)/(E*I))
    T.append(EIt(x)/(E*I))

Ymax = max(Y, key=abs)

print("Условие жесткости {}выполняется для сечения: {}"\
      .format("НЕ " if abs(Ymax)>Ydop else "", items[0].name))


# Сосчитаем коэффициент экономии


if len(items) > 1:
    econom = items[0].A/items[1].A
    print("Коэффициент экономии: {:.4f}".format(econom))


# Наконец, построим эпюры поперечных сил, изгибающих моментов, углов
# закручивания и прогибов.

fig1 = plt.figure(1, figsize=(10,15))
axs1 = fig1.subplots(4, 1)

for ax, s in zip(axs1, [(Q, "Q / Н"), (M, "М / (Н*м)"),\
                        (T, "θ"), (Y, "y / м")]):
    ax.plot(X, s[0])
    ax.set_ylabel(s[1])
    ax.set_xlabel("x / м")
    ax.axhline(color="black")
    ax.grid()

fig1.savefig("balka.png")


# Построим таблицу значений поперечных сил, изгибающих моментов, углов
# поворота сечения и прогибов в зависимости от координаты:

with open("x_data.txt", "w") as f:
    f.writelines("     x       Q          M       theta       y\r\n")
    for x, q, m, t, y in zip(X, Q, M, T, Y):
        msg = ""
        if q == Qmax:
            msg += "Qmax "
        if m == Mmax:
            msg += "Mmax"
        if y == Ymax:
            msg += "Ymax"
        f.writelines("{:8.4f}{:8.0f}{:12.2f}{:10.5f}{:10.5f} {}\r\n"\
                    .format(x, q, m, t, y, msg))


# Выведем таблицы распределения касательных напряжений по высоте
# сечения. Сначала идет таблица для двутавра (если его удалось найти),
# потом для прямоугольника.

with open("y_data.txt", "w") as f:
    for prof, item in zip(p_plots, items):
        name = item.name
        f.writelines("     y       tau       === {} \r\n".format(name))
        for y, tau in zip(prof[0], prof[3]):
            f.writelines("{:8.4f}{:12.4e}\r\n".format(y, tau))

print("Результаты расчета сохранены в {}".format(os.getcwd()))