Исключение: поверхность 1 не может быть построена с использованием трансфинитного алгоритма (несовпадающее количество узлов на противоположных сторонах 1524 != 1473 или 99 != 99)
Ошибка: поверхность 1 не может быть построена с помощью трансфинитного алгоритма (несовпадающее количество узлов на противоположных сторонах 1524 != 1473) или 99 != 99)
Я понятия не имею, что происходит! Я только начинаю работать с gmsh, а с python все становится гораздо более абстрактно. Если кто-нибудь знает, как ее решить, буду бесконечно рад!
Код: Выделить всё
import gmsh
import sys
import numpy as np
import matplotlib.pyplot as plt
def neg(lst):
return [ -i for i in lst ]
def getS(n,r):
if r!=1:
s=(r-1)/(r**n-1)
elif r==1:
s=1/n
return s
def getPoints(start,stop,n,r):
s=getS(n,r)
x=np.zeros(n+1)
x[0]=start
x[1:]=start+(stop-start)*s*np.cumsum(np.array([r**i for i in range(n)]))
return x
def chamber(Rc,step_count,Lc,progression):
x = getPoints(0,Lc,step_count,progression)#np.linspace(0,Lc,step_count)
y = Rc*np.ones(len(x))
return [x,y]
def chamber_exit_circ(Rc,Rt,Lc,theta_t_i,step_count):
chamber_e_l = Lc
angles = np.linspace(90, 90 - theta_t_i, step_count+1)
Xc = chamber_e_l
Yc = Rc -1.5*Rt
x = 1.5*Rt*np.cos(angles*np.pi/180) + Xc
y = 1.5*Rt*np.sin(angles*np.pi/180) + Yc
return [x,y]
def throat_inlet(Rt,theta_t_i,step_count,progression):
throat_i_l = 1.5*Rt*np.sin(theta_t_i*np.pi/180)
angles = getPoints(270 - theta_t_i, 270, step_count+1,progression)
Xc = throat_i_l
Yc = 2.5*Rt;
x = 1.5*Rt*np.cos(angles*np.pi/180) + Xc
y = 1.5*Rt*np.sin(angles*np.pi/180) + Yc
return [x,y]
def throat_exit(Rt,theta_N,step_count):
throat_e_l = 0.382*Rt*np.sin(theta_N*np.pi/180)
angles = np.linspace(270 , 270 + theta_N ,step_count+1)
Xc = throat_e_l
Yc = 1.382*Rt
x = 0.382*Rt*np.cos(angles*np.pi/180) + Xc
y = 0.382*Rt*np.sin(angles*np.pi/180) + Yc
return[x,y]
def bell_curve(Rt,Re,theta_N,theta_e,Ln_ratio,step_count,progression):
# bezler quadratic curve equation %
expansion_ratio = (Re/Rt)**2
Nx = 0.382*Rt*np.cos((theta_N - 90)*np.pi/180)
Ny = 1.382*Rt + 0.382*Rt*np.sin((theta_N - 90)*np.pi/180)
Ex = Ln_ratio*( (np.sqrt(expansion_ratio) -1)*Rt + 1.5*Rt*(1/np.cos(15*np.pi/180) -1))/np.tan(15*np.pi/180) # Ln
Ey = Re
m1 = np.tan(theta_N*np.pi/180)
m2 = np.tan(theta_e*np.pi/180)
C1 = Ny - m1*Nx
C2 = Ey - m2*Ex
Qx = (C2 - C1)/(m1-m2)
Qy = (m1*C2 -m2*C1)/(m1-m2)
x = np.zeros(step_count+1)
y = np.zeros(step_count+1)
i=0
for t in getPoints(0,1,step_count,progression):#np.linspace(0,1,step_count):
x[i] = (1-t)**2*Nx + 2*(1-t)*t*Qx + t**2*Ex
y[i] = (1-t)**2*Ny + 2*(1-t)*t*Qy + t**2*Ey
i=i+1
return [x,y]
# Rocket Dimensions
Rc = 71 # [mm]
Rt = 35.5 # [mm]
Re = 275 # [mm]
Lc = 100 # [mm]
Ln_ratio = 0.8; # Ratio of conical nozzle
theta_t_i = 30 # [deg]
theta_N = 30 # [deg]
theta_e = 7 # [deg]
step_count_c_e = 40 #20;
step_count_t_i = 60
step_count_t_e = 60 #40
step_count_r = 100 #80
progression_i=1/1.1
progression_bell=1.01
progression_r=1/1.05
# throat exit section
[x_t_e,y_t_e]=throat_exit(Rt,theta_N,step_count_t_e)
# throat inlet section
deltaN_t_i=x_t_e[1]-x_t_e[0]
throat_i_l = 1.5*Rt*np.sin(theta_t_i*np.pi/180)
delta0_t_i=deltaN_t_i*progression_i-throat_i_l *(progression_i-1)
step_count_t_i=int(np.ceil(np.log(deltaN_t_i/delta0_t_i)/np.log(progression_i)))
[x_t_i,y_t_i]=throat_inlet(Rt,theta_t_i,step_count_t_i,progression_i)
# chamber exit section
[x_c_e,y_c_e]=chamber_exit_circ(Rc,Rt,Lc,theta_t_i,step_count_c_e)
delta0_c_e=x_c_e[1]-x_c_e[0]
# chamber exit to throat inlet line
deltaN_ce_ti=x_t_i[1]-x_t_i[0]
L_c_e_ti=np.abs(np.tan((90 + theta_t_i)*np.pi/180)*(y_c_e[-1]-y_t_i[0]))
delta0_ce_ti=x_c_e[-1]-x_c_e[-2]
progression_ce_ti=(L_c_e_ti-delta0_ce_ti)/(L_c_e_ti-deltaN_ce_ti)
step_count_ce_ti=int(np.ceil(np.log(deltaN_ce_ti/delta0_ce_ti)/np.log(progression_ce_ti)))
y_ce_ti = getPoints(y_c_e[-1],y_t_i[0],step_count_ce_ti,progression_ce_ti)
x_ce_ti = np.tan((90 + theta_t_i)*np.pi/180)*y_ce_ti
# chamber section
step_count_c=int(np.ceil(np.log(Lc*(1/progression_i-1)/delta0_c_e+1)/np.log(1/progression_i)))-1
[x_c,y_c]=chamber(Rc,step_count_c,Lc,progression_i)
delta0_c=Lc*getS(step_count_c,progression_i)
## chamber exit section (translate)
x_c_e = x_c_e + x_c[-1] - x_c_e[0]
#translate chamber exit to throat inlet, throat inlet and throat exit section sections
x_ce_ti = x_ce_ti + x_c_e[-1] - x_ce_ti[0]
x_t_i = x_t_i + x_ce_ti[-1] - x_t_i[0]
x_t_e = x_t_e + x_t_i[-1] - x_t_e[0]
# bell section
expansion_ratio = (Re/Rt)**2
delta0_t_e=x_t_e[-1]-x_t_e[-2]
Ln=Ln_ratio*( (np.sqrt(expansion_ratio) -1)*Rt + 1.5*Rt*(1/np.cos(15*np.pi/180) -1))/np.tan(15*np.pi/180)
step_count_nozzle=int(np.ceil(np.log(Ln*(progression_bell-1)/delta0_t_e+1)/np.log(progression_bell)))-1
[x_bell,y_bell]=bell_curve(Rt,Re,theta_N,theta_e,Ln_ratio,step_count_nozzle,progression_bell)
#translate bell section
x_bell = x_bell + x_t_e[-1] - x_bell[0]
gmsh.initialize()
gmsh.model.add("RaoNozzle")
factory=gmsh.model.occ
lc = 1
#pos=np.loadtxt('/Users/jsalazar/matlab/RAO-Bell-Nozzle/nozzleContour.csv',delimiter=',', skiprows=1)
x=np.concatenate((x_c,x_c_e[1:],x_ce_ti[1:],x_t_i[1:],x_t_e[1:],x_bell[1:]))
y=np.concatenate((y_c,y_c_e[1:],y_ce_ti[1:],y_t_i[1:],y_t_e[1:],y_bell[1:]))
pos=np.concatenate((x,y))
pos=pos.reshape((len(x),2),order='F')
npts=len(pos)
erase=[]
for i in range(0,npts-1):
if np.linalg.norm(pos[i]-pos[i+1])==0:
erase.append(i)
pos=np.delete(pos,erase,0)
npts=len(pos)
np.savetxt('nozzleContour.csv', pos/1e3,delimiter=',', header='x [m] y [m]')
pAxis=[[]]*npts
pBell=[[]]*npts
lAxis=[[]]*(npts-1)
lBell=[[]]*(npts-1)
lRadial=[[]]*2
for i in range(0,npts):
pBell[i]=factory.addPoint(pos[i,0],pos[i,1],0,lc)
pAxis[i]=factory.addPoint(pos[i,0],0,0,lc)
for i in range(0,npts-1):
lBell[i]=factory.addLine(pBell[i],pBell[i+1])
lAxis[i]=factory.addLine(pAxis[i],pAxis[i+1])
j=0
for i in [0,npts-1]:
print(i)
lRadial[j]=factory.addLine(pAxis[i],pBell[i])
factory.synchronize()
gmsh.model.mesh.setTransfiniteCurve(lRadial[j], step_count_r, meshType="Progression", coef=progression_r)
j=j+1
lAxis.reverse()
loop=lBell+neg([lRadial[1]])+neg(lAxis)+[lRadial[0]]
cLoop=[]
cLoop.append(factory.addCurveLoop(loop))
surface=[]
surface.append(factory.addPlaneSurface(cLoop))
factory.synchronize()
gmsh.model.mesh.setRecombine(2, surface[0])
factory.synchronize()
gmsh.model.mesh.setTransfiniteSurface(surface[0],cornerTags=[pAxis[0],pAxis[-1],pBell[0],pBell[-1]])
factory.synchronize()
v=factory.extrude([(2,surface[0])],0,0,1,numElements=[1],recombine=True)
factory.synchronize()
gmsh.model.addPhysicalGroup(2,[surface[0],v[0][1]], name="frontAndBack")
gmsh.model.addPhysicalGroup(2,[v[-1][1]], name="inlet")
gmsh.model.addPhysicalGroup(2,[x[1] for x in v[2:npts+1]], name="wall")
gmsh.model.addPhysicalGroup(2,[v[npts+1][1]], name="outlet")
gmsh.model.addPhysicalGroup(2,[x[1] for x in v[npts+2:2*npts+1]], name="axis")
gmsh.model.addPhysicalGroup(3,[v[1][1]], name="internal")
factory.synchronize()
gmsh.model.mesh.generate(3)
gmsh.write("raoNozzle.msh2")
if '-nopopup' not in sys.argv:
gmsh.fltk.run()
gmsh.finalize()
Подробнее здесь: https://stackoverflow.com/questions/788 ... ve-command
