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