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