## Dymola中计算Fan2ndOrder
## 选取额定工况为1600转
import copy
import pandas as pd
import numpy as np
from scipy.optimize import curve_fit as cf
from matplotlib import pyplot as plt
## 将不同转速下的流量-功率曲线和流量-静压曲线拟合成三次型函数
## 系统的静压曲线使用指数型关系式代替
ALPHA = 0.2
BETA = 0.5
class FanCurve:
def __init__(self):
self.curve_power = {
1280: [np.float64(-1.0939824906995037e-09), np.float64(2.333022185311005e-05), np.float64(-0.06796629659302983), np.float64(2055.3002100486133)],
1440: [np.float64(-1.136696332147217e-09), np.float64(2.7660598448342913e-05), np.float64(-0.09613228325870021), np.float64(2922.59538328202)],
1500: [np.float64(-1.7603856529124879e-09), np.float64(4.4241229797522305e-05), np.float64(-0.14983705834350444), np.float64(5160.7745241685625)],
1600: [np.float64(-1.1701720186561717e-09), np.float64(3.0605089172086764e-05), np.float64(-0.0989935007927685), np.float64(3869.5082357198503)]
}
self.curve_dp = {
1280: [np.float64(-3.267975749307142e-10), np.float64(5.2912277771354885e-06), np.float64(-0.03403187464621037), np.float64(741.9041721827268)],
1440: [np.float64(-2.422194407553435e-10), np.float64(3.7698924924279502e-06), np.float64(-0.02559742162953265), np.float64(931.5430549263344)],
1500: [np.float64(-2.3878745827961e-10), np.float64(3.921270379695461e-06), np.float64(-0.02640953385995719), np.float64(1002.1261618905066)],
1600: [np.float64(-2.2348965671136672e-10), np.float64(3.886724444389961e-06), np.float64(-0.027754898177200903), np.float64(1147.392301537067)]
}
self.speed = [1280, 1440, 1500, 1600] # 有顺序的 不要改
self.alpha = ALPHA
self.beta = BETA
self.q_min = 0
self.q_max = 24000
def pressure_q(self, q):
return self.alpha * q ** self.beta
@staticmethod
def curve_function(q, a, b, c, d):
return a * q **3 + b * q **2 + c * q + d
@staticmethod
def mk_function(a, b, c, d):
# 此函数适用于不插值的计算
def curve_function(q):
return a * q ** 3 + b * q ** 2 + c * q + d
return curve_function
@staticmethod
def mk_inter_function(percent, a, b, c, d, e, f, g, h):
## 如果使用插值,就用这个方法构造函数
## abcd是左侧值的参数,efgh是右侧的
def inter_curve_function(q):
l_value = a * q ** 3 + b * q ** 2 + c * q + d
r_value = e * q ** 3 + f * q ** 2 + g * q + h
return (1 - percent) * l_value + percent * r_value
return inter_curve_function
@staticmethod
def mk_outside_function(coefficient, a, b, c, d):
def curve_function(q):
return (a * q ** 3 + b * q ** 2 + c * q + d) * coefficient
return curve_function
def get_functions(self, n):
"""
使用转速和系统的阻力特性计算体积流量
转速在性能曲线范围内的采用插值处理
转速在性能曲线外的采用风机相似定律
:param n:
:return:
"""
if 1280 <= n <= 1600: ## 插值得出风机的静压,电功率
if n in self.speed: ## 给定的转速在风机性能曲线内,不需要插值
power_params = self.curve_power[n]
power_function = self.mk_function(*power_params)
dp_params = self.curve_dp[n]
dp_function = self.mk_function(*dp_params)
else: ## 给定的转速不在风机性能曲线上,则需要进行插值
new_speed = copy.deepcopy(self.speed)
new_speed.append(n)
new_speed.sort()
idx = 0
for i in range(len(new_speed)):
if n == new_speed[i]:
idx = i
break
l_speed = new_speed[idx - 1]
r_speed = new_speed[idx + 1]
l_power_params = self.curve_power[l_speed]
r_power_params = self.curve_power[r_speed]
l_dp_params = self.curve_dp[l_speed]
r_dp_params = self.curve_dp[r_speed]
percent = (n - l_speed) / (r_speed - l_speed)
power_function = self.mk_inter_function(percent, *l_power_params, *r_power_params)
dp_function = self.mk_inter_function(percent, *l_dp_params, *r_dp_params)
else: ## 使用风机相似定律求出风机电功率、静压
if n < 1280:
volume_coefficient = n / 1280
else:
volume_coefficient = n / 1600
pressure_coefficient = volume_coefficient ** 2
power_coefficient = volume_coefficient ** 3
dp_params = self.curve_dp[1280]
power_params = self.curve_power[1280]
dp_function = self.mk_outside_function(pressure_coefficient, *dp_params)
power_function = self.mk_outside_function(power_coefficient, *power_params)
return {"dp_function": dp_function, "power_function": power_function}
def solve(self, f, g, q_init=10000, convergence = 1e-06, max_iter = 10000):
## 暂时不对无解情况进行规避
count = 0
def difference(q):
return f(q) - g(q)
l, r = self.q_min, self.q_max
if difference(self.q_min) == 0:
return self.q_min
if difference(self.q_max) == 0:
return self.q_max
while True:
if difference(l) * difference(r) > 0:
raise ValueError("使用二分法时,似乎存在多解,无法求得确切值")
if 0 < difference(q_init) <= convergence:
break
if difference(q_init) * difference(l) > 0: ## 跟左边的值同号
l = q_init
q_init = (l + r) / 2
else:
r = q_init
q_init = (l + r) / 2
count += 1
if count > max_iter:
raise TimeoutError("迭代次数已经超过了最大迭代次数")
return q_init
def run(self, n):
functions = self.get_functions(n)
q = self.solve(functions["dp_function"], self.pressure_q)
return q
if __name__ == '__main__':
# ## 读取数据
# dataframe_dp = pd.read_excel(r"E:\projects\air_conditioner\fan\v_dp_1600.xlsx")
#
# # dataframe_power = pd.read_excel(r"E:\projects\air_conditioner\fan\power_1600.xlsx")
#
# ## 拟合
#
# def f_power(x, a, b, c, d):
# return a * x **3 + b * x **2 + c * x + d
#
#
# for file_name in ["power_1280", "power_1440", "power_1500", "power_1600"]:
# file_path = fr"E:\projects\air_conditioner\fan\{file_name}.xlsx"
# dataframe_power = pd.read_excel(file_path)
#
# x_data = np.array(dataframe_power['x'])
# y_data = np.array(dataframe_power['y'])
#
# a, b, c, d = cf(f_power, x_data, y_data)[0]
# xx = [a, b, c, d]
# print(file_name, f"x={xx}")
n_set = 1400
solution = FanCurve()
q = solution.run(n_set)
p = solution.pressure_q(q)
functions = solution.get_functions(n_set)
print("性能曲线压力:", functions["dp_function"](q))
print(q)
print("阻力压力:", p)
请帮我看看这段代码中,是如何取不同转速下性能曲线函数中的几个参数值的