1.增加与手动微分公式比较。

2.牛顿法对初值比较敏感。
This commit is contained in:
facat 2020-12-14 10:31:29 +08:00
parent b4db8b612e
commit 2750267b94
4 changed files with 122 additions and 12 deletions

View File

@ -126,6 +126,21 @@ def evaluate_d_delta_l_i_sigma_i(val_delta_l_li, val_sigma_i):
(beta_i, math.atan(h_array[i] / l_array[i])),
]
)
manual_val = exp.manual_diff_delta_li_sigma_i(
h_array[i],
l_array[i],
lambda_m_data,
lambda_i_array[i],
val_sigma_i[i],
sigma_m_data,
alpha_data,
t_data,
t_e_data,
t_m_data,
elastic_data,
)
if math.fabs(val - manual_val) > 1e-5:
raise Exception("d_delta_l_i_sigma_i 自动和手动微分不匹配")
val_list.append(val)
return val_list
@ -201,6 +216,21 @@ def evaluate_d_sigma_i1_d_delta_l_i(val_delta_l_li, val_sigma_i):
for index, li in enumerate(delta_Li):
_val = _val.subs(li, _val_delta_l_li[index])
pass
manual_val = exp.manual_diff_sigma_i1_d_l_i(
h_array[i],
l_array[i],
h_array[i + 1],
l_array[i + 1],
string_g / conductor_n,
area_data,
lambda_i_array[i],
lambda_i_array[i + 1],
val_sigma_i[i],
string_length,
math.fsum(val_delta_l_li[0 : i + 1]),
)
if math.fabs(manual_val - _val) > 1e-5:
raise Exception("d_sigma_i1_delta_L_i 自动和手动微分不匹配")
col.append(_val)
row.append(col)
return sympy.Matrix(row)
@ -237,11 +267,11 @@ def evaluate_d_sigma_i1_d_delta_sigma_i(val_delta_li):
def solve():
# 初始化
val_delta_li = [0.1 for i in range(span_count)]
# val_delta_li = [0.15864687475316822, -0.1935189734784845, 0.03478489898855073]
val_delta_li = [0.1 for _ in range(span_count)]
# val_delta_li = [0.15931589580330385, -0.19294219226439696, 0.035236744603670586, -0.03124795962518869, 0.034422847061933395, 0.1639551737321582, -0.18876605469917845, 0.03862307099222713, -0.01858088432390902]
val_sigma_i = [sigma_m_data for _ in range(span_count)]
# val_sigma_i = [175.38451579479482, 176.01015153076614, 175.88355419459572]
val_sigma_i = [175.44277576372207, 176.07105062365383, 175.94897602628689, 175.95574563415994, 175.83116207025824, 175.84778616876434, 176.5132398968081, 176.4405408454816, 176.51871673188586]
# val_sigma_i = [sigma_m_data for _ in range(span_count)]
loop = 0
while True:
@ -333,7 +363,7 @@ def solve():
x_list = list(x)[0]
abs_min = [math.fabs(_x) for _x in x_list]
abs_min.sort()
if abs_min[-1] < 1e-5:
if abs_min[-1] < 1e-6:
break
print("最大偏差{max_dx}".format(max_dx=abs_min[-1]))
# 更新变量
@ -346,7 +376,7 @@ def solve():
print(loop)
print(val_delta_li)
print(val_sigma_i)
verify(val_delta_li,val_sigma_i)
verify(val_delta_li, val_sigma_i)
def verify(val_delta_li, val_sigma_i):

32
data.py
View File

@ -1,4 +1,4 @@
#loop_end = 10000000 # 最大循环次数
# loop_end = 10000000 # 最大循环次数
loop_end = 1000000 # 最大循环次数
# 架线时的状态
# 取外过无风
@ -10,13 +10,35 @@ alpha = 0.0000155 # 导线膨胀系数 1/°C
elastic = 95900 # 弹性系数 N/mm2
area = 154.48 # 导线面积 mm2
lambda_m = 14.8129 / area # 导线比载 N/(m.mm)
lambda_i_array = [lambda_m*0.9,lambda_m*1.3,lambda_m,lambda_m,lambda_m]
lambda_i_array = [
lambda_m * 0.9,
lambda_m * 1.3,
lambda_m,
lambda_m,
lambda_m,
lambda_m * 0.9,
lambda_m * 1.3,
lambda_m,
lambda_m,
lambda_m,
]
# 取400m代表档距下
sigma_m = 28517 / area # 架线时初伸长未释放前的最低点水平应力。单位N/mm2
span_count = 3 # 几个档距
span_count = 9 # 几个档距
# n个档距,n-1个直线塔
h_array = [0, 0, 0, 0, 0]
l_array = [400, 300, 300, 500, 300]
h_array = [
0,
0,
0,
0,
0,
0,
0,
0,
0,
0,
]
l_array = [400, 300, 300, 500, 300, 400, 300, 300, 500, 300]
t = 15
epsilon = 1e-4 # 收敛判据
conductor_n = 6 # 导线分裂数

56
exp.py
View File

@ -27,7 +27,7 @@ def delta_li():
t_m,
sigma_m,
sigma_i,
beta_i
beta_i,
) = sympy.symbols(
"""
delta_l_i,
@ -106,6 +106,7 @@ def fun_sigma_i1(delta_Li):
for f in delta_Li:
_t += f
return _t
_sigma_i1 = sigma_i1 - (
(
G_i / 2 / A # G_i传入时已考虑导线分裂数
@ -116,3 +117,56 @@ def fun_sigma_i1(delta_Li):
+ sigma_i / b_i() * sympy.sqrt(stringlen_i ** 2 - b_i() ** 2)
) / (sympy.sqrt(stringlen_i ** 2 - b_i() ** 2) / b_i() + h_i1 / l_i1)
return _sigma_i1
def manual_diff_delta_li_sigma_i(
h_i, l_i, lambda_m, lambda_i, sigma_i, sigma_m, alpha, t_i, t_e, t_m, E
):
beta_i = math.atan(h_i / l_i)
A = (
(l_i * math.cos(beta_i)) ** 2
/ 24
* ((lambda_m / sigma_m) ** 2 - (lambda_i / sigma_i) ** 2)
+ (sigma_i - sigma_m) / E / math.cos(beta_i)
+ alpha * (t_i + t_e - t_m)
)
B = 1 + (lambda_i * l_i / sigma_i) ** 2 / 8
dA_dsigma_i = ((l_i * math.cos(beta_i)) ** 2) / 24 * 2 * lambda_i ** 2 * (
sigma_i ** (-3)
) + 1 / E / math.cos(beta_i)
dB_dsigma_i = -2 * (lambda_i * l_i) ** 2 / 8 * (sigma_i ** (-3))
_t = -l_i / (math.cos(beta_i) ** 2) * (dA_dsigma_i * B - A * dB_dsigma_i) / (B ** 2)
return _t
def manual_diff_sigma_i1_d_l_i(
h_i, l_i, h_i1, l_i1, Gi, A, lambda_i, lambda_i1, sigma_i, stringlen, b_i
):
beta_i = math.atan(h_i / l_i)
beta_i1 = math.atan(h_i1 / l_i1)
A = (
Gi / 2 / A
+ lambda_i * l_i / 2 / math.cos(beta_i)
+ lambda_i1 * l_i1 / 2 / math.cos(beta_i1)
+ sigma_i * h_i / l_i
+ sigma_i / b_i * ((stringlen ** 2 - b_i ** 2) ** 0.5)
)
B = ((stringlen ** 2 - b_i ** 2) ** 0.5)/b_i + h_i1 / l_i1
dA_d_delta_L1 = (
sigma_i
* (
-((stringlen ** 2 - b_i ** 2) ** -0.5) * (b_i ** 2)
- (stringlen ** 2 - b_i ** 2) ** 0.5
)
/ (b_i ** 2)
)
dB_d_delta_L1 = (
1
/ (b_i ** 2)
* (
-((stringlen ** 2 - b_i ** 2) ** -0.5) * (b_i ** 2)
- (stringlen ** 2 - b_i ** 2) ** 0.5
)
)
_t = -(dA_d_delta_L1 * B - A * dB_d_delta_L1) / (B ** 2)
return _t

View File

@ -257,6 +257,10 @@ def verify(
print("!!!应力等式不满足")
return
print("等式满足。")
print('sigma')
print(sigma_array)
print('delta_li')
print(delta_l_i_array)
if __name__ == "__main__":