Начало
Всё началось с того, что мне подкинули интересный проект с разработкой тороидальных винтов.
Всё началось в 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", "# t\tCd\n"); //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, а так в планах сделать, что питон с нуля строил модель, но идей по этому поводу пока нет.
ссылка на оригинал статьи https://habr.com/ru/articles/839136/
Добавить комментарий