Conference Agenda

We live in a world of exploding complexity with enormous challenges. Digital twins, tightly integrating the real and the digital world, are a key enabler for decision making in the context of complex systems. While the digital twin has become an intrinsic part of the product creation process, its true power lies in the connectivity of the digital representation with its physical counterpart.

To be able to use a digital twin scalable in this context, the concept of an executable digital twin has been proposed. An executable digital twin is a stand-alone and self-contained executable model for a specific set of behaviors in a specific context. It can be leveraged by anyone at any point in lifecycle. To achieve this, a broad toolset of mathematical technologies is required - ranging from model order reduction, calibration to hybrid physics- and data-based models.

In this presentation, we review the concept of executable digital twins, address mathematical key building blocks such as model order reduction, real-time models, state estimation, and co-simulation and detail its power along a few selected use cases.

The numerical simulation of nonlinear problems in contact mechanics and mixed-dimensional coupling poses significant challenges for high-performance computing (HPC). The considerable computational effort introduced by the evaluation of discrete contact or mixed-dimensional coupling operators requires an efficient framework that scales well on parallel hardware architectures and is, thus, suitable for the solution of high-fidelity models with potentially many million degrees of freedom. Another unsolved task is the efficient solution of the arising linear systems of equations on parallel computing clusters. In this study, we investigate the scalability of computational strategies and algebraic multigrid (AMG) solvers to address these challenges effectively and ultimately target a faster time-to-solution. Two problem classes of utmost practical relevance serve as prototypes for this endeavor.

First, we investigate the finite element analysis of nonlinear contact problems based on non-matching mortar interface discretizations. Mortar methods enable a variationally consistent imposition of almost arbitrarily complex coupling conditions but come with considerable computational effort for the evaluation of the discrete coupling operators (especially in 3D). We identify bottlenecks in parallel data layout and domain decomposition that hinder a truly efficient evaluation of the mortar operators on modern-day HPC systems and then propose computational strategies to restore optimal parallel communication and scalability. In particular, we suggest a dynamic load balancing strategy in combination with a geometrically motivated reduction of ghosting data. Using increasingly complex 3D examples, we demonstrate strong and weak scalability of the proposed algorithms up to 480 parallel processes. In addition to the computational strategies, we investigate the integration of tailored AMG preconditioners for the resulting saddle point-type linear systems of equations into a then fully scalable simulation pipeline.

Second, we present the mixed-dimensional interaction of slender fiber- or rod-like structures with surrounding solid volumes, thus leading to a 1D-3D approach that we refer to as beam-to-solid coupling. In terms of practical applications, natural and artificial fiber-reinforced materials (e.g., biological tissue, composites) come to mind. Again, not so different from contact problems, the main challenges on the way to a scalable computational framework lie in the efficient and parallelizable evaluation of the underlying mixed-dimensional coupling operators and the design of tailored AMG preconditioners. Here, we will focus on a novel physics-based block preconditioning approach based on AMG that uses multilevel ideas to approximate the block inverses appearing in the system. Eventually, we will assess the performance and the weak scalability of the proposed block preconditioner using large-scale numerical examples.

The treatment of cardiovascular diseases most often involves coronary stents. Even with drug-eluting stents, implantation can give rise to in-stent restenosis: endothelial denudation and overstretch injuries may result in uncontrolled tissue growth and formation of obstruction to the blood flow. Critical areas where such side effects occur highly depend on the shear stresses and drug distribution inside the artery. For this reason, the analysis of blood flow dynamics in stented arteries is of great interest. The current work is aimed at coupling hemodynamics and tissue growth to include the fluid-structure interaction of pharmacokinetics at the interface between artery and lumen.

Navier Stokes equations and Newtonian constitutive model are used to simulate blood in a stented artery. Wall shear stress (WSS) related quantities are analyzed as indicators of the possible areas of inflammation and thrombosis. Drug elution and deposition on the vessel wall is modeled by means of an advection-diffusion equation and tailored boundary conditions [1]. The convective field is obtained coupling the drug equation to a steady averaged blood flow over three heart beats. Since the healing process and drug elution span a time frame of weeks, a staggered approach is derived to simulate the drug release into the blood stream. Advection-diffusion-reaction equations form the basis of modeling the transport and interaction of species in the vessel wall. The corresponding equations for PDGF, TGF-ß, ECM and SMC can be found in [2]. The drug concentration field is coupled at the interface between the arterial wall and the lumen to account for downstream deposition of the drug. All governing equations for the wall species are coupled to a continuum mechanical description of volumetric growth.

In this work, we test our method on a simplified ring stent geometry with matching interface between the artery wall and the blood domain. We compare the effects of drug coupling and WSS on the endothelium and volumetric growth. All simulations are performed by means of finite element method using FEAP and the in-house code XNS.

[1] Hassler S, Ranno AM, Behr M. Finite-element formulation for advection–reaction equations with change of variable and discontinuity capturing. Computer Methods in Applied Mechanics and Engineering, 2020; 369: 113171.

[2] Manjunatha K, Behr M, Vogt F, Reese S. A multiphysics modeling approach for in-stent restenosis: Theoretical aspects and finite element implementation. Computers in Biology and Medicine, 2022; 150: 106166.

In-stent restenosis is one of the main adverse events after initially successful percutaneous coronary interventions (PCI) with stent implantation. Comprehensive statistical analyses of large clinical datasets identified several independent risk factors for restenosis occurrence, such as patient- or lesion-specific factors, which include small vessel size or the extended length of the stented section. However, it is widely accepted that the local mechanical state within the vessel wall strongly affects vascular growth mechanisms. Nevertheless, these biomechanical factors are currently not integrated into the predictive assessment of lesions at risk. For instance, high intramural stresses and overstretch of healthy vascular tissue during PCI may disturb the natural homeostasis and thus promote excessive tissue growth. Additionally, insufficient stent expansion and incomplete stent apposition reduce the long-term success rate of the procedure. We propose an individualized biomechanical model to study the influence of specific plaque characteristics on the mechanical state of the artery wall during loading conditions experienced in PCI and the final stent placement. In this work, we employ patient-specific artery models based on coronary computed tomography angiography data combined with resolved models of the stent delivery system for physics-informed PCI simulations. We define the system as a computational structural mechanics problem with large deformations and a nonlinear, viscoelastic material formulation for the artery considering the plaque constituents in a heterogeneous manner. The stent structure is resolved and is discretized with reduced-dimensional 1D Cosserat continua with an elastoplastic material formulation. An idealized inflatable balloon model governs the stent expansion. The interaction between balloon catheter and artery is modeled with computational contact mechanics using mortar methods; for the stents, we utilize a beam-to-solid contact approach. All simulations are performed with our in-house multiphysics high-performance code BACI, which uses finite element methods for all problem types considered here. We assess the local stresses and strains within the vessel wall during and after the stent implantation and collate cases with different lesion characteristics. We evaluate the contact between stent struts and endothelium for lesions at risk of incomplete stent apposition. Additionally, we compare the results of our resolved approach to a simplified model, where we model the stent as a pure cylinder with similar mechanical characteristics. In the future, insights from such modeling may inform the clinical assessment of lesions considered for stent implantation.

Neuroblastoma (NB) is the most frequent solid cancer of early childhood. It is a type of cancer that is highly representative of the cancer disease itself, since NB is strongly heterogeneous with very diverse clinical courses that may vary from an indolent disease causing little or no harm and exhibiting spontaneous regression, to an aggressive disease with fatal progression. For these reasons, NB is considered a paradigm of cancer disease and an excellent context of application for the validation of novel developments which have the ambition to be of potential application in a large variety of solid cancers.

NB tumours consist of two main cell populations, neuroblasts and Schwann cells, and the current neuroblastoma classification is based on histological criteria, e. g. the quantity of Schwannian stroma. Neuroblasts and Schwann cells are primary interest herein for contribute directly to the mechanical properties of the tissue through the proliferation and death processes. Extracellular matrix also have a principal role in the cell-microenvironmental cross-talk therefore the tumour can promote to a better stage or keep growing.

We here present a phenomenological model which takes into account as detail as possible to better mimic the real tumour behaviour. Our hypothesis proposes that tumour evolution can be attributed to three distinct processes: growth, shrinkage, and remodelling. The biomechanical model is based on the mass and cellular balance equations coupled with elasticity. The multispecies model simulates the effect of the cellular processes that occur during tumour growth and shrinkage, namely proliferation and death.

The biomechanical finite element model of NB tumour growth starts from imaging data derived mainly from MRI sequences. This data comprises the geometry, the initial cellularity distribution and the tumour vasculature evaluation. At the end of the simulation, the results obtained are validated with a second set of imaging data obtained after treatment.

The study simulates three-month chemotherapy using real patient cases, and presents two distinct outcomes: in one of them, the tumour volume was reduced 20% and in the other one, the volume decreased 90%. One of the patients was classified as low-risk, following the International Neuroblastoma Risk Group (INRG) system, whereas the other was classified as intermediate-risk. Differences appeared in the histology analysis, which reveal one tumour with a higher concentration of tumoural cells, and in the radiomic data obtained after image analysis. The model effectively reproduces these varying outcomes following the application of chemotherapy, facilitating the identification of cases in which the treatment may be effective.

Biomechanical models typically contain numerous parameters. Global sensitivity analysis helps identify the most influential and the non-influential parameters, as well as interactions between the parameters.

We show how to apply variance-based global sensitivity analysis to complex biomechanical models. As the method necessitates numerous model evaluations, we utilize Gaussian process metamodels [1] to lessen the computational burden. The approach is illustrated for models of active biomechanical systems by applying it to nanoparticle-mediated drug delivery in a multiphase tumour-growth model [2] and the formation of aneurysms in a model of aortic growth and remodelling [3].

We discover that a small number of full model evaluations suffices to effectively differentiate influential from non-influential parameters, while further evaluations enable the estimation of higher-order interactions. From a biomechanical modeling standpoint, we observe that often a few influential parameters predominantly govern the model output variance. Simultaneously, substantial parameter interactions can exist, emphasizing the necessity for global methods.

Gaussian process-based global sensitivity analysis proves feasible and beneficial for intricate, computationally demanding biomechanical models. Specifically, it can serve as a foundational building block for parameter identification.

[1] Le Gratiet L, Cannamela C, Iooss B. A Bayesian Approach for Global Sensitivity Analysis of (Multifidelity) Computer Codes. J Uncertain Quantif. 2, 336–363. DOI: 10.1137/130926869 (2014).

[2] Wirthl B, Brandstaeter S, Nitzler J, Schrefler BA, Wall WA. Global sensitivity analysis based on Gaussian-process metamodelling for complex biomechanical problems. Int J Numer Meth Biomed Engng. e3675. DOI: 10.1002/cnm.3675 (2023).

[3] Brandstaeter S, Fuchs SL, Biehler J, Aydin RC, Wall WA, Cyron CJ. Global Sensitivity Analysis of a Homogenized Constrained Mixture Model of Arterial Growth and Remodeling. J Elast. 145, 191–221 DOI: 10.1007/s10659-021-09833-9 (2021).

Flows with both compressible and (nearly) incompressible species appear in many engineering applications, such as cavitation near propellers or oil-dragging action by blow-by gases in internal combustion engine piston sealing rings[1].

Simulation of these systems can offer insight into the behavior of these systems in the early design stages of products or when investigating a large design space, where building prototypes is more expensive than simulating them.

Nevertheless, simulation requires accurate and efficient models for these complex flow phenomena.

The level-set method has been used for incompressible two-phase flow with success[2].

It offers intuitive computation of surface curvature based on the signed distance field.

For systems where one of the phases is compressible while the other is incompressible, material models have to be used that capture the behavior in both cases.

Also, the level-set method does not guarantee mass conservation, requiring special treatment of advection and re-initialization[3].

The models and numerical methods used for this will be presented, as well as

results for test and simple application cases.

[1] Benoist Thirouard and Tian Tian, “Oil transport in the piston ring pack (Part I): identification and characterization of the main oil transport routes and mechanisms”, Technical report, SAE Technical Paper, 2003

[2] Violeta Karyofylli, Loïc Wendling, Michel Make, Norbert Hosters, and Marek Behr, “Simplex space-time meshes in thermally coupled two-phase flow simulations of mold filling”, Computers & Fluids, 192 (2019) 104261.

[3] Elin Olsson and Gunilla Kreiss, “A conservative level set method for two-phase flow”, Journal of computational physics, 210 (2005) 225–246.

Facing phase-change systems ubiquitous in engineering and geophysical applications, we today leverage a range of problem-tailored numerical techniques. Monolithically formulated mathematical models, also referred to as fixed grid techniques [1], consider a single domain containing two or more phases with distinct material parameters, but impose the same governing equations everywhere. It is well known that numerical difficulties arise from the steeply varying material behaviour and the resulting strong non-linearities at the phase-change interface (PCI) [2].

At sufficiently large spatial scales, mixture solidification problems, such as the solidification of alloys or saltwater systems, feature a continuous mushy-layer transition from pure solid to pure liquid. At thermodynamic equilibrium, fixed grid techniques considering just one energy conservation equation for all phases are the de facto standard. While these models theoretically reduce to the classical Stefan problem as the mixture impurity tends to zero, they are often challenging to solve numerically. This is because the implicitly assumed continuity of the mixture PCI, acting as a regularization of the interface non-linearity, is lost in the limit of a sharp Stefan-problem type interface.

In this contribution, we contextualize numerical approaches to single-equation models with applications to solidification of pure substances and mixtures. In particular, we discuss the choice of primary variables in the energy equation, providing a comparative study of enthalpy-based and temperature-based formulations [3]. We examine the two approaches in terms of ease of implementation, accuracy and computational effort and provide a reference simulation based on a level-set method [4]. The methods will be applied both to a solidifying pure substance, where the interface propagation can be interpreted as a Stefan problem, and a binary mixture with a continuous mushy type phase-transition region. Based on our numerical experiments, we will conclude on guidelines to picking the most efficient formulation for the considered problem.

[1] Voller et al., "Fixed grid techniques for phase change problems: A review", 1990

[2] Krabbenhoft et al., "An implicit mixed enthalpy–temperature method for phase-change problems", 2007

[3] Terschanski et al.,: "Reactive Transport Models for Ice-Ocean Interfaces", Poster presented at SIAM CSE Amsterdam, Feb. 28, 2023

[4] Boledi et al., "A level-set based space-time finite element approach to the modelling of solidification and melting processes", 2022

Phase transition problems in the setting of non-equilibrium thermodynamics appear in various industrial and academical problems. More specifically, sublimation and deposition phenomena involving rarefied gases are important processes in freeze drying in the pharmaceutical area or in the behavior of planetary atmospheres. In these cases, the low-pressure regime is the reason for the rarefaction and consequently for the non-equilibrium behavior.

Modeling non-equilibrium behavior of rarefied gases is challenging as default continuum models can no longer accurately describe them. Instead, the kinetic theory is consulted to derive suitable descriptions.

The Stefan problem is a classical model for phase change problems, originally designed for solid-liquid interactions. We are generalizing this model for solid-gas interaction, so for sublimation and deposition problems. The gas phase is treated rarefied, where the non-equilibrium effects are introduced. A dependency on the well-studied Knudsen and Mach numbers, defined at the phase transition interface, is established. Depending on the level of the rarefaction, the resulting differences from the classical model are significant.

Motivated by the vast number of applications of higher-order partial differential equations representing gradient flow like the Cahn-Hilliard equation, the Allen-Cahn equation, the phase-field crystal equation, and the Swift–Hohenberg equation, among others, in this contribution we propose an efficient combination of the Preconditioned Conjugated Gradient (PCG) solver and the recently proposed Iterative Sherman-Morrison Inversion (ISMI) to solve the systems of linear equations that arise during an implicit-time integration of these equations if Fourier-spectral methods are employed for the spatial discretization. PCG is computationally expensive when compared to ISMI, which has a superior convergence, especially during the first few iterations but its computational edge over PCG is lost if too many solver iterations are carried out due to a higher storage demand. Therefore, in this work we propose to selectively combine PCG and ISMI solvers such that the advantages of both solvers are exploited at different stages of the solution scheme which improves the convergence of the residual error of the linear system and thereby the computational efficiency considerably in comparison to standalone versions of the solvers. Some numerical examples are presented in the context of all the aforementioned types of gradient flow in 2D and 3D to demonstrate the benefits of the selective combination in the context of different phase distributions.

Every material in nature exhibits heterogeneous behaviour at a certain scale. In a system, defects such as pores, grain boundaries, phase boundaries, secondary phases and particles can be the reasons for heterogeneity. The effective behaviour of the materials is significantly influenced by the underlying microstructure. Interfaces, such as grain boundaries, can affect the overall response of the material under consideration. Experimental findings shows that grain boundaries have a critical influence on electrical properties [1] and in order to model the macroscopic behaviour realistically, interfaces at the microscale should be taken into account.

Motivated by the change of effective electrical properties due to interfaces, e.g. microcracks or grain boundaries, a computational multiscale framework for continua with interfaces at the microscale is proposed in this contribution. More specifically speaking, the computational multiscale formulation for electrical conductors [2] is extended to account for interfaces at the microscale. Cohesive-type interfaces are considered at the microscale, such that displacement and electrical potential jumps can be accounted for. The governing equations for the materials with interfaces under mechanical and electical loads are provided. Based on these, a computational multiscale formulation is established. In particular, averaging theorems for kinematic quantities and for their energetic duals are discussed and their consistency with an extended Hill-Mandel condition for suitable boundary conditions is shown. The coupling between the electrical and mechanical subproblem is established by the constitutive equations at the material interface. In order to investigate deformation-induced property changes at the microscale, evolution of interface damage is elaborated.

To show the capabilities of the proposed framework, different representative simulations are selected. In particular, the calculation of effective macroscopic conductivity tensors for given two-dimensional microstructures is discussed and the fully coupled effective electro-mechanical material response due to the damage evolution is presented.

References

[1] H. Bishara, S. Lee, T. Brink, M. Ghidelli, and G. Dehm, “Understanding grain boundary electrical resistivity in Cu: The effect of boundary structure,” ACS Nano, vol. 15, no. 10, pp. 16607–16615, 2021.

[2] T. Kaiser and A. Menzel, “An electro-mechanically coupled computational multiscale formulation for electrical conductors,” Arch. Appl. Mech., vol. 91, pp. 1–18, 2021.

The technical relevance of small-scale electromechanical systems is rapidly increasing today. For this reason, the flexoelectric effect, which occurs in all dielectrics, is increasingly getting into the focus of research. This size-dependent effect describes the linear coupling between the electric polarization in the material and an occurring strain gradient in, for example, bent cantilever beams. There also exists a converse flexoelectric effect defined as a mechanical stress response under the action of an electric field gradient especially noticeable at sharp electrode tips in microelectromechanical systems (MEMS). In order to make these coupling effects technically usable, suitable models are required to predict the resulting system response.

A continuum-based model approach that takes into account elastic, dielectric, piezoelectric and flexoelectric effects is presented. Different model variants will be discussed and suitable finite element formulations for solving the electromechanical boundary value problem will be presented. A mixed variation formulation is used here in order to reduce the higher continuity requirements due to the occurring gradient fields. When considering ferroelectric materials (e.g. PZT), microstructural domain switching processes must be taken into account in order to be able to predict the behavior realistically. A microscopically motivated material model representing these dissipative processes is introduced and fitted into the flexoelectric continuum approach. The influence of acting strain and electric field gradients on the domain switching processes in ferroelectrics when considering the flexoelectric effect is studied by numerical experiments.

Inspired by lattice structures that can be observed in nature, periodic unit cells and their mechanical properties have caused an ever increasing interest in recent years due to the growing performance of additive manufacturing methods. In order to incorporate cells with optimal properties into printed high-performance structures and devices that can respond to given macroscopic stress-strain states in an optimal manner, one has to provide anisotropic properties that can respond to these individual loads.

We discuss the performance of a phase-field approach for optimizing periodic micro-structures based on triply periodic minimal surface problems (TPMS) to obtain unit cells with an optimal homogenized stiffness response in the direction of the maximal principal stress direction. We show that different TPMS-types exhibit fundamental differences in the way they can respond to uni-axial or shear-dominated loads. An essential aspect in optimizing cells is, on the one hand, to maximize the compliance with external loads and, on the other hand, to limit the danger of failure due to local buckling which is achieved by preserving the connectivity of the cell grid.

Further aspects that are discussed include numerical strategies to handle linear systems of such high-resolution optimization problems in an efficient manner as well as strategies to verify the gain of the homogenized stiffness experimentally.

The Bouligand structure comprises twisted parallel fibers arranged in a helical pattern, which enables greater energy dissipation and fracture toughness, mainly through a large crack surface area and crack-bridging phenomenon compared to regular fiber-reinforced composites. Considering the complex nature of this structure, numerical models that accurately capture the propagation of cracks through its twisted fiber arrangement are limited. This is due to the most popular simulating approach, the finite element method (FEM), is based on the classic continuum mechanics that uses spatial differential equations to describe continuous material behaviors. In contrast, Peridynamics is a computational framework that has been developed to overcome the limitations of classical continuum mechanics in describing crack propagation. Unlike FEM, Peridynamics is a non-local continuum theory that utilizes integral equations instead of differential equations in space to simulate material behaviors. This characteristic makes it highly suitable for modeling the complex crack propagation behaviors in the Bouligand structure. In this study, we present a bond-based peridynamics model to accurately describe the fiber-reinforced composites with a small angle mismatch between adjacent layers in Bouligand structures. To investigate the fracture mechanisms of such a structure, we conduct comprehensive numerical simulations, including 3-point bending and low-velocity impact tests, to obtain detailed information on its deformation and failure behavior. This information is difficult to achieve solely through experimental and theoretical studies. Based on our insights into the toughening mechanisms of the Bouligand structure, we propose a novel approach to further enhance the material’s fracture toughness by combining the Bouligand structure with other toughening mechanisms. Overall, the current study provides important insights into the fracture behavior of Bouligand structures and presents new avenues for designing advanced materials with superior mechanical properties.

Reduced Order Models (ROMs) have been widely used to efficiently solve large-scale problems in many fields including computational fluid dynamics (CFD) [1]. ROMs techniques allow to replace the expansive Full Order Model (FOM), by a ROM that captures the essential features of the system while significantly reducing the computational cost. In this work, we draw inspiration from [2] to implement a reduced basis (RB) method for model reduction of the Shallow Water Equations (SWEs) using Proper Orthogonal Decomposition (POD) and Artificial Neural Networks (ANNs). This method, referred to as POD-NN, starts with the POD technique to construct a reduced basis, and then makes us of an ANN to learn the associated coefficients in the reduced basis. It follows an offline-online strategy: the POD reduced basis along with the training of the ANN are performed in an offline stage, and then the surrogate model can be used for hyper-fast predictions. The process is non-intrusive since it does not require opening the black box of the FOM. The developed method is tested [3] on a real data set aiming at simulating an inundation of the Aude river (Southern France). The results show that the proposed method can achieve significant computational savings while maintaining satisfactory accuracy on the hydraulic variables of interest compared to the full-order hydraulic model. The proposed method is able to capture the key features of the SWEs in particular the wave propagation. Overall, the proposed non-intrusive POD-NN method offers a promising approach for ROM of SWEs while being affordable in view of fast real time inundation simulations.

[1] Benner, P., Schilders, W., Grivet-Talocia, S., Quarteroni, A., Rozza, G., & Miguel Silveira, L. (2020). Model Order Reduction: Volume 2: Snapshot-Based Methods and Algorithms (p. 348). De Gruyter.

[2] Hesthaven, J. S., & Ubbiali, S. (2018). Non-intrusive reduced order modeling of nonlinear problems using neural networks. Journal of Computational Physics, 363, 55-78.

[3] IMT-INSA, INRAE et al., “DassFlow: Data Assimilation for Free Surface Flows”, Open source computational software. https://www.math.univ-toulouse.fr/DassFlow

Let us consider wave propagation problems over 2-dimensional domains with piecewise-linear boundaries, possibly including scatterers. Under some assumptions on the initial conditions and forcing terms, we have proposed an approximation of the propagating wave as the sum of some special nonlinear space-time functions. Each term in this sum identifies a particular ray, modeling the result of a single reflection or diffraction effect. In this talk I will describe an algorithm for identifying such rays automatically, based on the domain geometry.

To showcase our proposed method, I will present several numerical examples, such as waves scattering off wedges and waves propagating through a room in presence of obstacles.

This contribution is motivated by the combined advantages of an integrated framework from CAD geometries to simulation in real time. In a typical workflow for design and shape optimization, multiple simulations are required for all possible designs represented by different geometrical parameters. This might entail a high computational cost in particular for real world, engineering applications. The development of efficient reduced order models that enable fast parametric analysis is essential for such applications. At the same time, the capabilities of splines and isogeometric analysis allow for flexible geometric design and higher-order continuity in the analysis. In CAD design, trimmed multi-patch geometries are widely used to represent complex shapes. The presence of geometric parameters introduces challenges for efficient reduced order modeling of problems formulated on such unfitted geometries. We propose a localized reduced basis method to circumvent the shortcomings of standard reduced order models in this context [1]. In this talk we present the developed strategy and address fast parametric analysis of problems in structural mechanics. The construction of efficient reduced order models for geometries described by multiple trimmed patches as well as their use in parametric shape optimization will be discussed. Numerical examples will be presented to demonstrate the accuracy and computational efficiency of the method.

[1] M. Chasapi, P. Antolin, A. Buffa, A localized reduced basis approach for unfitted domain methods on parameterized geometries, Comput. Methods Appl. Mech. Engrg. 410 (2023) 115997.

In many applications in Computer Aided Engineering, like parametric studies, structural optimization or virtual material design, a large number of almost similar models have to be simulated. Although the individual scenarios may differ only slightly in both space and time, the same effort is invested for every single new simulation with no account for experience and knowledge from previous simulations. Therefore, we have developed a method that combines data-based Model Order Reduction (MOR) and reanalysis, thus exploiting knowledge from previous simulation runs to accelerate computations in multi-query contexts. While MOR allows reducing model fidelity in space and time without significantly deteriorating accuracy, reanalysis uses results from previous computations as a predictor or preconditioner.

The workflow of our method, named Reduced Model Reanalysis (RMR), is divided into an offline and online phase. In the offline phase, data are generated to cover a wide range of the parameter space. From this data a surrogate model is learned in a reduced space using regression algorithms from the field of machine learning. Depending on the requirements of the system, different regression algorithms are favorable, e.g. linear regression, a k-nearest neighbor algorithm, a neural network, or a Gaussian process. The models are learned in the reduced space due to the prohibitively large number of degrees of freedom of the full finite element model. The reduced subspaces are obtained via a snapshot POD (proper orthogonal decomposition). In the online phase, approximations of all relevant solution quantities are obtained from the surrogate model. Their projection to the full space provides predictors that allow for an accelerated solution of the system in comparison to a standard structural mechanics computation.

In the case of nonlinear stability analysis this method can for example be used to accelerate the exact computation of critical points by the method of extended systems. Data generation in the offline phase is also accelerated by a newly developed adaptive time stepping scheme. With this scheme the number of steps to approach critical points with a path following scheme can be significantly reduced. Further potential fields of application of RMR are general nonlinear static and transient problems, with particular challenges as soon as path-dependence comes into play.

Bones have the ability to adapt their structure and thus their density to external loads. Cancellous bone, which forms the spongy interior of bones, is a microstructural network of rods. Under- or overloading strengthens or narrows these rods, altering the microstructural pattern. In this adaption process, osteocytes act as mechanosensors, activated by mechanical signals and regulating the mechanical adaptation of bone. That is, they communicate with bone-forming or bone-resorbing cells. Thus, bone remodelling at a particular point is triggered by non-local mechanosensors in its vicinity, i.e. the sensors involved act in a specific sphere of influence and not only locally.

In this work, we present a micromorphic approach that extends the established concept of local bone adaption to account for both the non-locality of bone remodelling and the heterogeneous structure of the material without explicitly resolving it within a two-scale approach. Our approach enables a simple implementation in the open source finite element environment deal.II and avoids the need for laborious neighborhood sampling, as is the case with integral approaches, or for higher continuity requirements, as is the case with higher gradient approaches.

Our approach is phenomenological in nature and refers to nominal bone density to be interpreted as a macroscopic measure of the ratio of bone mass to pore volume in the underlying trabecular microstructure. This way, we account for the heterogeneous microstructure of bone by capturing its effect on nominal bone density, but without actually resolving individual trabeculae. Since bone is a living material, in the continuum approach to bone remodelling we apply the theory of open-system thermodynamics, which assumes that there is a mass source corresponding to the change in nominal density over time. The mass source is equated with a mechanical stimulus, comparing the stored energy to an attractor. The attractor can be interpreted as a biological stimulus that drives remodelling. In the local case, the stored energy is a purely local quantity that depends on the macroscopic deformation. In our novel non-local approach, we now extend this by adding a micromorphic and a scale-bridging component to the stored energy. This allows us to account for non-locality with a characteristic length scale, which acts as a measure for the heterogeneous microstructure and a scale-bridging parameter that penalizes the deviation of the micromorphic from a higher gradient model.

The approach is illustrated in depth and its implications are discussed using benchmark examples. In addition, the modeling approach is discussed using long tubular bones and compared with CT images in health and osteoporosis.

Gastric peristalsis refers to the coordinated contraction and relaxation of the muscles in the stomach wall that mixes and grinds food and propels chyme down the digestive tract. Gastric peristalsis is realized by an intricate electromechanical system. We present a computational multiphysics framework for its simulation on patient-specific stomach geometries. It combines a robust gastric electrophysiology model with an active-strain finite elasticity model for the tissue mechanics [1,2,3]. The patient-specific spatially varying parameter distributions are determined by a novel algorithm mapping a two-dimensional parameter distribution function onto a general tube-like surface. The capability of the proposed computational framework for large-scale in silico analyses of gastric electromechanics is demonstrated on patient-specific human stomach models derived from magnetic resonance images. The proposed framework can reproduce essential phenomena on patient-specific stomach geometries, including the entrainment and propagation of stable ICC slow waves as well as the propagation of physiological ring-shaped peristaltic contraction waves. In summary, the presented framework provides a powerful tool to study gastric electromechanics in health and disease. This can enable optimized patient-specific diagnosis and therapy planning.

References

[1] Brandstaeter, S., et al., Computational model of gastric motility with active-strain electromechanics. ZAMM - Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik, 2018.

[2] Djabella, K., M. Landau, and M. Sorine, A two-variable model of cardiac action potential with controlled pacemaker activity and ionic current interpretation, in 2007 46th IEEE Conference on Decision and Control. 2007, Institute of Electrical & Electronics Engineers (IEEE). p. 5186 - 5191.

[3] Ruiz-Baier, R., et al., Mathematical modelling of active contraction in isolated cardiomyocytes. Mathematical Medicine and Biology, 2014. 31: p. 259-283.

Human brain tissue shows complex, nonlinear, and time dependent mechanical behavior, and thus presents a significant challenge to those interested in developing accurate constitutive models. Research in our group is focused on better understanding the factors that influence the mechanical response of the tissue. To this end, we combine the mechanical testing of tissue samples from different brain regions under finite deformation in compression, tension and shear with microstructural analyses, continuum mechanics modeling, and finite element simulations. The application of an inverse parameter identification allows us to determine material parameters with a subsequent statistical analysis revealing their regional dependence. Here, we find that the corpus callosum and corona radiata in particular have to be considered as regions with distinct mechanical properties when modeling the whole brain. Furthermore, we analyze the protein content of the tested specimens by enzyme-linked immunosorbent assays and show their correlation with the identified material parameters. These results may motivate and guide the development of microstructurally informed constitutive models that may enable patient-specific predictions.

Traumatic brain injury (TBI) brought on by a severe head impact in a car accident, a fall, or a sports injury results in internal tissue damage beyond recovery. The human brain mainly has two vital tissues; gray matter and white matter. During accidental impact, forces and torques are imparted in the brain tissues to trigger significant local damage. Although the brain can recover from a TBI, the force necessary to cause permanent brain damage is still not fully understood. One aspect of investigating TBI is to provide a mathematical model and a computational framework to identify the level of injury. Mechanical characterization of the brain tissue is essential to understand brain damage caused by TBI. Since 1960, many studies have been done to understand the brain’s mechanical behavior. It is found that brain tissue’s behavior is an incompressible, viscoelastic material and anisotropic material. The human brain’s finite element (FE) models have been utilized to investigate the risk and mechanisms of traumatic brain injuries. Many human brain FE models have been developed. Many different constitutive models have been used for different parts of the human head. Still, there is scope for improvement in constitutive modeling and its finite element implementation. In this work, we present an anisotropic viscoelastic constitutive model and essential equations for finite element simulations. We implemented the constitutive model by ABAQUS UMAT for doing finite element simulations of the human head FE-model for real-life loading cases. Our uniaxial and cyclic loading simulation using UMAT agree with experimental and MATLAB results.

Ferroelectric as well as ferromagnetic materials are widely used in smart structures and devices as actuators, sensors etc. Regarding their nonlinear behavior, a variety of models has been established in the past decades. Investigating hysteresis loops or electromechanical/magnetoelectric coupling effects, only simple boundary value problems (BVP) are considered. In [1] a new scale–bridging approach is introduced to investigate the polycrystalline ferroelectric behavior at a macroscopic material point (MMP) without any kind of discretization scheme, the so–called Condensed Method (CM). Besides classical ferroelectrics, other fields of application of the CM have been exploited, e.g. [2, 3, 4]. Since just the behavior at a MMP is represented by the CM, the method itself is unable to solve complex BVP, which is technically disadvantageous if a structure with e.g. notches or cracks shall be investigated.

In this paper, a concept is presented, which integrates the CM into a Finite Element (FE) environment. Considering the constitutive equations of a homogenized MMP in the weak formulation, the FE framework represents the polycrystalline behavior of the whole discretized structure, which finally enables the CM to handle arbitrary BVP. A more sophisticated approach completely decouples the constitutive evolution from the FE discretization, by introducing an independent material grid. Furthermore, energetic consistencies of scale transitions from grain to MMP and MMP to macroscale are investigated. Numerical examples are finally presented in order to verify the approach.

References

[1] Lange, S. and Ricoeur, A., A condensed microelectromechanical approach for modeling tetragonal ferroelectrics, International Journal of Solids and Structures 54, 2015, pp. 100 – 110.

[2] Lange, S. and Ricoeur, A., High cycle fatigue damage and life time prediction for tetragonal ferroelectrics under electromechanical loading, International Journal of Solids and Structures 80, 2016, pp. 181 – 192.

[3] Ricoeur, A. and Lange, S., Constitutive modeling of polycrystalline and multiphase ferroic materials based on a condensed approach, Archive of Applied Mechanics 89, 2019, pp. 973 – 994.

[4] Warkentin, A. and Ricoeur, A., A semi-analytical scale bridging approach towards polycrystalline ferroelectrics with mutual nonlinear caloric–electromechanical couplings, International Journal of Solids and Structures 200 – 201, 2020, pp. 286 – 296.

Components used in the aerospace or automotive industries are often exposed to multi-physical loading conditions and thus may simultaneously be subjected to high stresses and strains as well as temperature changes. Therefore, high-strength and high-temperature resistant materials such as metals are commonly used for applications in this field. Since the overall material behavior is directly influenced by the distribution, size and morphology of the individual grains of the underlying polycrystalline microstructure, detailed knowledge of this microstructural behavior is required in order to accurately predict the macroscopic material response. Hence, multi-scale simulation approaches have been developed. Considering a two-scale finite element (FE) and fast Fourier transform (FFT)-based simulation approach [1, 2], the macroscopic and microscopic boundary value problems are first solved individually by assuming scale separation. In this context, the homogeneous macroscale is subdivided into a discrete number of finite elements. The microscopic boundary value problem is attached to each macroscopic integration point and solved using the FFT-based simulation approach. The scale transition is then performed by defining the macroscopic quantities as the average value over the corresponding local fields. This simulation approach is an efficient alternative to the classical FE² method for the simulation of periodic unit cells [3]. To illustrate the applicability of our model, we will present several numerical examples.

[1] J. Spahn, H. Andrä, M. Kabel, and R. Müller. A multiscale approach for modeling progressive damage of composite materials using fast Fourier transforms. Computer Methods in Applied Mechanics and Engineering, 268, 871–883, 2014

[2] J. Kochmann, S. Wulfinghoff, S. Reese, J. R. Mianroodi, and B. Svendsen. Two-scale FE–FFT- and phase-field-based computational modeling of bulk microstructural evolution and macroscopic material behavior. Computer Methods in Applied Mechanics and Engineering, 305, 89–110, 2016

[3] C. Gierden, J. Kochmann, J. Waimann, B. Svendsen, and S. Reese. A review of FE-FFT-based two-scale methods for computational modeling of microstructure evolution and macroscopic material behavior. Archives of Computational Methods in Engineering, 29(6), 4115-4135, 2022.

Numerical simulation of complex geometries and microstructures can be costly and time consuming, in particular due to the long process of preparing the geometry for meshing and the meshing process itself [1]. Several methods were proposed to overcome this issue, such as the extended finite element, meshless, Fourier transform and immersed boundary methods. Immersed boundary methods rely on embedding the physical domain into a Cartesian grid of finite elements and resolving the geometry only by adaptive numerical integration schemes. For instance, the isogeometric finite cell method (FCM) exploits the accuracy of higher-order, smooth B-Spline basis functions for the discretization and employs an octree scheme in order to refine the quadrature rule in trimmed elements. FCM has been applied successfully to various problems in solid mechanics, including linear and nonlinear elasticity, elasto-plasticity, and thermo-elasticity [2]. However, FCM typically requires several levels of refinement of the quadrature rule in order to deliver accurate results, which may lead to high computation times, especially for nonlinear, internal variable, and coupled multiphysics problems.

In this work, we adopt a novel algorithm for boundary-conformal quadrature based on a high-order reparameterization of trimmed elements [3] to solve small and large deformation thermo-elastic problems using spline-based immersed isogeometric analysis (IGA) without the need for a body conformal finite element mesh. In particular, the Gauss points on trimmed elements are obtained by a NURBS reparameterization of the physical subdomains of the cut elements of the Cartesian grid. This ensures an accurate integration with a minimal number of quadrature points. Furthermore, using periodic B-Spline discretizations, periodic boundary conditions for homogenization can be automatically fulfilled. Several numerical examples are presented to show the accuracy and efficacy of the boundary-conformal quadrature algorithm.

REFERENCES

[1] T. Hughes, J. Cottrell, and Y. Bazilevs. Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering, 194:4135– 4195, 2005.

[2] Schillinger, D. and Ruess, M., 2015. The Finite Cell Method: A review in the context of higher-order structural analysis of CAD and image-based geometric models. Archives of Computational Methods in Engineering, 22(3), pp.391-455.

[3] Wei, X., Marussig, B., Antolin, P. and Buffa, A., 2021. Immersed boundary-conformal isogeometric method for linear elliptic problems. Computational Mechanics, 68(6), pp.1385-1405.

Metamaterials are attracting growing attention in industry and academia due to their unique mechanical behaviour. However, when the scale separation does not hold, they show size-effects. Generalized continua can model such materials as a homogeneous continuum with capturing the size-effects.

The relaxed micromorphic model [1] describes the kinematics of each material point via a displacement vector and a second-order micro-distortion field. It has demonstrated many advantages over other higher-order continua such as using fewer material parameters and the drastically simplified strain energy compared to the classical micromorphic theory. Moreover, the relaxed micromorphic model operates between two bounds; linear elasticity with the micro and macro elasticity tensors. The strain energy function in the relaxed micromorphic model employs the Curl of the micro-distortion field and therefore H(Curl)-conforming FEM implementation is necessary [2-3].

In our talk, we will present our recent results in identifying the material parameters and boundary conditions in the relaxed micromorphic model [4].

REFERENCES

[1] P. Neff, I.D. Ghiba, A. Madeo, L. Placidi and G. Rosi. A unifying perspective: the relaxed linear micromorphic continuum. Continuum Mechanics and Thermodynamics 26,639-681(2014).

[2] J. Schröder, M. Sarhil, L. Scheunemann and P. Neff. Lagrange and H(curl,B) based Finite Element formulations for the relaxed micromorphic model, Computational Mechanics 70, pages 1309–1333 (2022).

[3] A. Sky, M. Neunteufel, I. Muench, J. Schöberl, and P. Neff. Primal and mixed finite element formulationsfor the relaxed micromorphic model. Computer Methods in Applied Mechanics and Engineering 399, p. 115298 (2022).

[4] M. Sarhil, L. Scheunemann, J. Schröder, P. Neff. Size-effects of metamaterial beams subjected to pure bending: on boundary conditions and parameter identification in the relaxed micromorphic model. https://arxiv.org/abs/2210.17117 (2022).

In the present contribution, we deal with a second-order computational homogenization of fluid flow in porous materials. Similar to the first-order computational homogenization in [1], the microscopic problem is formulated employing a minimization-type variational formulation at small strains; see also [2]. While a first-order Darcy-Biot-type fluid transport is considered at the microscale [2], the macroscopic problem is characterized by a second-order material response [3]. Hence, the present formulation allows the relaxation of the scale-separation assumption and the incorporation of the macroscopic second-order terms associated with deformation and fluid-flux fields at the microscale. The macro- and microscale boundary value problems are then bridged via an extended form of the Hill-Mandel condition, which results in suitable boundary conditions at the microscale and a set of constraints [4,5]. Finally, we present numerical examples that provide further insights into the presented formulation.

References:

[1] E. Polukhov and M.-A. Keip. Computational homogenization of transient chemo-mechanical processes based on a variational minimization principle. Advanced Modeling and Simulation in Engineering Sciences, 7, 1-26 (2020).

[2] C. Miehe, S. Mauthe, and S. Teichtmeister. Minimization principles for the coupled problem of Darcy--Biot-type fluid transport in porous media linked to phase field modeling of fracture. Journal of the Mechanics and Physics of Solids, 82, 186-217 (2015).

[3] G. Sciarra, F. dell'Isola, and O. Coussy. Second gradient poromechancis. International Journal of Solids and Structures, 44, 6607-6629 (2007).

[4] V.G. Kouznetsova, M.G.D. Geers and W.A.M. Brekelmans. Multi-scale second-order computational homogenization of multi-phase materials: a nested finite element solution strategy. Computer methods in Applied Mechanics and Engineering, 193, 5525-5550 (2020).

[5] I. A. Rodrigues Lopez, and F. M. Andrade Pires. Unlocking the potential of second-order computational homogenisation: An overview of distinct formulations and a guide for their implementation. Archives of Computational Methods in Engineering, 1-55 (2021).

Rayleigh-Benard convection (RBC) at high Rayleigh (Ra) numbers represents one of the most important model systems to study turbulent convection [1]. Experiments reaching very high Ra, approaching values relevant for convective systems in Nature, like the atmospheric convection, have been performed using various working uids, prominently with cryogenic helium 4He [2] and sulphur hexafluoride SF6 [3]. The goal of attaining high Ra often comes at the cost of breaking the Oberbeck-Boussinesq (OB) conditions at the phase boundaries or near the critical points of the working fluids. In particular, the recent analysis [4] of RBC experiments performed near the saturated vapor curves (SVC) in 4He and SF6 indicates that the heat transport measurements of the Nusselt number Nu(Ra), which apparently show the transition of RBC to the ultimate Kraichnan regime, are significantly affected by non-OB (NOB) effects, thus keeping the question of experimental observation of the ultimate regime open. The present study investigates the NOB effects which arise due to the temperature dependence of material properties in cryogenic helium experiments of turbulent RBC. The material properties such as specific heat at constant pressure, dynamic viscosity, thermal conductivity, the isobaric expansivity, and the mass density are expanded into power series with respect to temperature up to the quadratic order with coeffcients obtained from the software package HEPAK. A subsequent nonlinear regression that uses deep convolutional networks delivers a dependence of the strength of NOB effects in the pressure{temperature parameter plane. Strength of the NOB effects is evaluated via the deviation of the mean temperature profile $\xi_{NOB} \equiv T_m-T_c$ from the top/bottom-symmetric OB case $\xi_{NOB} = 0$. Training data for the regression task are obtained from 236 individual long-term laboratory measurements at different Rayleigh numbers which span 8 orders of magnitude. The work has been supported by the Czech Science Foundation project No. 21-06012J.

References

[1] G. Ahlers, S. Grossmann and D. Lohse, Rev. Mod. Phys., 81, 503, 2009.

[2] L. Skrbek and P. Urban, Journal of Fluid Mechanics, 785, 270, 2015.

[3] X. He, D. Funfschilling, H. Nobach, E. Bodenschatz and G. Ahlers, Phys. Rev. Letters, 108, 024502, 2012.

[4] P. Urban, P. Hanzelka, T. Kralk, M. Macek, V. Musilova and L. Skrbek, Phys. Rev. E, 99, 011101(R), 2019.

[5] D. D. Gray and A. Giorgini, Int. J. Heat Mass Transfer, 19, 545, 1976.

Stratified turbulent flows abound in environmental and industrial settings. Examples are atmospheric boundary layer flows, the transport of nutrients and organisms and the mixing of heat and salinity in the oceans, fluid flow in heat exchangers, and the transport of reactants and products in chemical reactions. These examples and many others consider stratified wall-bounded turbulence, in which the creation of turbulence by mechanical processes contends with its dissipation due to buoyancy effects. The buoyancy effects alter the structure of the flow, and consequently the dynamics of mass, heat, and momentum transport. As density fluctuations become more severe, the Oberbeck-Boussinesq approximation becomes inaccurate and the resulting dynamics are not correctly predicted. In the current work, we developed and validated a numerical solver for direct numerical simulations (DNS) of turbulent flows featuring strong property variations. More precisely, we solve the Navier-Stokes equations in the limit of vanishing Mach number (so-called low-Mach number limit), with the fluid density given by the ideal gas law, and the dynamic viscosity and thermal conductivity given by Sutherland's law.

Our numerical solver was then used to study stably-stratified turbulent channel flow under non-Oberbeck-Boussinesq conditions. The simulations will be carried out at friction Reynolds number of 395, Prandtl number of 0.71, and shear Richardson number in the O(10), where the friction Reynolds, Prandtl, and friction Richardson numbers are governing parameters defined based on the prescribed pressure drop and properties of the fluid at the reference temperature. Stratification is achieved by imposing constant temperature boundary conditions, with a high upper-to-lower wall temperature ratio (larger than 2), resulting in strong density variations in the flow. We will vary the temperature ratios and adjust gravity to maintain a similar Richardson number between cases, thereby isolating the effects of strong property variations in the flow dynamics. In the presentation, we will analyze the dynamics of heat and momentum transport under strong stratification for these conditions, also in light of DNS data of the same system under the Oberbeck-Boussinesq regime.

Transport phenomena in high Reynolds number wall-bounded stratified flows are dominated by the interplay between the turbulence structures generated at the wall and the buoyancy-induced large scale waves populating the channel core. In this study, we want to investigate the flow physics of wall-bounded stratified turbulence at relatively high shear Reynolds number Ret and for mild to moderate stratification level (quantified here by the shear Richardson number varying in the range 0<Rit<300). By increasing stratification, active turbulence is sustained only in the near-wall region, whereas intermittent turbulence, modulated by the presence of non-turbulent wavy structures (Internal Gravity Waves, IGW), is observed at the channel core. In such conditions, the wall-normal transport of momentum and heat is considerably reduced compared to the case of non-stratified turbulence. A careful characterization of the flow-field statistics shows that, despite temperature and wall-normal velocity fluctuations are very large at the channel center, the mean value of their product (buoyancy flux) vanishes for Rit>200. We show that this behavior is due to the presence of a pi/2 phase delay between the temperature and the wall-normal velocity signals: when wall-normal velocity fluctuations are large (in magnitude), temperature fluctuations are almost zero, and viceversa. This constitutes a blockage effect to the wall-normal exchange of energy. In addition, we show that the friction factor scales as a power of the Richardson number (-1/3), and we propose a new scaling for the Nusselt number (as a function of Reynolds and Richardson number). These scaling laws, which seem robust over the explored range of parameters, are expected to help the development of improved models and parametrizations of stratified flows at large Re.

Multiscale methods are often employed to study the effect of microstructure on macroscopic behaviour. For non-linear problems, these usually result in a two-scale formulation, where macro- and microstructure are simultaneously solved and coupled. If the microstructural features are much smaller compared to the macrostructural size, its effective behavior can be sufficiently predicted with first-order computational homogenization. However, if scale separation cannot be assumed or non-local effects due to buckling, softening, etc., emerge, higher-order methods, such as second-order homogenization [1], need to be considered. This formulation contains the second gradient of the displacement field, giving rise to a length-scale associated with the length-scale of the underlying unit cell, thus making it possible to capture size and non-local effects. Solving such problems is currently computationally expensive and typically infeasible for realistic applications, which limits the applicability of this method.

In this work, we address this issue by developing a reduced order model for second-order computational homogenization scheme based on Proper Orthogonal Decomposition. We consider different numerical examples and discuss different training strategies, computational savings and accuracy of the surrogate model.

Acknowledgements: This result is part of a project that has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 Research and Innovation Programme (Grant Agreement No. 818473).

[1] Kouznetsova, V., Geers, M.G.D. and Brekelmans, W.A.M. (2002), Multi-scale constitutive modelling of heterogeneous materials with a gradient-enhanced computational homogenization scheme. Int. J. Numer. Meth. Engng., 54: 1235-1260. https://doi.org/10.1002/nme.541

Manifold learning techniques such as Laplacian Eigenmaps (LE) [1] are commonly applied in fields like image and speech processing, to extract nonlinear trends from large sets of high-dimensional data [2]. Such techniques are also intriguing as model order reduction methods in multiscale solid mechanics: LE can capture nonlinearities in the solution manifolds of discretised physical problems [3]. Compared to POD-based algorithms, this may result in reduced order models yielding more accurate results with fewer parameters and lower computational effort [3]. Consequently, manifold learning techniques have been applied successfully to problems in fluid mechanics [3] and elastodynamics [4].

In the framework of the FE² method – in which computations on microscale representative volume elements (RVEs) are performed at each Gauss point of a macroscopic problem [5] – the payoff of such computational cost reduction may also be significant.

This contribution discusses the application of LE to model order reduction for RVE computations. Nonlinear, hyperelastic and elastoplastic behaviour is considered. The area of application comes with unique challenges and opportunities: for example, the mapping between reduced and original spaces and the projection of residuals onto reduced bases is not trivial [3]. On the other hand, the underlying PDEs [4] and the parametrisation of the RVE problem via a macroscopic deformation gradient and history variables [5] imply a strong (nonlinear) correlation in the unknown displacement degrees of freedom to be reduced. This talk will explore some of these challenges and opportunities.

[1] Belkin, Mikhail, and Partha Niyogi. "Laplacian eigenmaps for dimensionality reduction and data representation." Neural computation 15, no. 6 (2003): 1373-1396.

[2] Lee, John A., and Michel Verleysen. Nonlinear dimensionality reduction. Vol. 1. New York: Springer, 2007.

[3] Pyta, Lorenz Matthias. "Modellreduktion und optimale Regelung nichtlinearer Strömungsprozesse." PhD diss., Dissertation, RWTH Aachen University, 2018.

[4] Millán, Daniel, and Marino Arroyo. "Nonlinear manifold learning for model reduction in finite elastodynamics." Computer Methods in Applied Mechanics and Engineering 261 (2013): 118-131.

[5] Schröder, Jörg. "A numerical two-scale homogenization scheme: the FE 2-method." Plasticity and beyond: microstructures, crystal-plasticity and phase transitions (2014): 1-64.

Plasticity is the result of the motion and interaction of discrete dislocations in a crystalline material. Modelling plasticity at the crystal level based on discrete dislocation dynamics (DDD) is challenging due to the complexities associated with the dislocation activities of different slip planes. A data-driven approach provides an alternative method for simulating the complex behavior associated with plasticity at a small scale. We use methods based on dynamic mode decomposition1 (DMD) to analyze the DDD data2. We built reduced-order models for describing system dynamics with a few dominant modes. The models are built upon datasets of different physical resolution, e.g. dislocation density information resolved on the slip system level based on total dislocation densities, dislocation density vectors, or second-order dislocation alignment tensors. Different levels of spatial resolution are used to evaluate the effectiveness of models in reconstruction of the analysed data.

The modelling approach is then extended to forecast material response beyond the training dataset, for which we adopt more general (non-linear) Koopman operator theory and advanced stabilized DMD schemes, like shift invariant (physically informed) DMD or optimized DMD3. The different schemes are compared in their ability to predict the nonlinear behaviour in crystal plasticity from the DDD data.

REFERENCES

[1] Schmid, Peter J. "Dynamic mode decomposition of numerical and experimental data." Journal of fluid mechanics 656 (2010): 5-28.

[2] Akhondzadeh, Sh, Ryan B. Sills, Nicolas Bertin, and Wei Cai. "Dislocation density-based plasticity model from massive discrete dislocation dynamics database." Journal of the Mechanics and Physics of Solids 145 (2020): 104152.

[3] Askham, Travis, and J. Nathan Kutz. "Variable projection methods for an optimized dynamic mode decomposition." SIAM Journal on Applied Dynamical Systems 17, no. 1 (2018): 380-416.

The availability of high quality µ-CT images of materials allows for detailed multiscale simulation workflows in digital material characterization. In this case, data driven hybrid machine learning approaches are used to speed up full field approaches. Efficient and performance implementations of such data driven methods are essential for them being used for industrial applications.

This work concerns DMN (Deep Material Network) whose potential applications were exploited recently [1,2,3]. They only need linear elastic training data to identify equivalent laminate microstructure, which can be used to predict nonlinear behavior.

The industrial applicability of the DMN for short fiber reinforced plastic is investigated by comparing its speed and accuracy against direct numerical simulation results[4,5] on RVEs[6] using different physically nonlinear material behavior.

[1]- Liu, Z., Wu, C., & Koishi, M. (2019). A deep material network for multiscale topology learning and accelerated nonlinear modeling of heterogeneous materials. Computer Methods in Applied Mechanics and Engineering, 345, 1138–1168.

[2]- Liu, Z., & Wu, C. (2019). Exploring the 3D architectures of deep material network in data-driven multiscale mechanics. Journal of the Mechanics and Physics of Solids, 127, 20–46.

[3]- Gajek, S., Schneider, M., & Böhlke, T. (2020). On the micromechanics of deep material networks. Journal of the Mechanics and Physics of Solids, 142, 103984.

[4]- Matthias Kabel, Dennis Merkert, & Matti Schneider (2015). Use of composite voxels in FFT-based homogenization. Computer Methods in Applied Mechanics and Engineering, 294, 168-188.

[5]- Matthias Kabel, Andreas Fink, & Matti Schneider (2017). The composite voxel technique for inelastic problems. Computer Methods in Applied Mechanics and Engineering, 322, 396-418.

[6]- Schneider, M. (2022). An algorithm for generating microstructures of fiber-reinforced composites with long fibers. International Journal for Numerical Methods in Engineering, 123(24), 6197-6219.

Modeling for the description and prediction of processes in nature often leads to partial differential equations. Solving these field equations can only be done analytically in very few cases, so that in practice numerical approximation methods are often used. Variational methods like the Galerkin method have proven to be very effective and are widely used in industry and research. To set up the system of equations, integration over the area to be calculated is necessary. For more complex geometries or nonlinear equations, analytical integration becomes difficult or even infeasible, so that integration is also often performed numerically in the form of weighted evaluations of the integrand, the Gauss quadrature. In order to benefit from the quasi-optimal accuracy of the Galerkin method according to Cea’s lemma in the linear case, the quadrature scheme must also be of sufficient accuracy. On the contrary, for more complex constitutive laws, under-integration is often used in engineering to save computational time. Based on a split of geometric and material nonlinearities, the present talk introduces a one-point integration scheme that is able to integrate polynomial shape functions of arbitrary order geometrically accurate. The material nonlinearity can be captured with the desired accuracy via a Taylor series expansion from the nonlinear state. As a demonstration the integration scheme is applied to two-dimensional polygonal shaped second order virtual elements where the quadratic projection is integrated via a single integration point.

The development of effective and locking free shell elements is intensive topic of research since several decades. Recently, the Hellan-Herrmann-Johnson (HHJ) method for linear Kirchhoff-Love plates has been extended to nonlinear Koiter shells. Therein, the bending moment tensor is introduced as additional unknown to rewrite the fourth order as a second order mixed saddle point problem circumventing the necessity of C1-conforming finite elements. Via hybridization techniques the saddle point translates into a minimization problem again.

The tangential-displacement and normal-normal-stress continuous (TDNNS) method has successfully been applied to linear Reissner-Mindlin plates leading to a shear locking free formulation.

In this talk we present a shear locking free extension of the TDNNS method from linear Reissner-Mindlin plates to nonlinear Naghdi shells by means of a hierarchical approach. Therefore, the HHJ method for Koiter shells is enriched with shearing degrees of freedom, discretized by H(curl)-conforming Nedelec elements. We discuss the small-strain regime leading to the HHJ and TDNNS method for linear Koiter and Naghdi shells. We show how the so-called Regge interpolant can be used in all methods to avoid membrane locking by inserting into the membrane energy term.

Several benchmark examples, implemented in the open-source finite element software NGSolve (www.ngsolve.org), are presented to demonstrate the excellent performance of the proposed shell elements.

We model the formation of wrinkles of an elastic substrate coated with a thin film. The elastic substrate is first stretched, then the film is attached to a part of the substrate boundary in the deformed state. Once the the external force is released, wrinkles form due to the stress mismatch between the two materials. The elastic substrate is modeled using a hyperelastic, homogeneous and isotropic material. The film is modeled using a geometrically exact Cosserat shell. The resulting deformation and microrotation $(varphi, R)$ are a minimizing pair of the combined energy functional

$$

J(varphi, R) = int_{Omega} W_textup{bulk}(nablavarphi) : dV + int_{Gamma_c}W_textup{coss}(nablavarphi_{|_{Gamma_c}}, R) : dS

$$

in the admissible set

begin{align*}

mathcal{A} = Big{&(varphi, R) in W^{1,q}(Omega, mathbb{R}^3) times H^1(Gamma_c,textup{SO(3)}) : Big| :

varphi textnormal{ is a deformation function, }& (varphi, R) textnormal{ fullfill the Dirichlet boundary conditions} Big}

end{align*}

with $q > 3$.

We discretize the problem using Lagrange finite elements for the substrate displacement. For the numerical treatment of the microrotation field, standard Lagrange finite elements cannot be used, as the microrotation field maps to the nonlinear manifold $textnormal{SO}(3)$. We present a generalization of Lagrange finite elements that is suitable for such manifold-valued functions: geometric finite elements.

The resulting finite element spaces are complete and invariant under isometries of the manifold. The best approximation error depends on the mesh size h. We prove the existence of solutions of the discrete coupled model. We compare two Newton-type methods to solve the resulting discrete problem: a Riemannian trust-region method and a Riemannian proximal Newton method.

Numerical experiments show that we can efficiently reproduce wrinkling patterns of coupled systems. Our approach works as well for more complex scenarios like multi-layer systems or systems involving various stress-free configurations.

In many different fields of engineering beam models play a significant role in the efficient simulation of slender structures. The most important model for large deformations is the so-called geometrically exact beam also often referred to as Simo-Reissner beam. The configuration manifold of the beam model is given by special Euclidian group as it describes the position of the centerline as well as the orientation of the beam's cross-section. The partial differential equations describing the behavior of the beam is usually solved with the help of the Finite Element Method (FEM). So it becomes necessary to discretize the special orthogonal group in a finite element sense.

A finite element discretization of the special orthogonal group is rather difficult as the special orthogonal group is not an abelian, additive group but a matrix group under multiplication. Though there exist parametrizations of the orthogonal group, which have an additive structure, they result in path-dependency. This can be overcome by discretizing the group directly by using so-called directors. The directors can be discretized additively, so in a classical finite element sense. This, however, leads to an increase in the number of degrees of freedom

if Lagrange multipliers are used to ensure the orthonormality of the directors. Further, this formulation does not conserve the structure of the manifold at every point of the discretization. A possible remedy could be a projection method via the polar decomposition, which is very costly in numerical terms.

The use of unit quaternions for the parametrization presents an interesting alternative. Even though unit quaternions have a complex mathematical structure, it can easily be ensured that their unit length is conserved after a finite element discretization by normalizing the discretized quaternions. This still allows for a classical additive discretization technique in a finite element sense.

In the literature, it is often shown that the Isogeometric Analysis (IGA) is advantageous over the classical FEM with Lagrangian elements, especially for dynamic problems. We thus apply the IGA to the quaternion formulation of the geometrically exact beam.

Thermal convection is a phenomenon seen in almost all facets of life, ranging from planetary convection to ocean currents and convection inside the earth. The physics of thermal convection complicates when a porous wall layer is present. Flow over urban canopies, forest canopies or flow in underground aquifers are classic examples of such scenarios where thermal buoyancy-driven convection occurs in the presence of turbulent flow over a porous wall layer. The present research work focuses on simulating pressure-driven turbulent flow over a simplified, ordered porous medium consisting of a regular array of cubes. The work further couples it with natural convection arising due to unstable stratification, to provide insight into the momentum and heat transfer characteristics of such a flow scenario.

Direct numerical simulations (DNS) have been performed with a finite-difference solver to validate the model for buoyancy-driven convection and the classical Rayleigh-B´enard convection. Further, we extended the solver with an Immersed Boundary Method (IBM) to model the ordered porous medium, which was validated against reference data. The focal point of the present research, analyzing mixed convection over a porous wall layer, brings into the picture a large number of dimensionless control parameters. The bulk Reynolds in the overlying free channel region is fixed at 5500, the Prandtl number at 0.71. We impose an adiabatic boundary condition on the surface of the cubes. We varied the flux Richardson number to cover different flow scenarios between pure shear and purely buoyancy-driven flows. The porous and free regions are expected to show different convective patterns and different critical flux Richardson numbers for the transition to natural convection cells. Further, the interface regime dynamics should provide insight into the heat transfer characteristics, since the heat transfer timescales vary drastically between the porous region and the turbulent flow region.

We use direct numerical simulation (DNS) to study the problem of drag reduction in a lubricated channel, a flow instance in which two thin layers of a lubricating fluid (e.g. water) are injected in the near-wall region of a plane channel, so to favor the transportation of a primary fluid (e.g. oil). All DNS are run within the constant power input (CPI) approach, which prescribes that the flow-rate is adjusted according to actual pressure gradient so to keep constant the power injected into the flow. A phase-field method (PFM) is used to describe the dynamics of the liquid-liquid interface and when prescribed, also the presence of surfactants/contaminants. As this technique is tailored toward the transport of very viscous fluids like oils, we study the drag reduction performance of the system by keeping fixed the lubricating fluid properties (water) and by considering two different types of oil characterized by different viscosities, 10 and 100 times larger than that of water, respectively. As these systems are also characterized by the presence of contaminants and surfactants – which act by locally reducing the local value of the surface tension – for each type of transported oil, we consider a clean and a surfactant-laden interface. For all the four tested configurations, we unambiguously show that a significant drag reduction (DR) can be achieved. Upon a detailed analysis of the turbulence activity in the two lubricating layers, the interfacial wave dynamics and their interplay, we are able to characterize the effects of surface tension forces, surfactant concentration and viscosity contrast on the drag reduction performance.

We use pseudo-spectral Direct Numerical Simulation (DNS), coupled with a Phase Field Method (PFM), to investigate the turbulent Poiseuille flow of two immiscible liquid layers inside a channel. The two liquid layers, which have the same thickness (h1 = h2 = h), are characterised by the same density (ρ1 = ρ2 = ρ) but different viscosities (η1≠ η2), so mimicking a stratified oil-water flow. This setting allows for the interplay between inertial, viscous and surface tension forces to be studied in the absence of gravity. We focus on the role of turbulence in initially deforming the interface and on the subsequent growth of capillary waves. Capillary wave propagation and interaction is studied by means of a spatiotemporal spectral analysis and compared with previous theoretical and experimental results. Wave propagation is found in agreement with the theoretical dispersion relation. At wave scales larger than the turbulent forcing range the observed scaling of the one-dimensional wavenumber spectrum suggests an energy equipartition regime, which is predicted by theory and recently has been observed in experiments with capillary wave turbulence in microgravity. At wave scales directly forced by hydrodynamic turbulence an initially mild slope of the wavenumber spectrum is succeeded by a sharp decay of wave energy at larger wavenumbers, with the transition taking place near the Kolmogorov-Hinze critical scale, where surface tension forces and turbulent inertial forces are balanced.

The motion of dislocations has been determined to be one of the main mechanisms leading to inelastic deformation in crystalline materials. Their motion is affected by other crystal imperfections, e.g., at grain boundaries their advancement is hindered due to misalignment between the crystals' slip systems. The pile-up that occurs at the boundaries can lead to yielding inside the adjacent grains or intergranular fracture. The damage induced by the latter acts as a precursor to failure at the macroscopic scale. As such, a formulation capable of describing the interaction between the aforementioned crystal imperfections could provide a feasible tool to predict failure of components made from crystalline materials.

To this end, a gradient crystal plasticity formulation which accounts for grain misorientation is enhanced by considering the grain boundary as a cohesive interface and by introducing a damage variable influencing the interaction between adjacent grains. Numerical examples demonstrating the material response based on the proposed formulation are presented and discussed.

Heat treatment plays an essential role in the production of cold-work steel parts. While the material properties are adjusted by the heat treatment, side effects like distortion and residual stresses have to be controlled. A good prediction of the heat treatment plays a major role in reducing the necessary grinding time in subsequent finishing operations. Optimising the heat treatment process has the potential to save energy in the furnaces. This presentation discusses the application of simplified material models for the finite element (FE) simulation of quenching. The formation of martensite is covered by a purely temperature dependent Koistinen-Marburger model, whereas the diffusive formation of Bainite is modelled with an incrementally isothermal Johnson-Mehl-Avrami-Kolmogorov relation [1]. Both models are used in rate format and solved monolithically. The thermal-mechanical-microstructural-coupling implemented in the FE-software Abaqus is presented alongside numerical examples.

[1] de Oliveira, W.P., Savi, M.A., Pacheco, P.M.C.L., 2013. Finite element method applied to the quenching of steel cylinders using a multi-phase constitutive model. Arch Appl Mech 83, 1013–1037. https://doi.org/10.1007/s00419-013-0733-x

Modeling single crystal plasticity is essential for understanding the behavior of polycrystalline materials such as metals and alloys. The mechanical properties of such materials depend on the microstructure of individual grains and their interaction through grain boundaries. Single crystal plasticity aims to model the behavior of an individual grain based on the microscopic lattice structure. It can be expressed mathematically using the concept of multisurface plasticity. Applying the principle of maximum plastic dissipation leads to an optimization problem where the individual slip systems of the crystal, represented by yield criteria, define the constraints of the optimization problem.

In the framework of rate-independent crystal plasticity models, the set of active slip systems is possibly non-unique, which makes the algorithmic treatment challenging. Typical approaches are either based on an active set search using various regularization techniques [3] or simplifying the problem in such a way that it becomes unique [1]. In computationally extensive simulations, the problem needs to be evaluated multiple times. Therefore, a stable, robust, and efficient algorithm is required to obtain satisfactory results.

Recently, an alternative strategy based on the infeasible primal-dual interior point method (IPDIPM [2] has been presented in [4], which handles the ill-posed problem without perturbation techniques. Through the introduction of slack variables, a stabilization of the conventional active set search approach is reached. The introduction of barrier terms with related barrier parameters continuously penalizes the violation of the feasibility of the intermediate solution. This talk especially focuses on the treatment of the barrier parameter and the related speed of convergence.

[1] M. Arminjon. A Regular Form of the Schmid Law. Application to the Ambiguity Problem. Textures and Microstructures, 14:1121–1128, 1991.

[2] A. S. El-Bakry, R. A. Tapia, T. Tsuchiya, and Y. Zhang. Journal of Optimization Theory and Applications, 89(3):507–541, 1996.

[3] C. Miehe and J. Schr ̈oder. A comparative study of stress update algorithms for rate-independent and rate-dependent crystal plasticity. International Journal for Numerical Methods in Engineering, 50:273–298, 2001.

[4] L. Scheunemann, P. Nigro, J. Schröder, and P. Pimenta. A novel algorithm for rate independent small strain crystal plasticity based on the infeasible primal-dual interior point method. International Journal of Plasticity, 124:1–19, 2020.

Metal additive manufacturing (MAM) has received significant attention in recent years due to its significant advantages such as increased design flexibility for complex geometries, shorter production-cycle, and efficient use of raw materials. To fully realize the potential of MAM in the context of Industry 4.0, it is necessary to address challenges related to the mechanical reliability of printed parts and their associated costs. Currently, trial-and-error methods are the most common way of optimizing MAM process conditions for achieving the desired printing quality. Meanwhile, numerical simulations can provide a more profound understanding of the physical phenomena involved in the build process, leading to a more systematic optimization of process conditions, and ultimately making the `first-time-right' high-quality production possible. Achieving a thorough quantitative understanding of the process requires insights from models covering various physical aspects including thermal, mechanical, metallurgical, and fluid-dynamics interactions. However, high-fidelity simulations of such models are accompanied by significant computational costs and therefore have limits in applications, particularly in sensitivity and optimization analyses where solutions for a wide range of scenarios are required.

To address this challenge, we initiated a project in 2021, with the support of the Swiss National Science Foundation (SNSF), to explore the feasibility of meaningful acceleration of these simulations without significant compromise in accuracy and reliability. Specifically, the project aims to develop solutions for thermal, microstructure, and residual stress simulations for the laser powder bed fusion (LPBF) process. To generate experimental validation data, Hastelloy X serve as the 'model material'. An overview of the results obtained so far, focusing on thermal and microstructure simulations, are presented.

Several techniques have been examined to reduce the computational cost of thermal simulation for LPBF, including a multi-scale simulation strategy, surrogate modelling, and physics-informed neural networks (PINNs), where the advantages and limitations of each approach are discussed. In the field of microstructure modelling, a 'Neural Cellular Automata' method has been developed, which outperforms the conventional Cellular Automata with up to 6 orders of magnitude acceleration in computation speed. Moving forward, the project will continue with a focus on the development of affordable high-fidelity models of residual stress development until 2025.

Additive Manufacturing (AM) and especially its metal variants constitute today a reality for the fabrication of high-performance industrial components. In particular, AM allows the construction of novel cellular structures, the so-called lattices, where well-designed unit cells are periodically repeated over a macro-shape to achieve exceptional specific performances, such as unprecedent stiffness-to-weight ratios. These structures, however, are very difficult to simulate numerically: on the one hand, the application of multiscale methods based on homogenization appears delicate due to an insufficient separation of scales (macro versus cell scales); on the other hand, solving directly the high-fidelity, fine-scale problem requires handling large numbers of complex cells which is often intractable if standard methods are blindly used. As a solution, immersed domain techniques have been applied, but such methods, generic in terms of applications, may not be optimal in the case of lattices.

In this context, the purpose of this work is to develop a HPC algorithm dedicated to lattices that takes advantage of the geometric proximity of the different cells in the numerical solution. In order to do so, we start by adopting the CAD paradigm based on spline composition along with its corresponding IGA framework. This offers (i) great flexibility to design any lattice geometry and (ii) fast multiscale assembly of the IGA system. Then, we resort to the family of Domain Decomposition solvers, and develop an inexact FETI based algorithm that avoids solving numerous local cell-wise systems. More precisely, we extract the “principal” local cell stiffnesses using a greedy approach, and use the latter as a reduced basis to efficiently solve all the cell-wise systems. It results in a scalable algorithm that tends to be matrix-free. During the talk, a range of numerical examples in 2D and 3D will be presented to account for the efficiency of our method both in terms of memory and computational cost reduction.

Discontinuous fiber reinforced composites are used in many application areas ranging from automotive to healthcare. Such parts are often manufactured in and injection molding process, as it is an economical process for high volume markets. The simulation of the injection molding process is well established and specific commercial tools have been developed for this task. However, the transient solution of the underlying non-linear multi-phase flow is computationally expensive and computation may take multiple hours for complex geometries. This computational time is prohibitively large for computational optimization of the product design or the process parameters. Hence, we propose a two-step process to accelerate the mold filling prediction: i) Solve a modified Eikonal equation to compute distance maps to the injection gate and nearest walls. This is computationally cheap, as it is only a stationary equation to solve. ii) Train feed forward neural networks to obtain a data-driven relation between the encoded distance maps and mold filling features, such as fill time and fiber orientation. We sample a set of geometries, automatically generate CAD models, and simulate these in a commercial injection molding solver to build a training data set. Subsequently, we apply different feed forward neural network architectures and evaluate their performance.