search

mo-hanxuan / FEMcy

a finite element solver based on Taichi, being parallel (CPU/GPU), portable and open-source
MIT License
127 stars 16 forks source link
finite-element gpu python taichi

readme

FEMcy

a finite element solver based on Taichi, being parallel (CPU/GPU), portable and open-source

FEMcy is a finite element solver for structural analysis in continuum mechanics, powered by cross-platform parallel (CPU/GPU) computing language of Taichi. FEMcy is flexible for customized needs by open-source. The mechanism behind computational structural analysis (CSD) is now opened for you, and can be manipulated by you to meet your customized needs. We provide the implementation on GPU parallel computing, meanwhile maintain the friendly readability by Python language.

Features

currently supported Elements

currently supported Boundary Conditions

Installation and Usage

  1. install Python (3.8+) and pip, then install Taichi, numpy and scipy

    pip install taichi
    pip install numpy
    pip install scipy

  2. git clone this project
  3. go to the current directory of this project, and run the code by:

    python ./main.py

  4. Pre-processing: choose an .inp file (the Abaqus-input-file-format) which defines the geometry, mesh, material and boundary condition for structual analysis. The .inp file can be obtained by Abaqus pre-processing GUI. For example, insert the path and inp file name to the command line:

    please give the .inp format's input file path and name: tests/beam_deflection/load800_freeEnd_largeDef/beamDeflec_quadPSE_largeD_load800.inp

    more examples of inp files can be found at ./tests folder

  5. after convergence, the deformed body colored by mises-stress (defaultly) is showed at the window.

    Examples (some benchmark problems)

    1. elliptic-membrane with constant pressure

Fig. 1 shows the geometric and consititutive model definition of elliptic membrane problem, which can be refered to CoFEA benchmark. The stress $\sigma_{yy}$ at point D is expected to be 92.7 MPa.

$u{x}$ = 0 for edge DC, $u{y}$ = 0 for edge AB, uniform normal pressure of 10 MPa on edge BC

property value unit
Young's modulus, E $2.1\times10^{5}$ MPa
Poisson's ratio, $\nu$ 0.3 -

Results are showed below:

image-20220909163740116
Fig. 1 Results of elliptic membrane under normal pressure. (a) geometric model definition; (b) results from CoFEA; (c, d) results of linear triangular element from Abaqus and FEMcy respectively; (e, f) results of quadratic triangular element from Abaqus and FEMcy respectively, shown by nodal stress extrapolated from integration points.

Results of simulation:

$\sigma_{yy}$ [MPa] at point D linear element quadratic element
$\sigma_{yy}$ at node (extrapolated) $\sigma_{yy}$ at integration point
Abaqus 93.45 93.34 84.42
FEMcy 93.56 93.32 84.40
relative error 0.12% 0.021% 0.024%

The relative error shown above is the error between results of Abaqus and FEMcy, indicating that the results of these two softwares are almost the same.

Another interesting point is that, the linear element needs a very locally-refined mesh to get close to the expected result (~ 92.7 MPa), whereas the quadratic triangular element can get the accurate result at a relative coarse mesh as shown in the above figure.

corresponding .inp file:

./tests/elliptic_membrane/element_linear/ellip_membrane_linEle_localVeryFine.inp

./tests/elliptic_membrane/element_quadratic/ellip_membrane_quadritic_trig_neumann.inp

results can be compared with https://cofea.readthedocs.io/en/latest/benchmarks/004-eliptic-membrane/model.html

2. Beam Deflection (comparison of small and large deformation)

For a cantilever beam as shown in Fig. 2, when imposing surface traction at the free end, the results show great differences between small and large deformation.

Small deformation (linear elastic) has no horizontal displacement at the free end, whereas large deformation (geometric nonlinear) shows a bending beam with horizontal displacement at the free end. Curves of vertical deflection vs. load are shown in Fig. 2 (d), which shows that the geometric nonlinear case has nonlinear deflection-load relation, and results of FEMcy are pretty close to Abaqus.

Results of large deformation (geometric nonlinear) is computed by Newton's method with adaptive time step and relaxation factor, where the Jacobian matrix $\partial\ nodal\ force\ /\ \partial\ nodal\ displacement$ at each Newton-step is approximated by the stiffness matrix [K].

Fig. 2 (a) Cantilever Beam, with surface traction at the free end. (b) Linear elastic case has no horizontal displacement at the free end, whereas (c) geometric nonlinear case shows a bending beam with horizontal displacement at the free end. (d) Curves of vertical deflection vs. load

(The material here has 2.e5 MPa of Young's modulus and 0.3 of Poisson's ratio)

3. Twist Plate

Run Python main.py and then insert the .inp file tests/twist/twist_plate_C3D4.inp or tests/twist/twist_plate_C3D10.inp to the command line.

The results are shown in Fig. 3 with both linear and quadratic tetrahedral elements. You can see that at the cross section perpendicular to rotation axis, stress is higher at the outer edge with larger radius at the cross section, which agrees well with the classical torsion problem.

C3D4 linear Tetrahedral C3D10 quadratic tetrahedral
Fig. 3 Twist plate with C3D4 and C3D10 element respectively, colored by Mises stress.

Future work

References