geometricallyExactBeamTest.py
You can view and download this file on Github: geometricallyExactBeamTest.py
1#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2# This is an EXUDYN example
3#
4# Details: Test models for GeometricallyExactBeam (2-node shear deformable beam,
5# Lie group formulation for work of elastic forces);
6# test models: cantilever beam with tip force and torque
7#
8# Author: Johannes Gerstmayr
9# Date: 2023-04-05
10#
11# Copyright:This file is part of Exudyn. Exudyn is free software. You can redistribute it and/or modify it under the terms of the Exudyn license. See 'LICENSE.txt' for more details.
12#
13#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
14
15import exudyn as exu
16from exudyn.utilities import *
17
18import numpy as np
19
20useGraphics = True #without test
21#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
22#you can erase the following lines and all exudynTestGlobals related operations if this is not intended to be used as TestModel:
23try: #only if called from test suite
24 from modelUnitTests import exudynTestGlobals #for globally storing test results
25 useGraphics = exudynTestGlobals.useGraphics
26except:
27 class ExudynTestGlobals:
28 pass
29 exudynTestGlobals = ExudynTestGlobals()
30#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
31
32SC = exu.SystemContainer()
33mbs = SC.AddSystem()
34
35compute2D = False
36compute3D = True
37
38#test examples
39#2011 MUBO, Nachbagauer Pechstein Irschik Gerstmayr (2D)
40#2013 CND, Nachbagauer Gruber Gerstmayr (static, 3D); "Structural and Continuum Mechanics Approaches for a 3D Shear Deformable ANCF Beam Finite Element: Application to Static and Linearized Dynamic Examples"
41### not yet: 2013 CND, Nachbagauer Gerstmayr (dynamic, 3D)
42cases = ['CantileverLinear2011', 'Cantilever2011', 'GeneralBending2013', 'PrincetonBeamF2', 'PrincetonBeamF3', 'Eigenmodes2013']
43nElementsList = [1,2,4,8,16,32,64,128,256,512,1024]
44# nElementsList = [8,32, 128]
45nElements = 8
46
47betaList = [0,15,30,45,60,75,90]
48betaDegree = 45
49
50caseList = [0,1,2,3,4] #case 0 not working for Geometrically exact beam
51#case=2
52
53useGraphics = False
54verbose = 1*0
55
56useEP = True #for geometrically exact beam node
57if useEP:
58 NodeClass = NodeRigidBodyEP
59 initialRotationsGE = eulerParameters0
60else: #does not work for static case, as static solver currently (2023-04) cannot solve for Lie group nodes
61 NodeClass = NodeRigidBodyRotVecLG
62 initialRotationsGE = [0,0,0]
63
64
65bodyFixedLoad = False
66testErrorSum = 0
67case=4
68
69printCase = True
70#for nElements in nElementsList:
71#for betaDegree in betaList:
72for case in caseList:
73#if True:
74 if printCase:
75 printCase=False
76 exu.Print('case=', case, cases[case])
77 mbs.Reset()
78
79 computeEigenmodes = False
80 csFact = 1
81 sectionData = exu.BeamSection()
82 fTip = 0
83 MxTip = 0
84 MyTip = 0
85
86 ks1=1 #shear correction, torsion
87 ks2=1 #shear correction, bending
88 ks3=1 #shear correction, bending
89 ff=1 #drawing factor
90
91 if case == 0 or case == 1:
92 caseName = cases[case]
93
94 L = 2 #length of beam
95 w = 0.1 #width of beam
96 h = 0.5 #height Y
97
98 fTip = 5e5*h**3
99 if case == 1:
100 fTip *= 1000
101
102 Em = 2.07e11
103 rho = 1e2
104
105 A=h*w
106 nu = 0.3 # Poisson ratio
107 ks2= 10*(1+nu)/(12+11*nu)
108 ks3=ks2
109
110 elif case == 2:
111 L = 2 #length of beam
112 h = 0.2 #height Y
113 w = 0.4 #width Z of beam
114 Em = 2.07e11
115 rho = 1e2
116
117 A=h*w
118
119 nu = 0.3 # Poisson ratio
120 ks1= 0.5768 #torsion correction factor if J=Jyy+Jzz
121 ks2= 0.8331
122 ks3= 0.7961
123
124 MxTip = 0.5e6
125 MyTip = 2e6
126
127 csFact = 10
128 elif case == 3 or case == 4: #Princeton beam example
129 L = 0.508 #length of beam
130 h = 12.3777e-3 #height Y; 12.3777e-3 with Obrezkov's paper
131 w = 3.2024e-3 #width Z of beam
132 Em = 71.7e9
133 ks1=0.198
134 nu = 0.31
135
136 ks2=1
137 ks3=1
138 # ks2=0.9
139 # ks3=0.9
140
141
142 rho = 1e2 #unused
143 A=h*w
144
145 MxTip = 0
146 MyTip = 0
147 if case == 3:
148 fTip = 8.896 #F2
149 elif case == 4:
150 fTip = 13.345 #F3
151 #if kk==0: exu.Print('load=', fTip)
152
153 beta = betaDegree/180*pi #beta=0 => negative y-axis
154 bodyFixedLoad = False
155
156 csFact = 10
157
158 Gm = Em/(2*(1+nu)) # Shear modulus
159
160 # Cross-section properties
161 Iyy = h*w**3/12 # Second moment of area of the beam cross-section
162 Izz = w*h**3/12 # Second moment of area of the beam cross-section
163 J = (Iyy+Izz) # approximation; Polar moment of area of the beam cross-section
164
165 sectionData.stiffnessMatrix = np.diag([Em*A, Gm*A*ks2, Gm*A*ks3, Gm*J*ks1, Em*Iyy, Em*Izz])
166
167
168 rhoA = rho*A
169
170 if False:
171 #linear solution:
172 uzTip = fTip*L**3/(3*Em*Iyy)
173 exu.Print('uz linear=',uzTip)
174 uyTip = fTip*L**3/(3*Em*Izz)
175 exu.Print('uy linear=',uyTip)
176
177 sectionData.inertia= rho*J*np.eye(3)
178 sectionData.massPerLength = rhoA
179
180 sectionGeometry = exu.BeamSectionGeometry()
181
182 #points, in positive rotation sense viewing in x-direction, points in [Y,Z]-plane
183 #points do not need to be closed!
184 lp = exu.Vector2DList()
185 if True:
186 lp.Append([h*ff,-w*ff])
187 lp.Append([h*ff,w*ff])
188 lp.Append([-h*ff,w*ff])
189 lp.Append([-h*ff,-w*ff])
190
191 sectionGeometry.polygonalPoints = lp
192 #exu.Print('HERE\n',sectionGeometry.polygonalPoints)
193 nGround = mbs.AddNode(NodePointGround(referenceCoordinates=[0,0,0])) #ground node for coordinate constraint
194 mnGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nGround, coordinate=0))
195
196
197 eY=[0,1,0]
198 eZ=[0,0,1]
199 lElem = L/nElements
200 useGeometricallyExact = True
201 if compute3D:
202 if useGeometricallyExact:
203 n0 = mbs.AddNode(NodeClass(referenceCoordinates=[0,0,0]+initialRotationsGE))
204 else:
205 initialRotations = eY+eZ
206 n0 = mbs.AddNode(NodePointSlope23(referenceCoordinates=[0,0,0]+initialRotations))
207 nInit = n0
208 for k in range(nElements):
209 if useGeometricallyExact:
210 n1 = mbs.AddNode(NodeClass(referenceCoordinates=[lElem*(k+1),0,0]+initialRotationsGE))
211
212 oBeam = mbs.AddObject(ObjectBeamGeometricallyExact(nodeNumbers=[n0,n1], physicsLength = lElem,
213 sectionData = sectionData,
214 visualization=VBeam3D(sectionGeometry=sectionGeometry)))
215 else:
216 n1 = mbs.AddNode(NodePointSlope23(referenceCoordinates=[lElem*(k+1),0,0]+initialRotations))
217 oBeam = mbs.AddObject(ObjectANCFBeam(nodeNumbers=[n0,n1], physicsLength = lElem,
218 #testBeamRectangularSize = [h,w],
219 sectionData = sectionData,
220 crossSectionPenaltyFactor = [csFact,csFact,csFact],
221 visualization=VANCFBeam(sectionGeometry=sectionGeometry)))
222 n0 = n1
223
224
225 mTip = mbs.AddMarker(MarkerNodeRigid(nodeNumber = n1))
226 if fTip != 0:
227 if case < 3:
228 mbs.AddLoad(Force(markerNumber=mTip, loadVector = [0,fTip,0], bodyFixed = bodyFixedLoad))
229 elif case >= 3:
230 mbs.AddLoad(Force(markerNumber=mTip, loadVector = [0,-fTip*cos(beta),fTip*sin(beta)], bodyFixed = bodyFixedLoad))
231
232 if MxTip != 0 or MyTip != 0:
233 mbs.AddLoad(Torque(markerNumber=mTip, loadVector = [MxTip, MyTip,0]))#, bodyFixed = True ))
234
235 if useGeometricallyExact:
236 nm0 = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nInit))
237 nmGround = mbs.AddMarker(MarkerNodeRigid(nodeNumber=nGround))
238 mbs.AddObject(GenericJoint(markerNumbers=[nmGround, nm0]))
239 else:
240 for i in range(9):
241 #if i != 4 and i != 8: #exclude constraining the slope lengths
242 if True:
243 nm0 = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nInit, coordinate=i))
244 mbs.AddObject(CoordinateConstraint(markerNumbers=[mnGround, nm0]))
245
246
247 # exu.Print(mbs)
248 mbs.Assemble()
249
250 tEnd = 100 #end time of simulation
251 stepSize = 0.5*0.01*0.1 #step size; leads to 1000 steps
252
253 simulationSettings = exu.SimulationSettings()
254 simulationSettings.solutionSettings.solutionWritePeriod = 2e-2 #output interval general
255 simulationSettings.solutionSettings.sensorsWritePeriod = 1e-1 #output interval of sensors
256 simulationSettings.timeIntegration.numberOfSteps = int(tEnd/stepSize) #must be integer
257 simulationSettings.timeIntegration.endTime = tEnd
258 #simulationSettings.solutionSettings.solutionInformation = "This is the info\nNew line\n and another new line \n"
259 simulationSettings.timeIntegration.generalizedAlpha.spectralRadius = 0.5
260 #simulationSettings.timeIntegration.simulateInRealtime=True
261 #simulationSettings.timeIntegration.realtimeFactor=0.1
262
263 simulationSettings.timeIntegration.verboseMode = verbose
264 simulationSettings.staticSolver.verboseMode = verbose
265
266 #simulationSettings.parallel.numberOfThreads = 4
267 simulationSettings.timeIntegration.newton.useModifiedNewton = True
268 #simulationSettings.timeIntegration.newton.numericalDifferentiation.minimumCoordinateSize = 1e0
269
270 #simulationSettings.timeIntegration.newton.numericalDifferentiation.relativeEpsilon = 1e-4
271 simulationSettings.timeIntegration.newton.relativeTolerance = 1e-6
272
273 # simulationSettings.displayComputationTime = True
274 simulationSettings.linearSolverType = exu.LinearSolverType.EigenSparse
275 # simulationSettings.parallel.numberOfThreads = 4
276
277 #simulationSettings.staticSolver.newton.numericalDifferentiation.relativeEpsilon = 5e-5
278 #simulationSettings.staticSolver.newton.numericalDifferentiation.forODE2 = True
279 #simulationSettings.staticSolver.newton.relativeTolerance = 1e-6
280 # simulationSettings.staticSolver.newton.numericalDifferentiation.relativeEpsilon = 1e-4
281
282 simulationSettings.staticSolver.numberOfLoadSteps = 5
283 simulationSettings.staticSolver.adaptiveStep = True
284 #simulationSettings.staticSolver.stabilizerODE2term = 100
285
286 if useGeometricallyExact:
287 # simulationSettings.staticSolver.newton.numericalDifferentiation.forODE2 = True
288 # simulationSettings.staticSolver.newton.numericalDifferentiation.relativeEpsilon = 1e-5
289 # simulationSettings.staticSolver.newton.relativeTolerance = 1e-5
290 # simulationSettings.staticSolver.newton.absoluteTolerance = 1e-5
291 if case == 0:
292 simulationSettings.staticSolver.newton.relativeTolerance = 1e-4
293 simulationSettings.staticSolver.newton.absoluteTolerance = 1e-5
294 simulationSettings.staticSolver.numberOfLoadSteps = 1 #otherwise makes problems
295
296 if nElements > 32 and case==0: #change tolerance, because otherwise no convergence
297 simulationSettings.staticSolver.newton.relativeTolerance = 1e-6
298 if case == 1: #tolerance changed from 1e-8 to 5e-10 to achieve values of paper (1024 has difference at last digit in paper)
299 simulationSettings.staticSolver.newton.relativeTolerance = 0.5e-9
300
301
302
303 #add some drawing parameters for this example
304 SC.visualizationSettings.nodes.drawNodesAsPoint=False
305 SC.visualizationSettings.nodes.defaultSize=0.01
306
307 SC.visualizationSettings.bodies.beams.axialTiling = 50
308 SC.visualizationSettings.general.drawWorldBasis = True
309 SC.visualizationSettings.general.worldBasisSize = 0.1
310 SC.visualizationSettings.openGL.multiSampling = 4
311
312
313 # [M, K, D] = exu.solver.ComputeLinearizedSystem(mbs, simulationSettings, useSparseSolver=True)
314 # exu.Print('M=',M.round(1))
315
316 if useGraphics:
317 exu.StartRenderer()
318 mbs.WaitForUserToContinue()
319
320 # if computeEigenmodes:
321 # nModes = 3*(1+int(compute3D))
322 # nRigidModes = 3*(1+int(compute3D))
323 # if compute2D:
324 # constrainedCoordinates=[0,1,mbs.systemData.ODE2Size()-2]
325 # else:
326 # constrainedCoordinates=[0,1,2,5,mbs.systemData.ODE2Size()-8,mbs.systemData.ODE2Size()-7]
327
328 # # constrainedCoordinates=[]
329
330 # compeig=mbs.ComputeODE2Eigenvalues(simulationSettings, useSparseSolver=False,
331 # numberOfEigenvalues= nRigidModes+nModes,
332 # constrainedCoordinates=constrainedCoordinates,
333 # convert2Frequencies= False)
334
335 # exu.Print('eigvalues=',np.sqrt(compeig[0][nRigidModes:]))
336
337 # if False: #show modes:
338 # for i in range(nModes):
339 # iMode = nRigidModes+i
340 # mbs.systemData.SetODE2Coordinates(5*compeig[1][:,iMode], exudyn.ConfigurationType.Visualization)
341 # mbs.systemData.SetTime(np.sqrt(compeig[0][iMode]), exudyn.ConfigurationType.Visualization)
342 # mbs.SendRedrawSignal()
343
344 # mbs.WaitForUserToContinue()
345
346 # else:
347 mbs.SolveStatic(simulationSettings)
348 # mbs.SolveDynamic(simulationSettings)
349 #mbs.SolveDynamic(simulationSettings, solverType = exu.DynamicSolverType.RK44)
350
351 #check jacobian
352 if False:
353 #%%+++++++++++++++++++++++++++++++++++
354 solver=mbs.sys['staticSolver']
355 solver.InitializeSolver(mbs, simulationSettings)
356 solver.ComputeJacobianODE2RHS(mbs)
357 J=solver.GetSystemJacobian()
358 print((1e-6*J[:14,:7]).round(3))
359 print((1e-6*J[:14,7:14]).round(3))
360
361 #%%+++++++++++++++++++++++++++++++++++
362 if useGraphics:
363 SC.WaitForRenderEngineStopFlag()
364 exu.StopRenderer() #safely close rendering window!
365
366 ##evaluate final (=current) output values
367 uTip = mbs.GetNodeOutput(n1, exu.OutputVariableType.Displacement)
368
369 errorFact = 1
370 if case != 1:
371 errorFact *= 100
372
373 testErrorSum += np.linalg.norm(uTip)
374
375
376
377 if case < 2:
378 pTip = mbs.GetNodeOutput(n1, exu.OutputVariableType.Position)
379 exu.Print('ne=',nElements, ', ux=',L-pTip[0], ', uy=',pTip[1])
380 elif case == 2:
381 rotTip = mbs.GetNodeOutput(n1, exu.OutputVariableType.Rotation)
382 exu.Print('ne=',nElements, ', u=',list(uTip))
383 # exu.Print('ne=',nElements, ', rot=',rotTip)
384 elif case == 3 or case == 4:
385 exu.Print('ne=', nElements, ', beta=', round(beta*180/pi,1), ', u=',uTip.round(7))
386
387
388exu.Print('Solution of geometricallyExactBeamTest=', testErrorSum)
389exudynTestGlobals.testError = testErrorSum - (1.012822053539261)
390exudynTestGlobals.testResult = testErrorSum
391
392
393#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
394#all results are taken from ANCFBeam (shear deformable 2-node 3D beam):
395#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
396
397# case= 0/CantileverLinear2011
398#NachbagauerPechsteinIrschikGerstmayrMUBO2011 (2D):
399# ne=1, 9.12273046e–8, 6.16666566e–4, 0.000193
400# ne=2, 1.61293091e–7, 7.61594059e–4, 4.831e–5
401# ne=4, 1.81763233e–7, 7.97825954e–4, 1.208e–5
402# ne=256, 1.88847418e–7, 8.09900305e–4, 2.945e–9
403#Exudyn: ksFact=1
404# ne= 1 , ux= 9.122730637578513e-08 , uy= 0.0006166665660910789
405# ne= 2 , ux= 1.612930911054633e-07 , uy= 0.0007615940599560586
406# ne= 4 , ux= 1.8176323512975046e-07 , uy= 0.0007978259537503566
407# ne= 8 , ux= 1.8706537496804287e-07 , uy= 0.0008068839288072378
408# ne= 16 , ux= 1.8840244964124508e-07 , uy= 0.0008091484226773518
409# ne= 32 , ux= 1.887374359021976e-07 , uy= 0.0008097145461515286
410# ne= 64 , ux= 1.888212299849812e-07 , uy= 0.0008098560770202866
411# ne= 128 , ux= 1.8884218011550047e-07 , uy= 0.000809891459736643
412# ne= 256 , ux= 1.8884741770364144e-07 , uy= 0.0008099003054122335
413
414
415# case= 1/Cantilever2011
416#NachbagauerPechsteinIrschikGerstmayrMUBO2011 (2D):
417# ne=1, 0.07140274, 0.54225823, 0.168310
418# ne=2, 0.12379212, 0.65687111, 0.053697
419# ne=4, 0.14346767, 0.69593561, 0.014633
420# ne=1024, 0.15097103, 0.71056837, 2.280e–7
421
422#Exudyn: ksFact=1
423# ne= 1 , ux= 0.07140273975041422 , uy= 0.5422582285449739
424# ne= 2 , ux= 0.12379212054619537 , uy= 0.6568711099777776
425# ne= 4 , ux= 0.14346766617229956 , uy= 0.695935613449867
426# ne= 8 , ux= 0.14904162148449163 , uy= 0.7068152604035266
427# ne= 16 , ux= 0.15048521526298897 , uy= 0.709623891842095
428# ne= 32 , ux= 0.15084943688011565 , uy= 0.7103320154655514
429# ne= 64 , ux= 0.15094070328691145 , uy= 0.7105094267817303
430# ne= 128 , ux= 0.15096353326024237 , uy= 0.7105538037895819
431# ne= 256 , ux= 0.15096924149743085 , uy= 0.7105648993600513
432# ne= 512 , ux= 0.15097066651939461 , uy= 0.7105676689547459
433# ne= 1024 , ux= 0.15097102364723924 , uy= 0.7105683631862169
434
435# case = 2:
436#2013 CND, Nachbagauer Gruber Gerstmayr (static, 3D); "Structural and Continuum Mechanics Approaches for a 3D Shear Deformable ANCF Beam Finite Element: Application to Static and Linearized Dynamic Examples"
437#Table 4:
438# SMF
439# 8, 1.0943e-4, 1.8638e-4, 1.8117e-2
440# 32, 1.0943e-4, 1.8625e-4, 1.8117e-2
441# ANSYS
442# 40, 1.0939e-4, 1.8646e-4, 1.8117e-2
443#Exudyn, ksFact=10:
444# ne= 8 , u= [-0.00010900977088157404, -0.0001902100873246334, -0.01811732779800177]
445# ne= 32 , u= [-0.00010941122286522997, -0.00018667435478355072, -0.01811739809277171]
446# ne= 128 , u= [-0.00010943631319815239, -0.000186451835025629, -0.018117402461210096]
447#==> in 2013 paper, element performed slightly better, especially in ux and uy terms
448
449# case = 3:
450#Princeton beam with ANSYS (Leonid Obrezkov / Aki Mikkola / Marko Matikainen et al.,
451# Performance review of locking alleviation methods for continuum ANCF beam elements,
452# Nonlinear Dynamics, Vol. 109, pp. 31–546, May 2022
453# beta=[0 15 30 45 60 75 90];
454if (case==3 or case == 4) and False:
455 # F2=8.896
456 # % ANSYS beam (10-199 el)
457 ANSYSF2y=np.array([1.071417630E-002, 1.061328706E-002, 1.011169630E-002, 8.837226265E-003, 6.604665004E-003, 3.538889001E-003, 0])
458 ANSYSF2z=np.array([0, 4.208232124E-002, 7.939482948E-002, 0.108987937, 0.129887616, 0.142194370, 0.146245978])
459 exu.Print('refsol ANSYS F2=8.896:\n',ANSYSF2y.round(6), '\n', ANSYSF2z.round(6))
460 # % ANSYS solid (el) (4x12x500) - finer mesh doesn't have much influence see in Size effect file
461 # ANSYS_solid_y=[1.069752828E-002 1.057180106E-002 9.938278402E-003 8.686786771E-003 6.500006282E-003 3.481999513E-003 0];
462 # ANSYS_solid_z=[0 4.101165651E-002 7.696749069E-002 0.105976311 0.127251299 0.139594740 0.143848652];
463
464 # F3=13.345
465 # % ANSYS beam (10-199 el)
466 ANSYSF3y=np.array([1.606423724E-002, 1.645825752E-002, 1.665873206E-002, 1.518618440E-002, 1.157837500E-002, 6.248967384E-003, 0])
467 ANSYSF3z=np.array([0, 6.435812858E-002, 0.117735994, 0.156467239, 0.181861627, 0.196097131, 0.200677707])
468 # % ANSYS solid (el) (4x12x500) - finer mesh see in Size effect file
469 #ANSYS_solid_y=[1.603700622E-002 1.637026068E-002 1.640440775E-002 1.485055210E-002 1.127173264E-002 6.062461977E-003 0])
470 #ANSYS_solid_z=[0 6.270699533E-002 0.113752002 0.153554457 0.179978534 0.192972233 0.197669499])
471 exu.Print('refsol ANSYS F3=13.345:\n',ANSYSF3y.round(6), '\n', ANSYSF3z.round(6))
472#Exudyn results for Princeton beam:
473#not exactly the same, but around the previous values with HOTINT
474#using 16 elements, csFact=10 (no influence)
475# F2=8.896
476# case= 3, PrincetonBeam
477# ne= 16 , beta= 0.0 , u= [-0.0001352 -0.0107023 0. ]
478# ne= 16 , beta= 15.0 , u= [-0.0022414 -0.0106295 0.0421374]
479# ne= 16 , beta= 30.0 , u= [-0.0076567 -0.0101861 0.0794434]
480# ne= 16 , beta= 45.0 , u= [-0.0143664 -0.0089529 0.1089703]
481# ne= 16 , beta= 60.0 , u= [-0.0204225 -0.0067182 0.1297877]
482# ne= 16 , beta= 75.0 , u= [-0.0245093 -0.0036079 0.1420319]
483# ne= 16 , beta= 90.0 , u= [-0.0259403 -0. 0.1460608]
484
485# F3=13.345
486# case= 4, PrincetonBeam
487# ne= 16 , beta= 0.0 , u= [-0.0003039 -0.0160454 0. ]
488# ne= 16 , beta= 15.0 , u= [-0.005319 -0.0165469 0.064622 ]
489# ne= 16 , beta= 30.0 , u= [-0.0171901 -0.0169316 0.1179818]
490# ne= 16 , beta= 45.0 , u= [-0.0303357 -0.0155488 0.1565214]
491# ne= 16 , beta= 60.0 , u= [-0.0411035 -0.0118996 0.1817173]
492# ne= 16 , beta= 75.0 , u= [-0.0479101 -0.0064334 0.1958343]
493# ne= 16 , beta= 90.0 , u= [-0.0502184 -0. 0.2003738]