Samples of various research activities in Computational Engineering
The tetrahedral finite cell method:Higher-order immersogeometric analysis on adaptive non-boundary-fitted meshes
in cooperation with
University of Minnesota (Dr. Vasco Varduhn, Dr. D. Schillinger)
Iowa State University (Dr. Ming-Chen Hsu)
The finite cell method (FCM) is an immersed domain finite element method that combines higher-order non-boundary-fitted meshes, weak enforcement of Dirichlet boundary conditions, and adaptive quadrature based on recursive subdivision. Here, we extend the FCM, so far only used with Cartesian hexahedral elements, to higher-order non-boundary-fitted tetrahedral meshes, based on a reformulation of the octree-based subdivision algorithm for tetrahedral elements. We show that the resulting TetFCM scheme is fully accurate in an immersogeometric sense, that is, the solution fields achieve optimal and exponential rates of convergence for h-refinement and p-refinement, if the immersed geometry is resolved with sufficient accuracy (Int. J. Numer. Meth. Engng (2016))
Figure: Thick plate benchmark with fictitious domain modeled cut-out. Symmtric part of a plate with circular hole under tension. Discretization and local refinement with high-order tetrahedral elements.
Figure: Analysis results: von Mises stress distribution around the hole and convergence plot for various tetrahedral elements and approximation spaces.
Figure: Circumferential stress around the hole for different levels of geometric resolution through means of elementwise refinement.
Figure: Image based bone geometry and tetrahedral discretization model with refined phase-field approximation of a crack. The femur is loaded on the femural head with a compression load and fully clamped at the distal face.
Figure: Deformed, cracked femur bone and principal stress concentrations.
Multi-level hp-adaptivity for cohesive fracture modeling
in cooperation with
TU München (Nils Zander et al.)
Discretization-induced oscillations in the load–displacement curve are a well-known problem for simulations of cohesive crack growth with finite elements. The problem results from an insufficient resolution of the complex stress state within the cohesive zone ahead of the crack tip. Here, the proposed method demonstrates that the hp-version of the finite element method is ideally suited to resolve this complex and localized solution characteristic with high accuracy and low computational effort. To this end, we use a local and hierarchic mesh refinement scheme that follows dynamically the propagating crack tip. Furthemore, the combination of the approach with the finite cell method allows for a simplified modeling of cut-outs, holes and other geometric features. (Int. J. Numer. Meth. Engng (2016)) .
Figure: Perforated double cantilever beam benchmark problem modeled with ficticious perforation cut-outs.
Figure: dynamically adapted hierarchical and local refinement along the crack tip at various time steps. The total numer of degrees-of-freedom changes with the crack-tip-adapted refinement.
Figure: Load-displacement curve for the perforated double cantilever beam using different discretizations. The local hierarchic refinement reduces the problem by a factor five in terms of degrees-of-freedom.
Weak coupling of thin shell models and blended shell models, weak enforcement of essential boundary conditions
in cooperation with
Nanjing University of Aeronautics and Astronautics (Dr. Yujie Guo)
A continuity preserving Nitsche-based coupling approach is introcuded for thin isogeometric Kirchhoff-Love shell structures and blended shell structures including the coupling of thin and solid NURBS shells. To model complex shell structures which need to be assembled from multiple NURBS patches, the bending stiffness should be maintained across the patch interfaces. We propose a variationally consistent weak coupling method for thin-walled shell patches. The method overcomes the need for C1-continuity along the patch interface to ensure a corresponding geometric continuity in the deformed configuration and a correct transfer of bending moments across the interface. The proposed approach retains the high level of accuracy of single patch solutions and reveals its potential for authentic multi-patch NURBS modeling. The enforcement of boundary conditions of thin shells and trimmed shell patches follows in accordance to the coupling approach considering additional external flux terms in the formulation.
Figure: (top) Model and NURBS patch structure of a pinched cylinder model, modeled as pure thin shell (KL shell theory) and as blended shell (KL and solid shell theory). (bottom) Isogeometric analysis results of the shell model. A comparison of the deformed structure shows virtually no difference in the analysis results maintaining the high solution quality of higher order NURBS isogeometric analysis.
Figure: (top) NURBS T-Beam models, conforming and non-conforming patch structure. (bottom) Isogeometric analysis results (corner single force loaded) of a weak coupling approach considering all shear terms to provide variational consistency and optimal convergence behavior. The approach is stabilized with a penalty-like stability term. The free stability parameter is either chosen ad-hoc or determined on the basis of suited bounds satisfying the governing inverse inequatlity condition of the problem (Comput. Methods Appl. Mech. Engrg. 284 (2015) 881–905).
Figure: (top) Scordelis-Lo benchmark test (self-weight loading) with weakly enforced boundary conditions using symmetriy properties. (bottom) Convergence plot for uniform refinement which shows an almost optimal behavior for an ad-hoc choice of the stability parameter. Improvement of the results can be expected for a refined stability paramter analysis as proposed [IJNME Ruess et al. 2013, CMAME Ruess et al. 2014].
Figure: (top) Weakly enforced boundary conditions of a trimmed circular NURBS plate model with uniform loading, trimming according to the principles of the finite cell method. (bottom) Convergence plot for a uniform patch refinement. The stability parameter was chosen empirically.
Figure: (top) Trimmed NURBS multi-patch structure with weak enforcement of coupling and essential boundary conditions. The structure is loaded with an inner pressure load. (bottom) Analysis results: displacements of the deformed structure and von Mises stress distribution, both showing a highest degree of smoothness and continuity across the coupling interfaces and boundaries.
Reduced order modeling in the framework of the Koiter-Newton method
in cooperation with
China Academy of Space Technology (Dr. Ke Liang)
Delft University of Technology (Dr. Mostafa Abdalla)
The Koiter-Newton method is a reduced order modeling technique which allows us to trace efficiently the entire equilibrium path of a non-linear structural analysis. In the framework of buckling the method is capable to handle snap-back and snap-through phenomena and in an extended version to predict reliably bifurcation branches along the equilibrium path.
Figure: Basic principle of the Koiter Newton approach in comparison to classical NR steps. Non-linear predictor steps significantly reduce the computational complexity as proposed by Liang (PhD-thesis K. Liang, TU Delft, 2013, Comput. Methods Appl. Mech. Engrg. 279 (2014) 440–468).
Figure: Flow chart of the extended reduced-order model Koiter Newton approach using a polynomial homotopy continuation method to improve the post-buckling predition.
Figure: Model and response curves about the lateral displacement at the loading point for a cylindrical roof stability analysis. The considered model variants include an isotropic model and two laminate composite models. The Koiter-Newton results are completely coincident with an ABAQUS derived reference solution. The single predictor-corrector steps of the KN-approach are indicated by triangle markers. Since the first buckling modes of this model are well-sparated a reduced order model of dimension two was sufficient. Each of the response curves required 89 solution steps of the full FE system of equations with ABAQUS and only 9 (!) corresponding solution steps with the KN-method (Finite Elements in Analysis and Design 116 (2016) 38–54).
Figure: Response curve for a cylinder shell roof with snap-back, snap-through behavior. (left) Original Koiter-Newton approach solved with 7 KN-steps (including 26 equilibrium iterations) compared to 63 solution steps with ABAQUS (including 123 equilibrium iterations). (right) Extended KN approach for post buckling analysis revealing primary and secondary buckling paths, solved in 6 KN-steps and 22 iterations compared to 116 iterations of a classical branch-switching approach. The KN method constructed and used a reduced order model of dimension four(International Journal of Non-Linear Mechanics 81 (2016) 129–138).
Figure: First buckling mode of a cylindrical stringer-stiffened shell structure under consideration of a single-perturbation load modeling a geometric imperfection. The stringer shell model was considered in terms of a smeared material model. All models, an ABAQUS reference model and two KN-based models, revealed 11 full waves in circumferential direction and a single full wave in axial direction for the first mode. The Koiter-Newton results were found with a reduced order model of dimension 20, each degree-of-freedom representing closely spaced modes (Aerospace Science and Technology 55 (2016) 103–110).
Laminate composites in the framework of isogeometric analysis
in cooperation with
Nanjing University of Aeronautics and Astronautics (Dr. Yujie Guo)
Layerwise theories provide an accurate prediction of the three-dimensional stress state of composite structures which is in many cases a prerequisite for a successful prediction of failure mechanism. We focus on such a model for laminate composite plate and shell models in the framework of isogeometric analysis. We apply an exact but dimensionally reduced geometry description based on a computer-aided design (CAD) derived mid-surface NURBS model and a layerwise through-thickness B-spline description. Using the unique and characteristic refinement capabilities of isogeometric analysis we simplify the modeling process while keeping the higher order and higher continuity properties of NURBS. Trimmed geometries, as common in CAD models are efficiently handled by an embedded domain method which treats trimmed areas as a fictitious extension domain in the computational model.
Figure: (top) Layerwise model of a simply supported square plate subjected to a parabolic loading. (bottom) Through-thickness stresses comparison, [0° - 90 ° - 0°] square plate: (a) in-plane normal stress S11, (b) transverse normal stress S33, (c) transverse shear stress S13, (d) transverse shear stress S23. All results fully coincide with the exact 3D solution.
Figure: Cylindrical layerwise shell model with internal sinusoidal pressure load.
Figure: Through thickness transverse stresses comparison of [90° - 0° - 90°] cylindrical shell: (a) transverse shear stress S13, (b) transverse shear stress S23, (c) transverse normal stress S33.
Samples of earlier research activities (-- 2013) :
Weak enforcement of essential boundary conditions in the framework of the p-version FEM and isogeometric analysis combined with the finite cell method
in cooperation with
ICES, University of Texas, Austin (Dr. Schillinger)
Jacobs School of Engineering, UCSD (Prof. Bazilevs)
Chair for Computation in Engineering, TUM (Prof. Rank)
Enforcing essential boundary conditions plays a central role in immersed boundary methods. Nitsche’s idea has proven to be a reliable concept to satisfy weakly boundary and interface constraints. We use an extension of Nitsche’s method for elasticity problems in the framework of higher order and higher continuity approximation schemes such as the B-spline and non-uniform rational basis spline version of the finite cell method or isogeometric analysis on trimmed geometries. Furthermore, we illustrate a significant improvement of the flexibility and applicability of this extension in the modeling process of complex 3D geometries.
Figure: modeling sequence of a trimmed NURBS geometry: (a) single patch NURBS model (b) trimming curve definition (c) finally trimmed visual model = fictitious domain suitable model.
Figure: trimmed shell model under self weight: partially clamped (weak BC) at the lower rim (a) deflection model (b) von Mises stress distribution with element structure
Figure: model of a Laplace problem with sine load along the lower edge and weakly enforced homogeneous boundary conditions along the remaining edges: embedded finite cell model
Figure: (a) finite cell solution with embedding domain (b) absolute error (log) with embedding domain.
Figure: surface model of a crankshaft (left) and voxelized model with 62.5 mio voxel.
Figure: displacement field: pressure loads at all 4 crank pins, clamped (weakly enforced BC) at the crank journals.
Weak enforcement of coupling constraints for multi-patch NURBS model in the framework of a higher order fictitious domain method
in cooperation with
ICES, University of Texas, Austin (Dr. Schillinger)
Chair for Computation in Engineering, TUM (Prof. Rank)
Nitsche’s method can be used as a coupling tool for non-matching discretizations by weakly enforcing interface constraints. We explore the use of weak coupling based on Nitsche’s method in the context of higher order and higher continuity B-splines and NURBS. We demonstrate that weakly coupled spline discretizations do not compromise the accuracy of isogeometric analysis. It is further shown that the combination of weak coupling with the finite cell method opens the door for a truly isogeometric treatment of trimmed B-spline and NURBS geometries that eliminates the need for costly re-parameterization procedures.
Figure: metall hook model consisting of NURBS patches and finite cell patches for a simplified modeling, all patches are weakly coupled.
Figure: displacement field of the deformed hook.
Figure: von Mises stress distribution on the deformed model.
Figure: plate with a hole under tension: model data
Figure: convergence of various coupling model for uniform p-refinement including patch refinement around the hole (left), von Mises stress distribution (right)
Numerical prediction of the elasticity behavior of patient-specific femur bones
in cooperation with
Dept. of Mechanical Engineering, Ben-Gurion University of the Negev, Beer Sheva, Israel (Prof. Yosibash)
Standard methods for predicting bone’s mechanical response from quantitative computer tomography (qCT) scans aremainly based on classical h-version finite element methods (FEMs). Due to the low-order polynomial approximation, the need for segmentation and the simplified approach to assign a constant material property to each element in h-FE models, these often compromise the accuracy and efficiency of h-FE solutions. Here, a non-standard method, the finite cell method (FCM), is proposed for predicting the mechanical response of the human femur.
Figure: Model & results of a human femur elasticity analysis: (a) embedding knot-span cell model (b) displacement field of the bottom clamped and top loaded femur (1000N) (c) von Mises equivalent strain distribution.
Figure: Femur validation results: correlation between in-vitro experiment and numerical prediction with p-version FCM and B-spline version FCM: surface measured micro-strains and vertical & horizontal deflection of the femur head.
Computational Steering for Orthopeadics
Int'l Graduate School of Science & Engineering, Technische Universität München.
Within this interdisciplinary research project a computational steering platform was developed which allows clinicians and surgeons to explore and comprehend interactively the elastic response of patient-specific bone models in preparation of e.g. a hip-joint replacement. The steering platform supports an immediate physical response to any user interaction such as changing loads or inserting/adjusting different implants. The high responsiveness of the platform allows the identification and localization of critical regions in a bone-implant scenario compared to the physiological state to detect e.g. stress shielding phenomena which often lead to post-surgical late effects (Computing and Visualization in Science (14(5) 207–216).
in cooperation with
Chair for Computation in Engineering, TUM (Prof. Rank)
TUM Clinics rechts der Isar - Biomechanik (Prof. Gradinger)
Chair for Computer Graphics & Visualization, TUM (Prof. Westermann)
Siemens AG
Figure: von Mises stress distribution of the femur under tension and compression loading (a), principal stress lines: green - compression, purple - tension (b), principal stress lines - bone/implant model (c).
Figure: comparative stress visualization: increase (green) / decrease (red) compared with the physiological state for various implant sizes.
NURBS-version finite cell method
in cooperation with
Chair for Computation in Engineering, TUM (Prof. Rank)
Numerical Structural Analysis with Application in Ship Technology, TUHH (Prof. Düster)
The advent of isogeometric analysis (IGA) using the same basis functions for design and analysis constitutes a milestone in the unification of geometric modeling and numerical simulation. However, an important class of geometric models based on the CSG (Constructive Solid Geometry) concept such as trimmed NURBS surfaces do not fully support the isogeometric paradigm, since basis functions do not explicitly represent the boundary.
The finite cell method (FCM) is a high-order fictitious domain method, which offers simple meshing of potentially complex domains into a structured grid of cuboid cells, while still achieving exponential rates of convergence for smooth problems. Here, a finite cell extension to isogeometric analysis achieves a truly straightforward transfer of a trimmed NURBS surface into an analysis suitable NURBS basis, while benefiting from the favorable properties of the high-order and high-continuity basis functions.
Figure: NURBS patches as embedding domain model for two design variants of a brake disk.
eigenfrequencies: 19.4Hz (left), 18.6Hz (right)
eigenfrequencies: 3.7Hz (left), 3.4Hz (right)
Figure: modified Sordelis Lo NURBS model with holes.
Figure: convergence of the modified Scordelis Lo model for uniform k-refinement, comparison with a p-version finite cell solution, various number of elements.
A structure preserving QR-based acceleration method for medium to large size symmetric profile matrices
The proposed method is based on the well-known QR-method for dense matrices. The extension includes a profile preserving implementation which shows a highest degree of flexibility and reliability and is highly suited for the independent computation of any set of eigenvalues for medium and large scale profiled matrices. Furthermore, two effective, stable and numerically cheap extensions overcome the troublesome stagnation of the convergence in presence of poorly separated eigenvalues. A repeated preconditioning process in combination with Jacobi rotations in the parts of the matrix with the strongest convergence is used to improve significantly both, local and global convergence (Int. J. Numer. Meth. Engng 2009; 79:1419–1442).
Figure: Real-symmetric convex profile matrix indicating the preserved profile structure.
Figure: Convergence plots for a square Kirchhoff-plate model (N=47) illustrating the stepwise extension of a standard QR approach by a novel pre-conditioning strategy and a Jacobi correction step. The extensions allow a redcution of the computational effort by >35%.
Figure: Numerical effort of the extensions, pre-conditioner and Jacobi correction, with reference to the standard QR approach (shift & deflation version). All results refer to finite element discretizations of a square Kirchhoff-plate problem which can be considered as a severe test case due to a large number of multiple and clustered eigenvalues. Even for large scale matrices a reliable reduction of the computational effort is observed.
Figure: Residual norm of the standard eigenvalue problem solved with the QR and extended QR approach. The extensions do not compromise the quality of the computational result and even contribute to a slight improvement of the residual norm.
