t.y. для всех комментариев.
Будучи новичком, еще раз прошу прощения за размещение такого материала.
Все хорошие идеи приходят утром.
Проблема заключалась в том, что код выполнялся даже без цикла FOR. Дело в том, что каждый раз, когда код достигает последней вызванной итерации.
Решение: я создал глобальную переменную для вновь рассчитанных значений после i-го расчета и добавил ее в новый глобальный пустой список. Затем во второй итерации он извлекал значения из глобальной переменной, то есть ему не нужно снова вычислять все вещи.
P.S. Если some1 хочет вычислить немного инженерных разработок для нелинейного железобетона, просто измените заданные глобальные переменные.
import matplotlib
import matplotlib.pyplot as plt
import numpy as np
#Needed variables and constants
n = 2 #Degree of nonlinearity(natural number ex.2,3,4 etc.)
a1 = 0.02 #Concrete cover of tensile reinf.[m]
a2 = 0.02 #Concrete cover of compressive reinf.[m]
As1 = 20 * 10**(-4) #Area of tensile reinf.[m**2]
As2 = 20 * 10**(-4) #Area of compressive reinf.[m**2]
e0 = 0.04 #Initial excetricity in [m]
b = 0.3 #Width of the columns cross-section.[m]
h = 0.3 #Height of the columns cross-section.[m]
Ned = 2600 * 10**3 #Axial Load.[N]
fcd = 8 * 10**6 #Concrete strength.[Pa]
Ecm = 27.6 * 10**9 #Concrete modulus of elasticity.[Pa]
Es = 200 * 10**9 #Steel modulus of elasticity. [Pa]
d = h - a1 #Effective height of tensile reinf. [m]
e = (h/2) + e0 #Excentricity of Ned[m]
Eps2 = 0.002 #Deformaton when stresses in concrete reaches fcd.
Eps35 = 0.0035 #Limit deformation for concrete.
d_x0 = 5 * 10**(-3) #Increment for x. Needed when calculating derivatives of Mb and Nb.
d_r0 = 5 * 10**(-5) #Increment fr r. Needed when calculating derivatives of Mb and Nb.
fyd = 600 * 10**6 #Yield stress for steel.[Pa]
d_x_n = []
d_r_n = []
#Calculation of nonlinear reinforced concrete.
#This code calculates 2 unknonws with system made out of 2 equations.
#Unknownws are physicaly connected in nonlinear dependecies.
#Solving this nonlinear system of equations is done by making linear iterative calculation, based on Taylor series.
#Program stops, when 2 adjacent iterations give same values.(at least 0.9999 accuracy)
def x(i):
if i == 1:
return h
x_n = x(i-1) + d_x(i-1)
return x_n
def r(i):
if i == 1:
return (Ned * e) / ((Ecm * b * (h ** 3)) / 12)
r_n = r(i-1) + d_r(i-1)
return r_n
def d_r(i):
return d_r_n[i-1]
def d_x(i):
return d_x_n[i-1]
def A(i):
return np.array([[f1_x(i), f1_r(i)], [f2_x(i), f2_r(i)]])
def B(i):
return np.array([-f1(i), -f2(i)])
def C(i):
return np.linalg.solve(A(i), B(i))
def Nb(i):
return (b * fcd / r(i)) * (r(i) * h - (Eps2 / (n + 1)) * (1 - (r(i) * (x(i) - h)) / Eps2) ** (n + 1))
def Mb(i):
return (b * fcd / (r(i) ** 2)) * (-(((Eps2 - r(i) * (x(i) - h)) ** (n + 1)) * (r(i) * (x(i) - h) * (n + 1) + Eps2) / ((n + 2) * (n + 1) * (Eps2 ** n))) + (r(i) ** 2) * h * (x(i) - (h / 2)))
def Ns1(i):
return Es * r(i) * (x(i) - d) * As1
def Ns2(i):
return Es * r(i) * (x(i) - a2) * As2
def Ms1(i):
return Es * r(i) * ((x(i) - d) ** 2) * As1
def Ms2(i):
return Es * r(i) * ((x(i) - a2) ** 2) * As2
def f1(i):
return Nb(i) + Ns1(i) + Ns2(i) - Ned
def f2(i):
return Mb(i) + Ms1(i) + Ms2(i) - (Ned * (x(i) - (h / 2) + e0))
# Derivatives of all the required elements
def Ns1_x(i):
return Es * r(i) * As1
def Ns1_r(i):
return Es * As1 * (x(i) - d)
def Ns2_x(i):
return Es * r(i) * As2
def Ns2_r(i):
return Es * As2 * (x(i) - a2)
def Ms1_x(i):
return Es * r(i) * As1 * 2 * (x(i) - d)
def Ms1_r(i):
return Es * As1 * (x(i) - d) ** 2
def Ms2_x(i):
return Es * r(i) * As2 * 2 * (x(i) - a2)
def Ms2_r(i):
return Es * As2 * (x(i) - a2) ** 2
def Nb_x(i):
return (((b * fcd / r(i)) * (r(i) * h - (Eps2 / (n + 1)) * (1 - (r(i) * ((x(i) + d_x0) - h)) / Eps2) ** (n + 1))) -((b * fcd / r(i)) * (r(i) * h - (Eps2 / (n + 1)) * (1 - (r(i) * ((x(i) - d_x0) - h)) / Eps2) ** (n + 1)))) / (2 * d_x0)
def Nb_r(i):
return (((b * fcd / (r(i) + d_r0)) * ((r(i) + d_r0) * h - (Eps2 / (n + 1)) * (1 - ((r(i) + d_r0) * (x(i) - h)) / Eps2) ** (n + 1))) -((b * fcd / (r(i) - d_r0)) * ((r(i) - d_r0) * h - (Eps2 / (n + 1)) * (1 - ((r(i) - d_r0) * (x(i) - h)) / Eps2) ** (n + 1)))) / (2 * d_r0)
def Mb_x(i):
return (((b * fcd / (r(i) ** 2)) * (-(((Eps2 - r(i) * ((x(i)+d_x0) - h)) ** (n + 1)) * (r(i) * ((x(i)+d_x0) - h) * (n + 1) + Eps2) / ((n + 2) * (n + 1) * (Eps2 ** n))) + (r(i) ** 2) * h * ((x(i)+d_x0) - (h / 2))))-((b * fcd / (r(i) ** 2)) * (-(((Eps2 - r(i) * ((x(i)-d_x0) - h)) ** (n + 1)) * (r(i) * ((x(i)-d_x0) - h) * (n + 1) + Eps2) / ((n + 2) * (n + 1) * (Eps2 ** n))) + (r(i) ** 2) * h * ((x(i)-d_x0) - (h / 2))))) / (2 * d_x0 )
def Mb_r(i):
return (((b * fcd / ((r(i)+d_r0) ** 2)) * (-(((Eps2 - (r(i)+d_r0) * (x(i) - h)) ** (n + 1)) * ((r(i)+d_r0) * (x(i) - h) * (n + 1) + Eps2) / ((n + 2) * (n + 1) * (Eps2 ** n))) + ((r(i)+d_r0) ** 2) * h * (x(i) - (h / 2))))-((b * fcd / ((r(i)+d_r0) ** 2)) * (-(((Eps2 - (r(i)+d_r0) * (x(i) - h)) ** (n + 1)) * ((r(i)+d_r0) * (x(i) - h) * (n + 1) + Eps2) / ((n + 2) * (n + 1) * (Eps2 ** n))) + ((r(i)+d_r0) ** 2) * h * (x(i) - (h / 2))))) / (2 * d_r0)
def f1_x(i):
return Nb_x(i) + Ns1_x(i) + Ns2_x(i)
def f1_r(i):
return Nb_r(i) + Ns1_r(i) + Ns2_r(i)
def f2_x(i):
return Mb_x(i) + Ms1_x(i) + Ms2_x(i) - Ned
def f2_r(i):
return Mb_r(i) + Ms1_r(i) + Ms2_r(i)
# Results of iterations
def Eps_C2(i):
return x(i) * r(i)
def Eps_C1(i):
return (x(i) - h) * r(i)
def Sig_S1(i):
return Es * r(i) * (x(i) - d)
def Sig_S2(i):
return Es * r(i) * (x(i) - a2)
for i in range(1,20):
d_x_n.append((C(i)[0]))
d_r_n.append((C(i)[1]))
print('Iteration', i)
print(x(i))
if round(x(i)/x(i+1), 4) == 1 and round(r(i)/r(i+1), 4) == 1:
print('Ratio between two iterations ', round(x(i)/x(i+1), 4) )
print('Convergence reached in ', i, 'Iterations')
break
else:
print('Convergence was not reached in ', i, 'Iterations, calculation aborts')
print('Stress in tensile reinforcement,', Sig_S1(i))
print('Stress in compressive reinforcement,', Sig_S2(i))
print('Deformations in C1', Eps_C1(i))
print('Deformations in C2', Eps_C2(i))