Selecting the optimal geometry in compass 3d using fluid x3d
Start
It all started when I was given an interesting project to develop toroidal propellers.
It all started in Solidworks with the modeling of a standard two-bladed propeller with variable attack angles.
Next, the selection of parameters for launching the study began:
At first, we used a periodic calculation domain, since we read somewhere that it is necessary to use it for blowing the propellers, but after a week of active bewilderment about why the thrust changed when the calculation domain changed, we understood how the periodic calculation domain works.
That is, the flows of the calculation area are copied to neighboring cells with the same periodicity as the calculation area, and yes, all this is done in the rotation of the local Averaging area – this is where the model rotates in the grid, otherwise all the same settings as in all guides. There was also a misunderstanding why a 5-inch propeller at 25,000 rpm creates a thrust of 0.001 N.
Then, simply for scientific purposes, I set Sliding in the rotation of the local area – in it, the grid rotates around the model, that is, a grid cylinder is taken and begins to rotate in the grid. And, lo and behold, plausible results were obtained.
At the next stage, it became interesting what would happen if we tweaked the grid resolution, because the “Purge” time varied from a couple of seconds to a week on a not very weak PC with intel 9 10 940 and quadro a5000. For this, the screw was scanned with a shining 3d ue pro scanner. By the way, it has a scanning error of 0.01 mm and reverse engineering was performed to optimize the geometry. As a sign of the end of the purge, we looked at the thrust graph from the iterations and stopped the “purge” when the thrust reached the “shelf”.
The propeller thrust on the test bench was 6.75N and then questions arose about solidworks, in addition to this, it was said to do all this on domestic software, because of this we switched to compass, but its built-in cfd calculation module does not allow “blowing” moving parts. It was decided to use the open source project fluid x3d, which, unlike solidworks, performs calculations on a video card instead of a processor, which reduces the calculation time from a week to several minutes! So, after the first fluid tests, the following settings were written for blowing the propeller and exporting the blowing values to a txt file.
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 – link to fluid, if anyone is interested, dig around, I recommend
Now we need to import the blowdown data into python because I don't know how to program in anything else
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")Now to the compass. There are two options. Let's start with the simpler one.
We draw a standard two-bladed propeller by sections, and in the sketches we make angles that we take out into global variables.
Now it's time to work on the compass 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")There is nothing complicated here, just connect to the compass, change two variables and save/delete stl.
Now for complete happiness we need a genetic algorithm:
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)
All functions are labeled so that everything should be clear. I do not argue that it is not ideal – everything was done the night before the performance 🙂
So, now let's make a simple main for the genetic algorithm:
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])Up to line 35 we simply connect to the compass, then we create a list of populations, add the first population, perform blowdowns with it, calculate the thrust with the moment, then the survival rates—we select parents for the new generation and create a new generation, and all this is repeated several times
Now let's get back to toroidal screws.
The concept will be as follows: in Python we generate points with the required coordinates, redraw them in a compass, outline them with splines, stretch the surfaces, stitch them together – the blades are ready.
I suggest starting with profiles. To do this, rewrite the function in Python:
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,ylHere, a profile is simply generated without transformations and yes, later it turned out that more points are needed at the ends of the profile, so the distance between the points is sinusoidal, that is, at the beginning and end, there are more of them and this flows smoothly. So, we need to add profile movement, scaling, rotation, and so on. Here, I think, there is also nothing very complicated, I will simply insert the full function:
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_rotSo, when we have profiles, we need screw edges. Here I didn't invent anything much – I just draw points in a cylindrical coordinate system. Here too, no explanations.
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,zNow we need to draw all this in the compass: we simply go through the points in a loop, read their coordinates and draw:
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()Now we need to draw them once in the compass:
Drawn, connected with splines. Excellent (not so good, I didn't find how to automatically change the coordinates of the points, so we assign a global variable to each coordinate of each point and make it external).
So, now the functions for changing variables. Nothing complicated, just taken from the documentation:
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
Well, and a small main just to run the screw:
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)After making a slightly modified main, I got a table in Excel, from there I got the best screw in terms of thrust/torque ratio and printed it on photopolymer:
And so on, having printed such a propeller, I came to a slight surprise that fluid was mistaken in the thrust by only 2 grams and that such a propeller pulls 2 times less than a standard three-blade propeller. In terms of noise, it is about the same, only the noise frequency is lower.
Now I want to abandon the compass, make many points at once and make polygons between them and save in stl, and so in the plans to make Python build the model from scratch, but there are no ideas about this yet.
