Note
Go to the end to download the full example code.
Mesh10#
Simplified turbine mesh with data extraction in matlab.
============= Init GMSH interface =============
Link contours (55.733 ms)
Link contours (36.718 ms)
Link contours (45.157 ms)
Link contours (45.753 ms)
Link contours (63.360 ms)
Link contours (59.716 ms)
Link contours (62.979 ms)
Link contours (60.912 ms)
Link contours (51.439 ms)
Link contours (51.850 ms)
Meshing with gmsh (13.261 ms)
Saving .msh (6.070 ms)
Construct mesh object (1.545 ms)
Element type: PRISM6
Ne = 512, Nn = 539
12 from EasyFEA import Display, Folder, np, Mesher, ElemType, PyVista
13 from EasyFEA.Geoms import Point, Points, Contour, CircleArc, Line
14
15 import scipy.io
16
17 if __name__ == "__main__":
18 Display.Clear()
19
20 # ----------------------------------------------
21 # Contour
22 # ----------------------------------------------
23
24 # outputs
25 folder = Folder.Join(Folder.RESULTS_DIR, "Meshes", "Mesh10")
26 addCylinder = True
27 repeat = False
28
29 # geom
30 angleRev = 2 * np.pi / 20 # rad
31 saveToMatlab = False
32
33 # mesh
34 N = 4 # elements in the blade lenght l
35 elemType = ElemType.PRISM6
36 mesher = Mesher(False, True, True)
37 factory = mesher._factory
38
39 # ----------------------------------------------
40 # Geom
41 # ----------------------------------------------
42
43 l = 2
44 t = 0.3
45 H1 = 5
46 c1 = 1
47 H2 = 5
48 c2 = 2
49 H3 = 5
50 c3 = 2
51 R = 5
52 e = 1
53
54 rot = 15
55 coefRot = rot / (H1 + H2 + H3) # deg
56 rot0 = 0
57 rot1 = H1 * coefRot
58 rot2 = (H1 + H2) * coefRot - rot1
59 rot3 = (H1 + H2 + H3) * coefRot - rot2
60 assert t < R
61
62 mS = l / N
63
64 vols: list[tuple[Contour, Contour]] = []
65
66 def bladeSection(l: float, t: float, center: np.ndarray, angle=0.0, R: float = 0.0):
67 """Create a blade section
68
69 Parameters
70 ----------
71 l : float
72 blade length
73 t : float
74 blade thickness
75 center : np.ndarray
76 center of mass for the blade section
77 angle : float, optional
78 rotation in deg along z axis (0,0,1), by default 0.0
79 R : float, optional
80 radius used to project on the cylinder, by default 0.0
81
82 Returns
83 -------
84 Contour
85 created contour
86 """
87
88 # crete points coordinates
89 n = 11
90 y = np.linspace(0, l, n)
91 y = np.append(y, y[::-1])
92
93 mat = np.array([[0, 0, 1], [(l / 2) ** 2, l / 2, 1], [l**2, l, 1]])
94 a, b, c = tuple(np.linalg.inv(mat) @ [0, -t / 2, 0])
95 x: np.ndarray = a * y**2 + b * y + c
96
97 x[y.size // 2 :] *= -1
98
99 z = np.zeros_like(y)
100
101 coord = np.array([x, y, z]).T
102
103 # first center the surface on (0,0,0) to do the rotation
104 coord = coord - np.array([0, l / 2, 0])
105
106 # rotate along z axis
107 rot = angle * np.pi / 180
108 rotMat = np.array(
109 [[np.cos(rot), -np.sin(rot), 0], [np.sin(rot), np.cos(rot), 0], [0, 0, 1]]
110 )
111 coord = coord @ rotMat
112
113 # project on R radius the z coordinates
114 if R != 0:
115 coord[:, 2] = -(R - np.sqrt(R**2 - coord[:, 0] ** 2))
116
117 # center the surface
118 coord += center
119
120 # create points
121 points = [Point(*co) for co in coord]
122
123 # split points in 4 groups to construct 4 points list
124 split1, split2 = np.reshape(points, (2, -1))
125
126 idxM = len(split1) // 2
127 spline1 = Points(split1[: idxM + 1], mS)
128 spline2 = Points(split1[idxM:], mS)
129 spline3 = Points(split2[: idxM + 1], mS)
130 spline4 = Points(split2[idxM:], mS)
131
132 # construct contour
133 contour = Contour([spline1, spline2, spline3, spline4])
134
135 # contour.Plot()
136
137 return contour, np.asarray(coord)
138
139 def Cylinder(
140 R, angleRev: float, coord: np.ndarray, y1: float
141 ) -> tuple[Contour, Contour]:
142 s = np.sin(angleRev / 2)
143 c = np.cos(angleRev / 2)
144
145 P1 = Point(-s * R, coord[0, 1], c * R)
146 P2 = Point(*coord[0])
147 P3 = Point(s * R, P2.y, c * R)
148 P4 = Point(s * R, y1, c * R)
149 P5 = Point(P2.x, y1, P2.z)
150 P6 = Point(-s * R, y1, c * R)
151
152 C1 = Point(y=P2.y)
153 C2 = Point(y=y1)
154
155 L1_l = CircleArc(P1, P2, C1, meshSize=mS)
156 L2_l = Line(P2, P5, mS)
157 L3_l = CircleArc(P5, P6, C2, meshSize=mS)
158 L4_l = Line(P6, P1, mS)
159
160 contour_l = Contour([L1_l, L2_l, L3_l, L4_l])
161
162 L1_r = CircleArc(P2, P3, C1, meshSize=mS)
163 L2_r = Line(P3, P4, mS)
164 L3_r = CircleArc(P4, P5, C2, meshSize=mS)
165 L4_r = Line(P5, P2, mS)
166
167 contour_r = Contour([L1_r, L2_r, L3_r, L4_r])
168
169 return contour_l, contour_r
170
171 def CylinderBlade(R, angleRev: float, coord: np.ndarray) -> Contour:
172 s = np.sin(angleRev / 2)
173 c = np.cos(angleRev / 2)
174
175 y0 = coord[0, 1]
176 y1 = coord[-1, 1]
177
178 P1 = Point(-s * R, y0, c * R)
179 P2 = Point(*coord[0])
180 P3 = Point(*coord[-1])
181 P4 = Point(-s * R, y1, c * R)
182
183 C1 = Point(0, y0, 0)
184 C2 = Point(0, y1, 0)
185
186 L1 = CircleArc(P1, P2, C1, meshSize=mS)
187 L2 = Points([Point(*co) for co in coord], mS)
188 L3 = CircleArc(P3, P4, C2, meshSize=mS)
189 L4 = Line(P4, P1, mS)
190
191 contour = Contour([L1, L2, L3, L4])
192
193 return contour
194
195 def LinkEveryone(
196 vols: list[tuple[Contour, Contour, list[int]]], rot=0.0
197 ) -> list[tuple]:
198 """contour1, contour2, nLayers, numElems
199 return created entities"""
200 allEntities = []
201 for vol in vols:
202 contour1, contour2, nLayers, numElems = vol
203 contour1 = contour1.copy()
204 contour1.Rotate(rot, direction=(0, 1, 0))
205 contour2 = contour2.copy()
206 contour2.Rotate(rot, direction=(0, 1, 0))
207 ents = mesher._Link_Contours(
208 contour1, contour2, elemType, nLayers, numElems
209 )
210
211 allEntities.extend(ents)
212
213 return allEntities
214
215 # ----------------------------------------------
216 # Mesh
217 # ----------------------------------------------
218
219 bladeInf, coordInf = bladeSection(l, t, (0, 0, R), rot0, R) # blade in z = R
220 bladeSup, coordSup = bladeSection(
221 l, t, (0, 0, R + e), rot0, R + e
222 ) # blade in z = R+e
223
224 coordInfL, coordInfR = np.reshape(coordInf, (2, -1, 3))
225 coordSupL, coordSupR = np.reshape(coordSup, (2, -1, 3))
226
227 blade1, __ = bladeSection(l * c1, t, (0, -(c1 * l - l) / 2, R + e + H1), rot1)
228 blade2, __ = bladeSection(l * c2, t, (0, -(c2 * l - l) / 2, R + e + H1 + H2), rot2)
229 blade3, __ = bladeSection(
230 l * c3, t, (0, -(c3 * l - l) / 2, R + e + H1 + H2 + H3), rot3
231 )
232
233 elems = [N / 2] * 4 # number of elements for each lines for blades
234
235 # contour1, contour2, nLayers, numElems
236 vols.append((bladeInf, bladeSup, e / mS, elems))
237 vols.append((bladeSup, blade1, H1 / mS, elems))
238 vols.append((blade1, blade2, H2 / mS, elems))
239 vols.append((blade2, blade3, H3 / mS, elems))
240
241 if addCylinder:
242 elems = [N, l / 2 / mS] * 2
243
244 # cylinder y=-l
245 cylinInf1_l, cylinInf1_r = Cylinder(R, angleRev, coordInfL, -l)
246 cylinSup1_l, cylinSup1_r = Cylinder(R + e, angleRev, coordSupL, -l)
247
248 vols.append((cylinInf1_l, cylinSup1_l, e / mS, elems))
249 vols.append((cylinInf1_r, cylinSup1_r, e / mS, elems))
250
251 # cylinder y=l
252 cylinInf2_l, cylinInf2_r = Cylinder(R, angleRev, coordInfR, l)
253 cylinSup2_l, cylinSup2_r = Cylinder(R + e, angleRev, coordSupR, l)
254 vols.append((cylinInf2_l, cylinSup2_l, e / mS, elems))
255 vols.append((cylinInf2_r, cylinSup2_r, e / mS, elems))
256
257 elems = [N] * 4
258 # cylinder blade left
259 coordInfL = coordInf[: coordInf.shape[0] // 2]
260 coordSupfL = coordSup[: coordSup.shape[0] // 2]
261 contourInf = CylinderBlade(R, angleRev, coordInfL)
262 contourSup = CylinderBlade(R + e, angleRev, coordSupfL)
263
264 # cylinder blade right
265 coordInfR = coordInf[coordInf.shape[0] // 2 :]
266 coordSupfR = coordSup[coordSup.shape[0] // 2 :]
267 contourInf_c = CylinderBlade(R, -angleRev, coordInfR[::-1])
268 contourSup_c = CylinderBlade(R + e, -angleRev, coordSupfR[::-1])
269
270 vols.append((contourInf, contourSup, e / mS, elems))
271 vols.append((contourInf_c, contourSup_c, e / mS, elems))
272
273 # axGeom = Display.init_Axes(3)
274 # axGeom.axis('off')
275 # for v, vol in enumerate(vols):
276 # color = 'blue'
277 # vol[0].Plot(axGeom, color=color)
278 # vol[1].Plot(axGeom, color=color)
279 # # ax.legend()
280
281 firstPart = LinkEveryone(vols)
282
283 mesher._Set_mesh_algorithm(elemType)
284
285 # mesher._synchronize()
286 # parts = factory.getEntities()
287
288 if repeat:
289 na = int(np.pi * 2 / angleRev)
290
291 for i in range(na - 1):
292 print(f"{(i + 1) / (na - 1) * 100:2.0f} %")
293 LinkEveryone(vols, angleRev * (i + 1) * 180 / np.pi)
294
295 # # dont work
296 # partsC = factory.copy(parts)
297 # factory.rotate(partsC, 0,0,0,0,0,1, angleRev*(i+1))
298
299 mesher._Mesh_Generate(3, elemType, folder=folder, filename="blade")
300
301 mesh = mesher._Mesh_Get_Mesh()
302
303 nodesCircle = mesh.Nodes_Conditions(
304 lambda x, y, z: np.sqrt(x**2 + z**2) <= R + 1e-2
305 )
306 nodesUpper = mesh.Nodes_Conditions(
307 lambda x, y, z: z >= mesh.coord[:, 2].max() - 1e-2
308 )
309
310 nodesBlades = mesh.Nodes_Conditions(
311 lambda x, y, z: np.sqrt(x**2 + z**2) >= R + e + 1e-1
312 )
313 nodesCyl = mesh.Nodes_Conditions(lambda x, y, z: np.sqrt(x**2 + z**2) <= R + e)
314
315 mesh.Set_Tag(nodesCircle, "blades")
316 mesh.Set_Tag(nodesCyl, "cylindre")
317
318 if saveToMatlab:
319 matFile = Folder.Join(folder, "blade.mat")
320 msh = {
321 "connect": np.asarray(mesh.connect + 1, dtype=float),
322 "coord": mesh.coord,
323 "bladesNodes": np.asarray(nodesBlades + 1, dtype=float),
324 "bladesElems": np.asarray(
325 mesh.Elements_Nodes(nodesBlades) + 1, dtype=float
326 ),
327 "cylindreNodes": np.asarray(nodesCyl + 1, dtype=float),
328 "cylindreElems": np.asarray(mesh.Elements_Nodes(nodesCyl) + 1, dtype=float),
329 }
330 scipy.io.savemat(matFile, msh)
331
332 # ----------------------------------------------
333 # Plot
334 # ----------------------------------------------
335
336 PyVista.Plot_Mesh(mesh).show()
337
338 qual = mesh.Get_Quality("aspect")
339 PyVista.Plot(
340 mesh,
341 qual,
342 nodeValues=False,
343 plotMesh=True,
344 cmap="viridis",
345 clim=(0, 1),
346 nColors=11,
347 ).show()
348
349 Display.plt.show()
Total running time of the script: (0 minutes 1.120 seconds)