Начало
Всё началось с того, что мне подкинули интересный проект с разработкой тороидальных винтов.
Всё началось в solidworks с того, что был смоделирован стандартный двухлопастной винт с изменяемыми углами атаки через переменные.
Далее начался подбор параметров для запуска исследования:
Для начала использовали периодичную расчетную область, так как где-то вычитали то, что для продувки винтов необходимо использовать именно её, но через неделю активного недоумевания, почему при изменениях расчетной области менялась тяга, мы поняли, как работает периодическая расчетная область
То есть потоки расчетной области копируются в соседние ячейки с такой же периодичностью, как и область расчета и да всё это делается во вращении локальной области Averaging — это где вращается модель в сетке, в остальном всё такие же настройки, как и во всех гайдах. Также было непонимание почему 5 дюймовый винт на 25 000 об/мин создаёт тягу в 0,001Н.
Затем просто в научных интересах во вращении локальной области поставил Sliding — в нем сетка вращается вокруг модели, то есть берется цилиндр сетки и начинает вращаться в сетке. И, о чудо, были получены правдоподобные результаты.
На следующем этапе стало интересно, что будет, если покрутить разрешение сетки, ведь от него время «Продувки » варьируется от пары секунд до недели не на самом слабом ПК с intel 9 10 940 и quadro a5000, для этого был отсканирован винт сканером shining 3d ue pro. Так, к слову, у него погрешность сканирования 0,01 мм и произведен реверс инжиниринг для оптимизации геометрии. В качестве знака окончания продувки смотрели график тяги от итераций и останавливали «продувку», когда тяга выходила на «полку».
Тяга винта на испытательном стенде была 6,75Н и тут появились вопросы к solidworks, к этому прибавляется, что было сказано всё это делать на отечественном софте, из‑за этого пересели на компас, но в нем встроенный модуль cfd расчета не позволяет «продувать» подвижные детали. Было принято решение использовать open source проект fluid x3d, который в отличие от solidworks производит расчет на видеокарте вместо процессора, что сокращает время расчета с недели до нескольких минут! Итак, после первых тестов fluid были написаны следующие настройки для продувки винта и экспорта значений продувки в txt файл.
void main_setup() {
// ################################################################## define simulation box size, viscosity and volume force ###################################################################
const uint memory = 3500u; // available VRAM of GPU(s) in MB
const float lbm_u = 0.12f;
const float box_scale = 6.0f;
const float si_u = 0.17f;
const float si_nu = 0.17f, si_rho = 1.0f;
const float si_width = 0.3f, si_height = 0.3f, si_length = 0.3f;
const float si_A = si_width * si_height + 2.0f * 0.05f * 0.03f;
const float si_T = 0.25f;
const float si_Lx = units.x(box_scale * si_width);
const float si_Ly = units.x(box_scale * si_length);
const float si_Lz = units.x(0.5f * (box_scale - 1.0f) * si_width + si_height);
const uint3 lbm_N = resolution(float3(3.0f, 3.0f, 1.0f), memory); // input: simulation box aspect ratio and VRAM occupation in MB, output: grid resolution
units.set_m_kg_s((float)lbm_N.y, lbm_u, 1.0f, box_scale * si_length, si_u, si_rho);
print_info("Re = " + to_string(to_uint(units.si_Re(si_width, si_u, si_nu))));
const float lbm_nu = units.nu(si_nu);
const float si_omega = 1000u;
const float lbm_omega = units.omega(si_omega);
const float lbm_length = units.x(si_length);
LBM lbm(lbm_N, lbm_nu);
const uint lbm_T = 5000u;
const uint lbm_dt = 10u;
const float summa_cd = 0;
float cd[10];
float avergate_cd[10];
int a = 0;
float avergate = 0;
int i;
const float radius = 0.25f * (float)lbm_N.x;
const float3 center = float3(lbm.center().x, lbm.center().y, 0.36f * radius);
// ###################################################################################### define geometry ######################################################################################
Mesh* mesh = read_stl(get_exe_path() + "../stl/33.stl", lbm.size(), lbm.center(), float3x3(float3(0, 0, 1), radians(90.0f)), lbm_length);
mesh->translate(float3(0.0f, units.x(0.5f * (0.5f * box_scale * si_length - si_width)) - mesh->pmin.y, 1.0f - mesh->pmin.z));
lbm.voxelize_mesh_on_device(mesh, TYPE_S | TYPE_X); // https://github.com/nathanrooy/ahmed-bluff-body-cfd/blob/master/geometry/ahmed_25deg_m.stl converted to binary
const uint Nx = lbm.get_Nx(), Ny = lbm.get_Ny(), Nz = lbm.get_Nz(); parallel_for(lbm.get_N(), [&](ulong n) { uint x = 0u, y = 0u, z = 0u; lbm.coordinates(n, x, y, z);
//if(z==0u) lbm.flags[n] = TYPE_S;
//if(lbm.flags[n]!=TYPE_S) lbm.u.y[n] = lbm_u;
//if(x==0u||x==Nx-1u||y==0u||y==Ny-1u||z==Nz-1u) lbm.flags[n] = TYPE_E;
}); // ####################################################################### run simulation, export images and data ##########################################################################
//lbm.graphics.visualization_modes = VIS_FLAG_SURFACE | VIS_Q_CRITERION;
//lbm.graphics.set_camera_centered(20.0f, 30.0f, 0.0f, 1.648722f);
lbm.run(0u); // initialize simulation
#if defined(FP16S)
const string path = get_exe_path() + "FP16S/" + to_string(memory) + "MB/";
#elif defined(FP16C)
const string path = get_exe_path() + "FP16C/" + to_string(memory) + "MB/";
#else // FP32
const string path = get_exe_path() + "FP32/" + to_string(memory) + "MB/";
#endif // FP32
//lbm.write_status(path);
//write_file(path+"Cd.dat", "# ttCdn");
//std::ofstream app; // поток для записи
//app.open("force.txt"); // открываем файл для записи
while (true) { // main simulation loop
while (true) { // main simulation loop
//if (lbm.graphics.next_frame(units.t(si_T), 5.0f)) {
Clock clock;
lbm.run(lbm_dt);
lbm.calculate_force_on_boundaries();
lbm.F.read_from_device();
lbm.voxelize_mesh_on_device(mesh, TYPE_S, center, float3(0.0f), float3(0.0f, 0.0f, lbm_omega));
const float3 lbm_force = lbm.calculate_force_on_object(TYPE_S | TYPE_X);
const float Cd = units.si_F(lbm_force.y); // expect Cd to be too large by a factor 1.3-2.0x; need wall model
const float moment = units.si_F(lbm_force.x);
cd[a] = Cd;
a = a +1;
if (a == 10) {
int a = 0;
}
for (i = 0; i < 10; i++) {
avergate = (avergate + cd[i]) / 10;
}
avergate_cd[a] = avergate;
println("r" + to_string(Cd, 3u) + " " + to_string(moment, 3u) + " ");
// write_line(path+"Cd.dat", to_string(lbm.get_t())+"t"+to_string(Cd, 3u)+"n");
mesh->rotate(float3x3(float3(0, 0, 1), lbm_omega * (float)lbm_dt));// rotate mesh
lbm.update_fields();
lbm.run(lbm_dt);
std::ofstream ate;
ate.open("force.txt", std::ios_base::app);
ate << (to_string(Cd,3u) + " " + to_string(moment,3u)) << std::endl;
ate.close();
//if (a > 3) {
// if (avergate_cd[a] - avergate_cd[a - 3] < 0.1 && avergate_cd[a] - avergate_cd[a - 3] > -0.1) {
// break;
// }
//}
#if defined(GRAPHICS) && !defined(INTERACTIVE_GRAPHICS)
// lbm.graphics.write_frame(path+"images/");
#endif // GRAPHICS && !INTERACTIVE_GRAPHICS
}
lbm.run(1u);
//app.close();
}
//lbm.write_status(path);
} /**/
//}https://github.com/ProjectPhysX/FluidX3D - ссылка на fluid кому интересно будет поковыряйтесь советую
Теперь нам надо импортировать данные продувки в python ибо не на чем другом я программировать не умею
import os
import subprocess
def run():
command = ["путь к fluid.exe","/c","runas",'/user:administrator',"regedit"]
p = subprocess.Popen(command, stdin = subprocess.PIPE)
#p.stdin.write("1961")
p.communicate()
def get_result():
data = open("путь к force.txt","r")
a = data.read()
a = a[:-1]
l = a.split("n")
force = 0
moment = 0
#print(l)
for i in range(len(l)):
l[i] = l[i].split(" ")
force += abs(float(l[i][0]))
moment += abs(float(l[i][1]))
force = round(force/len(l),2)
moment = round(moment/len(l),2)
a = [force,moment]
return a
def remove_dat():
os.remove("путь к force.txt")Теперь к компасу. Там есть два варианта. Начнем с того, что попроще.
Рисуем стандартный двухлопастной винт по сечениям, причем в эскизах делаем углы, которые выносим в глобальные переменные.
Теперь пора заняться api компаса
# -*- coding: utf-8 -*-
#|1
import pythoncom
from win32com.client import Dispatch, gencache
import KompasAPI7
import LDefin3D
import MiscellaneousHelpers as MH
import os
def connect():
global iPart,iDocument3D,kompas_document
# Подключим константы API Компас
kompas6_constants = gencache.EnsureModule("{75C9F5D0-B5B8-4526-8681-9903C567D2ED}", 0, 1, 0).constants
kompas6_constants_3d = gencache.EnsureModule("{2CAF168C-7961-4B90-9DA2-701419BEEFE3}", 0, 1, 0).constants
# Подключим описание интерфейсов API5
kompas6_api5_module = gencache.EnsureModule("{0422828C-F174-495E-AC5D-D31014DBBE87}", 0, 1, 0)
kompas_object = kompas6_api5_module.KompasObject(Dispatch("Kompas.Application.5")._oleobj_.QueryInterface(kompas6_api5_module.KompasObject.CLSID, pythoncom.IID_IDispatch))
MH.iKompasObject = kompas_object
# Подключим описание интерфейсов API7
kompas_api7_module = gencache.EnsureModule("{69AC2981-37C0-4379-84FD-5DD2F3C0A520}", 0, 1, 0)
application = kompas_api7_module.IApplication(Dispatch("Kompas.Application.7")._oleobj_.QueryInterface(kompas_api7_module.IApplication.CLSID, pythoncom.IID_IDispatch))
MH.iApplication = application
Documents = application.Documents
# Получим активный документ
kompas_document = application.ActiveDocument
#kompas_document_3d = kompas_api7_module.IKompasDocument3D(kompas_document)
iDocument3D = kompas_object.ActiveDocument3D()
iPart = iDocument3D.GetPart(LDefin3D.pTop_Part)
return iPart,iDocument3D
def rebuild(iPart,iDocument3D):
iPart.RebuildModel()
iDocument3D.RebuildDocument()
def export_stl():
kompas_document.SaveAs("path\FluidX3D-master\stl\33.STL")
def edit_param(iPart,angle1,angle2):
global VariableCollection
VariableCollection = iPart.VariableCollection()
VariableCollection.GetByName("angle1",True,True).value = angle1
VariableCollection.GetByName("angle2",True,True).value = angle2
def remove_stl():
"""
Remove existing STL file
"""
# Delete file
os.remove("path\FluidX3D-master\stl\33.STL")Тут ничего сложного, просто подключаемся к компасу, меняем две переменные и сохраняем/удаляем stl.
Теперь для полного счастья нужен генетический алгоритм:
import random
import api_kompas
import api_fluid
def create_obj(size,obj): #создаем рандомный объект
obj = []
for i in range(size):
obj.append(random.randint(0,45))
return obj
def create_population(pop_size,obj_size): #создем рандомное поколени
population = []
obj = []
for i in range(pop_size):
population.append(create_obj(obj_size,obj))
return population
def calc_obj(iPart,iDocument,obj): #считаем тягу и момент для объекта
api_kompas.edit_param(iPart,obj[0],obj[1])
api_kompas.rebuild(iPart,iDocument)
api_kompas.export_stl()
api_fluid.run()
obj.append(api_fluid.get_result()[0])
obj.append(api_fluid.get_result()[1])
api_fluid.remove_dat()
api_kompas.remove_stl()
def calc_pop(iPart,iDocument,population): #считаем тягу и момент для популяции
for i in range(len(population)):
calc_obj(iPart,iDocument,population[i])
def calc_k_obj(obj): #считаем коофицент выживаемости для объекта
obj.append(obj[len(obj)-2]/obj[len(obj)-1])
def calc_k_pop(population): #считаем коофицент выживаемости для поколения
for i in range(len(population)):
calc_k_obj(population[i])
def sort_population(population):
for i in range(len(population)-1):
for g in range(len(population)-i-1):
if population[g][len(population[i])-1] < population[g+1][len(population[i])-1]:
population[g][len(population[i])-1],population[g+1][len(population[i])-1] = population[g+1][len(population[i])-1],population[g][len(population[i])-1]
def get_max_k(population):
sort_population(population)
m = population[0][len(population[0])-1]
return m
def get_weight(population): #получаем коофиценты выживаемости поколения в отдельный массив
weight = []
sort_population(population)
for i in range(len(population)):
weight.append(population[i][len(population[0])-1]/get_max_k(population))
#print(weight)
return weight
def get_parents_for_new_obj(population): #определяем родителей для объекта следующего поколения
parents = []
sort_population(population)
parents.append(random.choices(population,weights= get_weight(population))[0])
parent = random.choices(population,weights= get_weight(population))[0]
while parents[0] == parent:
parent = random.choices(population,weights= get_weight(population))[0]
parents.append(parent)
return parents
def get_parent_for_new_population(population): #определяем родителей для всего поколения
parents_population = []
for i in range(len(population)):
parents_population.append(get_parents_for_new_obj(population))
return parents_population
def select(parents): #размножаем
new_obj = [0,0]
new_obj[0] = parents[0][0]
new_obj[1] = parents[1][1]
print(new_obj)
return new_obj
def select_population(parents): #размножаем поколение
new_population = []
for i in range(len(parents)):
new_population.append(select(parents[i]))
return new_population
def mutate(population): #мутируем
for i in range(len(population)//2):
for g in range(len(population[0])):
gen = random.randint(0,len(population[0])-1)
population[i][gen] = random.randint(0,45)
Все функции подписаны так, что должно быть всё понятно. Не спорю, что он не идеален, - всё делалось за ночь до выступления :-)
Итак, теперь сделаем простенький main для генетического алгоритма:
import api_kompas
import api_fluid
import pythoncom
from win32com.client import Dispatch, gencache
import KompasAPI7
import LDefin3D
import MiscellaneousHelpers as MH
import exsel
import time
import gen_alg
# Подключим константы API Компас
kompas6_constants = gencache.EnsureModule("{75C9F5D0-B5B8-4526-8681-9903C567D2ED}", 0, 1, 0).constants
kompas6_constants_3d = gencache.EnsureModule("{2CAF168C-7961-4B90-9DA2-701419BEEFE3}", 0, 1, 0).constants
# Подключим описание интерфейсов API5
kompas6_api5_module = gencache.EnsureModule("{0422828C-F174-495E-AC5D-D31014DBBE87}", 0, 1, 0)
kompas_object = kompas6_api5_module.KompasObject(Dispatch("Kompas.Application.5")._oleobj_.QueryInterface(kompas6_api5_module.KompasObject.CLSID, pythoncom.IID_IDispatch))
MH.iKompasObject = kompas_object
# Подключим описание интерфейсов API7
kompas_api7_module = gencache.EnsureModule("{69AC2981-37C0-4379-84FD-5DD2F3C0A520}", 0, 1, 0)
application = kompas_api7_module.IApplication(Dispatch("Kompas.Application.7")._oleobj_.QueryInterface(kompas_api7_module.IApplication.CLSID, pythoncom.IID_IDispatch))
MH.iApplication = application
Documents = application.Documents
# Получим активный документ
kompas_document = application.ActiveDocument
kompas_document_3d = kompas_api7_module.IKompasDocument3D(kompas_document)
iDocument3D = kompas_object.ActiveDocument3D()
iPart = iDocument3D.GetPart(LDefin3D.pTop_Part)
api_kompas.connect()
a = 0
populations = []
populations.append(gen_alg.create_population(10,2))
for i in range(50):
parents = []
gen_alg.calc_pop(iPart,iDocument3D,populations[i])
gen_alg.calc_k_pop(populations[i])
parents = gen_alg.get_parent_for_new_population(populations[i])
print(parents)
populations.append(gen_alg.select_population(parents))
if i-1 < 50:
print(populations[i+1])
gen_alg.mutate(populations[i+1])До 35 строчки мы просто подключаемся к компасу, дальше создаем список популяций, добавляем первую популяцию, производим с ним продувки, считаем тягу с моментом, дальше коэффициенты выживаемости— подбираем родителей для нового поколения и создаем новое поколение и это всё повторяется несколько раз
Теперь вернёмся к тороидальным винтам.
Концепция будет следующей: в питоне генерируем точки с нужными координатами, перерисовываем в компас, обводим сплайнами, натягиваем поверхности, сшиваем — лопасти готовы.
Предлагаю начать с профилей. Для этого переписываем функцию на пайтон:
def naca4(number, chord=1, angle_of_attack=0, num_points=15, angle_z = 0,x0 =0 ,y0 = 0, z0 = 0, long = 1):
# Разбиваем число NACA на составляющие
m = int(number[0]) / 100.0
p = int(number[1]) / 10.0
t = int(number[2:]) / 100.0
# Угол атаки в радианах
alpha = np.radians(angle_of_attack)
beta = np.radians(angle_z)
# Создаем массив x от 0 до 1 с num_points точками
z = np.linspace(0, chord, num_points)
x =[0]
for i in range(num_points-2):
#if i<= 0.5*num_points:
x.append((1-np.cos(np.pi * z[i]))/2)
#else:
# x.append((1+(1-np.cos(np.pi*z[i])))/2)
x.append(1)
l = []
for i in x:
l.append(round(i, 2))
x = np.array(x)
# Толщина профиля
yt = 5 * t * chord * (
0.2969 * np.sqrt(x/chord)
- 0.1260 * (x/chord)
- 0.3516 * (x/chord)**2
+ 0.2843 * (x/chord)**3
- 0.1015 * (x/chord)**4
)
# Камбер линия
yc = np.where((x/chord) < p,
m * (x/chord) / p**2 * (2 * p - (x/chord)),
m * (1 - (x/chord)) / (1 - p)**2 * (1 + (x/chord) - 2 * p))
# Угол наклона камбер линии
dyc_dx = np.where((x/chord) < p,
2 * m / p**2 * (p - (x/chord)),
2 * m / (1 - p)**2 * (p - (x/chord)))
theta = np.arctan(dyc_dx)
# Координаты верхней и нижней поверхностей
xu = x - yt * np.sin(theta)
yu = yc + yt * np.cos(theta)
xl = x + yt * np.sin(theta)
yl = yc - yt * np.cos(theta)
return xu,yu,xl,ylТут просто генерируется профиль без преобразований и да, в последствии выяснилось, что на концах профиля нужно больше точек поэтому расстояние между точками идет по синусоиде, то есть в начале и конце их получается больше и это плавно переливается. Итак, нам надо добавить перемещение профиля, масштабирование, верчение и так далее. Тут, я думаю, тоже сильно ничего сложного нету просто вставлю полную функцию:
def naca4(number, chord=1, angle_of_attack=0, num_points=15, angle_z = 0,x0 =0 ,y0 = 0, z0 = 0, long = 1):
# Разбиваем число NACA на составляющие
m = int(number[0]) / 100.0
p = int(number[1]) / 10.0
t = int(number[2:]) / 100.0
# Угол атаки в радианах
alpha = np.radians(angle_of_attack)
beta = np.radians(angle_z)
# Создаем массив x от 0 до 1 с num_points точками
z = np.linspace(0, chord, num_points)
x =[0]
for i in range(num_points-2):
#if i<= 0.5*num_points:
x.append((1-np.cos(np.pi * z[i]))/2)
#else:
# x.append((1+(1-np.cos(np.pi*z[i])))/2)
x.append(1)
l = []
for i in x:
l.append(round(i, 2))
x = np.array(x)
# Толщина профиля
yt = 5 * t * chord * (
0.2969 * np.sqrt(x/chord)
- 0.1260 * (x/chord)
- 0.3516 * (x/chord)**2
+ 0.2843 * (x/chord)**3
- 0.1015 * (x/chord)**4
)
# Камбер линия
yc = np.where((x/chord) < p,
m * (x/chord) / p**2 * (2 * p - (x/chord)),
m * (1 - (x/chord)) / (1 - p)**2 * (1 + (x/chord) - 2 * p))
# Угол наклона камбер линии
dyc_dx = np.where((x/chord) < p,
2 * m / p**2 * (p - (x/chord)),
2 * m / (1 - p)**2 * (p - (x/chord)))
theta = np.arctan(dyc_dx)
# Координаты верхней и нижней поверхностей
xu = x - yt * np.sin(theta)
yu = yc + yt * np.cos(theta)
xl = x + yt * np.sin(theta)
yl = yc - yt * np.cos(theta)
# Поворачиваем профиль на угол атаки
xu_rot = xu * np.cos(alpha) - yu * np.sin(alpha)
yu_rot = xu * np.sin(alpha) + yu * np.cos(alpha)
xl_rot = xl * np.cos(alpha) - yl * np.sin(alpha)
yl_rot = xl * np.sin(alpha) + yl * np.cos(alpha)
#создаём 2 массива заполненые 0
zu_rot = []
zl_rot = []
for i in range(len(xu_rot)):
zu_rot.append(0)
zl_rot.append(0)
zu_rot = np.array(zu_rot)
zl_rot = np.array(zl_rot)
#маштабируем
xu_rot = xu_rot * long
xl_rot = xl_rot * long
yl_rot = yl_rot * long
yu_rot = yu_rot * long
zu_rot = zu_rot * long
zl_rot = zl_rot * long
x11 = xu_rot[0]
y11 = yu_rot[0]
z11 = 10
x12 = xu_rot[0] + 10
y12 = yu_rot[0]
z12 = 0
x13 = xu_rot[0]
y13 = yu_rot[0]+10
z13 = 0
x21 = xu_rot[-1]
y21 = yu_rot[-1]
z21 = 10
x22 = xu_rot[-1] + 10
y22 = yu_rot[-1]
z22 = 0
x23 = xu_rot[-1]
y23 = yu_rot[-1]+10
z23 = 0
#перемещаем профиль в x0,y0,z0
xu_rot = xu_rot + x0
xl_rot = xl_rot + x0
yl_rot = yl_rot + y0
yu_rot = yu_rot + y0
zu_rot = zu_rot + z0
zl_rot = zl_rot + z0
x11 = x11 + x0
y11 = y11 + y0
z11 = z11 + z0
x12 = x12 + x0
y12 = y12 + y0
z12 = z12 + z0
x13 = x13 + x0
y13 = y13 + y0
z13 = z13 + z0
x21 = x21 + x0
y21 = y21 + y0
z21 = z21 + z0
x22 = x22 + x0
y22 = y22 + y0
z22 = z22 + z0
x23 = x23 + x0
y23 = y23 + y0
z23 = z23 + z0
print(zl_rot)
#поворчиваем угол по оси z
xu_rot_z = xu_rot * np.cos(beta) - zu_rot*np.sin(beta)
zu_rot_z = xu_rot * np.sin(beta) + zu_rot*np.cos(beta)
xl_rot_z = xl_rot * np.cos(beta) - zu_rot*np.sin(beta)
zl_rot_z = xl_rot * np.sin(beta) + zl_rot*np.cos(beta)
x11_rot = x11* np.cos(beta) - z11*np.sin(beta)
z11_rot = x11* np.sin(beta) + z11*np.cos(beta)
x12_rot = x12* np.cos(beta) - z12*np.sin(beta)
z12_rot = x12* np.sin(beta) + z12*np.cos(beta)
x13_rot = x13* np.cos(beta) - z13*np.sin(beta)
z13_rot = x13* np.sin(beta) + z13*np.cos(beta)
x21_rot = x21* np.cos(beta) - z21*np.sin(beta)
z21_rot = x21* np.sin(beta) + z21*np.cos(beta)
x22_rot = x22* np.cos(beta) - z22*np.sin(beta)
z22_rot = x22* np.sin(beta) + z22*np.cos(beta)
x23_rot = x23* np.cos(beta) - z23*np.sin(beta)
z23_rot = x23* np.sin(beta) + z23*np.cos(beta)
#print(zu_rot_z)
#Отображаем профиль для отладки можно закоментировать
#plt.figure(figsize=(12, 6))
#plt.plot(xu_rot, yu_rot, 'b')
#plt.plot(xl_rot, yl_rot, 'b')
#plt.grid(True)
#plt.axis('equal')
#plt.show()
return xu_rot_z,xl_rot_z[1:15],yu_rot,yl_rot[1:15],zu_rot_z,zl_rot_z[1:15],x11_rot,y11,z11_rot,x12_rot,y12,z12_rot,x13_rot,y13,z13_rot,x21_rot,y21,z21_rot,x22_rot,y22,z22_rot,x23_rot,y23,z23_rotИтак, когда у нас профили, нам нужны кромки винта. Тут я сильно ничего не стал выдумывать — просто рисую точки в цилиндрической системе координат. Тут тоже без пояснений.
def point_edge(angle_z = 0, r = 1, angle_xy = 90, height = 1, x0 = 0, z0 = 0):
alpha = np.radians(angle_z)
beta = np.radians(angle_xy)
z = r
x = 0
z = z * np.cos(alpha) - x * np.sin(alpha)
x = z * np.sin(alpha) + x * np.cos(alpha)
z = z + z0
x = x + x0
#x = x * np.cos(beta) - z * np.sin(beta)
#z = x * np.sin(beta) + z * np.cos(beta)
y = height
return x,y,zТеперь нам надо нарисовать это всё в компасе: просто в цикле идем по точкам, читаем их координаты и рисуем:
def create_point_profl(xu_rot,xl_rot,yu_rot,yl_rot,zu_rot,zl_rot,kompas_document_3d):
kompas6_constants = gencache.EnsureModule("{75C9F5D0-B5B8-4526-8681-9903C567D2ED}", 0, 1, 0).constants
kompas6_constants_3d = gencache.EnsureModule("{2CAF168C-7961-4B90-9DA2-701419BEEFE3}", 0, 1, 0).constants
kompas_api7_module = gencache.EnsureModule("{69AC2981-37C0-4379-84FD-5DD2F3C0A520}", 0, 1, 0)
application = kompas_api7_module.IApplication(Dispatch("Kompas.Application.7")._oleobj_.QueryInterface(kompas_api7_module.IApplication.CLSID, pythoncom.IID_IDispatch))
MH.iApplication = application
for i in range(len(xu_rot)):
iPart7 = kompas_document_3d.TopPart
iModelContainer = iPart7._oleobj_.QueryInterface(kompas_api7_module.IModelContainer.CLSID, pythoncom.IID_IDispatch)
iModelContainer = kompas_api7_module.IModelContainer(iModelContainer)
iPoints3D = iModelContainer.Points3D
iPoint3D = iPoints3D.Add()
iPoint3D.ParameterType = kompas6_constants_3d.ksPParamCoord
iPoint3D.Symbol = kompas6_constants.ksDotPoint
iPoint3D.X = xu_rot[i]
iPoint3D.Y = yu_rot[i]
iPoint3D.Z = zu_rot[i]
iPoint3D.Update()
for i in range(len(xl_rot)):
iPart7 = kompas_document_3d.TopPart
iModelContainer = iPart7._oleobj_.QueryInterface(kompas_api7_module.IModelContainer.CLSID, pythoncom.IID_IDispatch)
iModelContainer = kompas_api7_module.IModelContainer(iModelContainer)
iPoints3D = iModelContainer.Points3D
iPoint3D = iPoints3D.Add()
iPoint3D.ParameterType = kompas6_constants_3d.ksPParamCoord
iPoint3D.Symbol = kompas6_constants.ksDotPoint
iPoint3D.X = xl_rot[i]
iPoint3D.Y = yl_rot[i]
iPoint3D.Z = zl_rot[i]
iPoint3D.Update()
def create_point_edge(x,y,z,kompas_document_3d):
kompas6_constants = gencache.EnsureModule("{75C9F5D0-B5B8-4526-8681-9903C567D2ED}", 0, 1, 0).constants
kompas6_constants_3d = gencache.EnsureModule("{2CAF168C-7961-4B90-9DA2-701419BEEFE3}", 0, 1, 0).constants
kompas_api7_module = gencache.EnsureModule("{69AC2981-37C0-4379-84FD-5DD2F3C0A520}", 0, 1, 0)
application = kompas_api7_module.IApplication(Dispatch("Kompas.Application.7")._oleobj_.QueryInterface(kompas_api7_module.IApplication.CLSID, pythoncom.IID_IDispatch))
MH.iApplication = application
iPart7 = kompas_document_3d.TopPart
iModelContainer = iPart7._oleobj_.QueryInterface(kompas_api7_module.IModelContainer.CLSID, pythoncom.IID_IDispatch)
iModelContainer = kompas_api7_module.IModelContainer(iModelContainer)
iPoints3D = iModelContainer.Points3D
iPoint3D = iPoints3D.Add()
iPoint3D.ParameterType = kompas6_constants_3d.ksPParamCoord
iPoint3D.Symbol = kompas6_constants.ksDotPoint
iPoint3D.X = x
iPoint3D.Y = y
iPoint3D.Z = z
iPoint3D.Update()Теперь нам надо один раз отрисовать их в компасе:
Нарисовали, соединили сплайнами. Отлично (не очень, я не нашел, как автоматически менять координаты точек, поэтому каждой координате каждой точки присваиваем глобальную переменную и делаем её внешней).
Итак, теперь функции для изменений переменных. Ничего сложного, просто взято из документации:
def edit_profl(iPart,num,xu,xl,yu,yl,zu,zl):
global VariableCollection
VariableCollection = iPart.VariableCollection()
for i in range(len(xu)):
VariableCollection.GetByName(f"px{num}1{i+1}",True,True).value = xu[i]
VariableCollection.GetByName(f"py{num}1{i+1}",True,True).value = yu[i]
VariableCollection.GetByName(f"pz{num}1{i+1}",True,True).value = zu[i]
for i in range(len(xl)):
VariableCollection.GetByName(f"px{num}2{i+1}",True,True).value = xl[i]
VariableCollection.GetByName(f"py{num}2{i+1}",True,True).value = yl[i]
VariableCollection.GetByName(f"pz{num}2{i+1}",True,True).value = zl[i]
print(f"pz{num}2{i+1}")
def edit_cross_points(iPart,num,x11,y11,z11,x12,y12,z12,x13,y13,z13,x21,y21,z21,x22,y22,z22,x23,y23,z23):
global VariableCollection
VariableCollection = iPart.VariableCollection()
VariableCollection.GetByName(f"xx{num}1",True,True).value = x11
VariableCollection.GetByName(f"xy{num}1",True,True).value = y11
VariableCollection.GetByName(f"xz{num}1",True,True).value = z11
VariableCollection.GetByName(f"yx{num}1",True,True).value = x12
VariableCollection.GetByName(f"yy{num}1",True,True).value = y12
VariableCollection.GetByName(f"yz{num}1",True,True).value = z12
VariableCollection.GetByName(f"zx{num}1",True,True).value = x13
VariableCollection.GetByName(f"zy{num}1",True,True).value = y13
VariableCollection.GetByName(f"zz{num}1",True,True).value = z13
VariableCollection.GetByName(f"xx{num}2",True,True).value = x21
VariableCollection.GetByName(f"xy{num}2",True,True).value = y21
VariableCollection.GetByName(f"xz{num}2",True,True).value = z21
VariableCollection.GetByName(f"yx{num}2",True,True).value = x22
VariableCollection.GetByName(f"yy{num}2",True,True).value = y22
VariableCollection.GetByName(f"yz{num}2",True,True).value = z22
VariableCollection.GetByName(f"zx{num}2",True,True).value = x23
VariableCollection.GetByName(f"zy{num}2",True,True).value = y23
VariableCollection.GetByName(f"zz{num}2",True,True).value = z23
def edit_edge_points(iPart,x1,y1,z1,x2,y2,z2,edge):
global VariableCollection
VariableCollection = iPart.VariableCollection()
VariableCollection.GetByName(f"{edge}x1",True,True).value = x1
VariableCollection.GetByName(f"{edge}y1",True,True).value = y1
VariableCollection.GetByName(f"{edge}z1",True,True).value = z1
VariableCollection.GetByName(f"{edge}x2",True,True).value = x2
print(f"{edge}x2")
VariableCollection.GetByName(f"{edge}y2",True,True).value = y2
VariableCollection.GetByName(f"{edge}z2",True,True).value = z2
Ну и небольшой main просто для прогонки винта:
import api_kompas
import api_fluid
import pythoncom
from win32com.client import Dispatch, gencache
import KompasAPI7
import LDefin3D
import MiscellaneousHelpers as MH
import exsel
import time
import points
iPart,iDocument = api_kompas.connect()
#xu_rot,xl_rot,yu_rot,yl_rot,zu_rot,zl_rot,x11,y11,z11,x12,y12,z12,x13,y13,z13,x21,y21,z21,x22,y22,z22,x23,y23,z23 = points.naca4("4812", angle_of_attack=-5,angle_z=-45,x0=-10,y0=0,z0=0,long=20)
#api_kompas.edit_profl(iPart,1,xu_rot,xl_rot,yu_rot,yl_rot,zu_rot,zl_rot)
#api_kompas.edit_cross_points(iPart,1,x11,y11,z11,x12,y12,z12,x13,y13,z13,x21,y21,z21,x22,y22,z22,x23,y23,z23)
#x1,y1,z1 = points.point_edge(45,20,90,3.75,xu_rot[0],zu_rot[0])
#x2,y2,z2 = points.point_edge(45,30,90,5,xu_rot[0],zu_rot[0])
#api_kompas.edit_edge_points(iPart,x1,y1,z1,x2,y2,z2,"kdf")
#x1,y1,z1 = points.point_edge(45,20,90,3.75,xu_rot[14],zu_rot[14])
#x2,y2,z2 = points.point_edge(45,30,90,5,xu_rot[14],zu_rot[14])
#api_kompas.edit_edge_points(iPart,x1,y1,z1,x2,y2,z2,"kdb")
#xu_rot,xl_rot,yu_rot,yl_rot,zu_rot,zl_rot,x11,y11,z11,x12,y12,z12,x13,y13,z13,x21,y21,z21,x22,y22,z22,x23,y23,z23 = points.naca4("4812", angle_of_attack=-5,angle_z=45,x0=-10,y0=30,z0=0,long=20)
#api_kompas.edit_profl(iPart,2,xu_rot,xl_rot,yu_rot,yl_rot,zu_rot,zl_rot)
#api_kompas.edit_cross_points(iPart,2,x11,y11,z11,x12,y12,z12,x13,y13,z13,x21,y21,z21,x22,y22,z22,x23,y23,z23)
#x1,y1,z1 = points.point_edge(-45,20,90,26.25,xu_rot[0],zu_rot[0])
#x2,y2,z2 = points.point_edge(-45,30,90,25,xu_rot[0],zu_rot[0])
#api_kompas.edit_edge_points(iPart,x1,y1,z1,x2,y2,z2,"kuf")
#x1,y1,z1 = points.point_edge(-45,20,90,26.25,xu_rot[14],zu_rot[14])
#x2,y2,z2 = points.point_edge(-45,30,90,25,xu_rot[14],zu_rot[14])
#print(x1,y1,z1,x2,y2,z2)
#xu_rot,xl_rot,yu_rot,yl_rot,zu_rot,zl_rot,x11,y11,z11,x12,y12,z12,x13,y13,z13,x21,y21,z21,x22,y22,z22,x23,y23,z23 = points.naca4("4812", angle_of_attack=-5,angle_z=0,x0=-10,y0=50,z0=15,long=20)
#api_kompas.edit_profl(iPart,3,xu_rot,xl_rot,zu_rot,zl_rot,yu_rot,yl_rot)
#api_kompas.edit_cross_points(iPart,3,x11,z11,y11,x12,z12,y12,x13,z13,y13,x21,z21,y21,x22,z22,y22,x23,z23,y23)
#api_kompas.edit_edge_points(iPart,x1,y1,z1,x2,y2,z2,"kub")
#api_kompas.rebuild(iPart,iDocument)
for q in range(0,45,5):
for w in range(0,45,5):
for e in range(0,15,3):
for r in range(0,15,3):
for t in range(0,45,5):
for y in range(0,45,5):
for u in range(15,30,3):
for i in range(15,30,3):
xu_rot,xl_rot,yu_rot,yl_rot,zu_rot,zl_rot,x11,y11,z11,x12,y12,z12,x13,y13,z13,x21,y21,z21,x22,y22,z22,x23,y23,z23 = points.naca4("4812", angle_of_attack=-q,angle_z=-w,x0=-10,y0=0,z0=0,long=20)
api_kompas.edit_profl(iPart,1,xu_rot,xl_rot,yu_rot,yl_rot,zu_rot,zl_rot)
api_kompas.edit_cross_points(iPart,1,x11,y11,z11,x12,y12,z12,x13,y13,z13,x21,y21,z21,x22,y22,z22,x23,y23,z23)
x1,y1,z1 = points.point_edge(w,20,90,e,xu_rot[0],zu_rot[0])
x2,y2,z2 = points.point_edge(w,30,90,r,xu_rot[0],zu_rot[0])
api_kompas.edit_edge_points(iPart,x1,y1,z1,x2,y2,z2,"kdf")
x1,y1,z1 = points.point_edge(w,20,90,e,xu_rot[14],zu_rot[14])
x2,y2,z2 = points.point_edge(w,30,90,r,xu_rot[14],zu_rot[14])
api_kompas.edit_edge_points(iPart,x1,y1,z1,x2,y2,z2,"kdb")
xu_rot,xl_rot,yu_rot,yl_rot,zu_rot,zl_rot,x11,y11,z11,x12,y12,z12,x13,y13,z13,x21,y21,z21,x22,y22,z22,x23,y23,z23 = points.naca4("4812", angle_of_attack=-t,angle_z=y,x0=-10,y0=30,z0=0,long=20)
api_kompas.edit_profl(iPart,2,xu_rot,xl_rot,yu_rot,yl_rot,zu_rot,zl_rot)
api_kompas.edit_cross_points(iPart,2,x11,y11,z11,x12,y12,z12,x13,y13,z13,x21,y21,z21,x22,y22,z22,x23,y23,z23)
x1,y1,z1 = points.point_edge(-y,20,90,u,xu_rot[0],zu_rot[0])
x2,y2,z2 = points.point_edge(-y,30,90,i,xu_rot[0],zu_rot[0])
api_kompas.edit_edge_points(iPart,x1,y1,z1,x2,y2,z2,"kuf")
x1,y1,z1 = points.point_edge(-y,20,90,u,xu_rot[14],zu_rot[14])
x2,y2,z2 = points.point_edge(-y,30,90,i,xu_rot[14],zu_rot[14])
xu_rot,xl_rot,yu_rot,yl_rot,zu_rot,zl_rot,x11,y11,z11,x12,y12,z12,x13,y13,z13,x21,y21,z21,x22,y22,z22,x23,y23,z23 = points.naca4("4812", angle_of_attack=-5,angle_z=0,x0=-10,y0=50,z0=15,long=20)
api_kompas.edit_profl(iPart,3,xu_rot,xl_rot,zu_rot,zl_rot,yu_rot,yl_rot)
api_kompas.edit_cross_points(iPart,3,x11,z11,y11,x12,z12,y12,x13,z13,y13,x21,z21,y21,x22,z22,y22,x23,z23,y23)
api_kompas.edit_edge_points(iPart,x1,y1,z1,x2,y2,z2,"kub")
api_kompas.rebuild(iPart,iDocument)После выполнения немного измененного main получил табличку в exsel, оттуда достал самый лучший винт по отношению тяга/момент и напечатал на фотополимере:
И так далее, напечатав такой вот винт, я пришел к небольшому удивлению, что fluid ошибся по тяге всего на 2 грамма и что такой винт тянет в 2 раза меньше стандартного трехлопастного винта. По шумности он примерно такой же, только частота шума ниже.
Сейчас хочу отказать от компаса, сразу делать много точек и между ними делать полигоны и сохранять в stl, а так в планах сделать, что питон с нуля строил модель, но идей по этому поводу пока нет.
Автор: nbvtim
