To capture and predict the chemo-thermo-mechanical behavior of ceramic foam filters, material models and simulation tools are required. The description of the thermo-mechanical inelastic behavior as well as the in-situ layer formation on reactive filters have been the aims of this subproject. Challenging aspects in the whole progress are the exact geometrical replication of the underlying foam structure of the filter and the lack of experimental data for many relevant loading cases. The software FoamGUI is developed to generate parametrized, periodic three-dimensional representative volume elements (RVE) of foam structures, which are used in continuum and fluid mechanical simulations as well as for 3D-printing. Calculation concepts are formulated to predict the inelastic deformation and failure behavior of ceramic open-cell foams under thermo-mechanical loading. First-order homogenization approaches are used to conclude from the mesoscopic behavior of the foam RVE to the macroscopic response of filter structures. A hybrid approach is developed in the established framework of rate-independent plasticity in combination with neural networks, which replace the plastic flow potential and the evolution equations of internal state variables. Another modeling aspect is motivated by the experimentally observed growth of an in situ layer during the so-called reactive phase of the filtration process. This phenomenon motivates the development of a model to describe diffusion, chemical reactions and phase transition processes of multi-phase/multi-component systems using the phase-field method. This allows the simulation of spatially and temporally resolved microstructure evolution leading to the layer formation.
16.1 Introduction
Metal melt filtration is a technological process used to clean the melt and to calm the melt flow during casting. The innovative concept of the multifunctional filters studied here centers on complementing the physical removal of impurities (particles) by additional cleaning mechanisms enabled by chemical reactions between melt and filter material as well as active coatings.
Figure 16.1 illustrates the principle of metal melt filtration and also shows the micrograph of an actual filter cross-section. To design and produce multifunctional filters that are able to withstand the very demanding thermo-mechanical conditions to be encountered during such filtration processes is a major challenge. These structures are, for instance, exposed to large temperature rates and gradients, when the hot molten metal hits the filter (thermal shock), which in turn result in thermally-induced stresses that can cause damage or, in extreme cases, complete failure. Additional stresses are caused by the melt flow, since the filters cause a pressure drop in the direction of flow. Even buoyancy related forces can play a role when ceramic filters are immersed in metal melt, due to the significant density differences of these materials. It must also be emphasized that instead of failing in a brittle manner, ceramic filter materials show non-negligible inelastic behavior (plasticity, creep, damage, etc.) at temperature levels (of up to 1.600 \(^\circ \)C) required for steel filtration.
A set of 2 schematics. The left schematic illustrates a central focus of liquid metal being poured onto a filter. The right schematic illustrates the components of a ceramic filter including pores, inclusions, micro-cracks, cavity, and interface cracks.
Fig. 16.1
Schematic of liquid metal filtration [1] and a detail of a ceramic filter with a coating
×
Anzeige
The foams studied in this research are reticulated foams, which are based on a polyurethane foam coated by a ceramic slurry. The coating is later dried and fired, whereas the polyurethane pyrolizes and sharp edged cavities remain in the strut [2]. The sharp edges are potential locations for the initiation of crack growth. The firing process also influences the thermo-mechanical properties of the foam. Additional (active) coatings can be applied, which allow specific chemical reactions between melt and filter surface to take place. These reactions can help to reduce certain unwanted elements in the melt.
In order to be able to provide simulation support for the filter design, a thermo-mechanical multiscale modeling approach has been pursued, based on state-of-the-art methods of theoretical and computational mechanics. The goal here was to enable virtual evaluations of the structural integrity and strength of the filters. These are of critical importance, since experimental investigations under realistic process conditions are often very difficult, expensive or even impossible to conduct. To this end, the modeling efforts have focused on three main aspects:
1.
The generation of open-cell foam-like Representative Volume Elements (RVEs), with direct control over geometrical features of the foam morphology
2.
Multiscale Finite Element Analysis (FEA) of the foam behavior that combines models relevant for different length-scales with scale-bridging concepts (homogenization)
3.
More recently, phase-field modeling of an in situ layer formation during the chemical reaction of steel melts and carbon-bonded aluminum oxide filters (reactive filtering phase).
Finite element analysis is generally suitable to model the structural behavior of the foam-like filters on all relevant length-scales. However, such filters consist of a large number of three-dimensionally interconnected struts. Even for macroscopically relatively small filters, the element numbers required for sufficiently accurate meshing under full spatial resolution would thus reach orders of magnitude that could no longer be reasonably handled, even with modern high-performance computers. Homogenization concepts must therefore be applied, if the goal is to predict macroscale properties that are directly influenced by micro- and mesoscale features and behavior, see Fig. 16.2. To this end, the proposed approach is to model the effective thermo-mechanical behavior of the filter material on the microscopic scale (within the struts) via elastoplastic and viscoplastic constitutive models, combined with simple damage approaches. These are established models that are available in many commercial FE programs. On the mesoscale, a representative volume element of the foam structure is simulated, in which the detailed geometrical structure is spatially discretized with finite elements. The RVE is then loaded via appropriate boundary conditions and the numerical analysis results in a model for the effective thermo-mechanical behavior of the foam-like mesoscale structure that can be evaluated at every integration point of a macroscale FE simulation.
A set of 3 micrographs. The left scan represents the micro-scale. The center scan represents the meso scale. The right scan represents the macro scale. From the left to right, homogenization takes place, indicated by an arrow.
Fig. 16.2
Relevant scales for the three-scale homogenization approach
×
In general, the effective inelastic behavior of such highly porous structures which the RVEs represent is very complex and may exhibit volumetric strain hardening or softening, and anisotropic non-associated plastic flow behavior, just to name a few phenomena typically encountered. Modeling approaches that try to formulate appropriate flow potentials and flow rules in phenomenological plasticity models are therefore limited in their applicability and accuracy. Classical models include the approach of Deshpande and Fleck [3] and the Cam-Clay model of Roscoe [4], and generalizations suggested in Bigoni-Piccolroaz [5] and Ehlers [6, 7]. There also exists a large body of work regarding the Theory of Porous Media (TPM) that is very relevant to this problem, see, e.g., [8] and the references therein.
Anzeige
The strategy to link the meso- and macro-scale behavior of the ceramic filters followed here, however, is different in nature and can be classified in two categories. Firstly, a novel hybrid approach has been developed, in which the flow potentials and flow rules are replaced by neural networks, while the classical structure of phenomenological elasto-visco-plasticity models known from continuum mechanics is retained. Secondly, direct numerical homogenization via the FE2 method has been conducted. In this approach, the explicit formulation of a macroscopic constitutive model is typically not necessary. Macroscopic stress-strain relations are instead established by simulating the response of the RVE at each integration point of the macroscale FE problem. Nested, but otherwise separate FE simulation are thus performed on two scales, hence the name FE2. However, while greatly reducing the degrees of freedom and memory required for the simulation compared to the full spatial resolution for a foam structure of technologically relevant size, this approach is computational still very extensive. Recently, a lot of research has therefore focused on significantly improving the computational efficiency of this method, see [9, 10].
The generated RVEs are, however, not only used for thermo-mechanical FE simulations. The underlying digital geometric models can be exported in different formats. Voxel structures are provided for simulations based on the Lattice-Boltzmann-Method (LBM), where the flow of the liquid melt through the filter can be analyzed [11, 12]. The geometric surface of the generic structures can be triangulated and exported as standard triangle language (STL) data structures, which are commonly used for 3D printing. The printed structures can be analyzed in various experiments to measure actual solid- or fluid-mechanical properties. Very importantly, the models can also be regarded as digital twins of real foam structures, which have the advantage that geometrical and topological properties can be easily varied and studied either using numerical or real world experiments.
Finally, recent continuum thermodynamics-based modeling approaches have been devoted to investigating the role of thermo-chemical phenomena as key mechanisms enabling the functionalization of carbon-bonded alumina filters for the purification of steel melt. In particular, the focus is on the formation of an in situ layer of secondary oxides during filtration, which has a high impact on the cleaning efficiency. During this reactive filtration phase, a dissolution of alumina from the filter material is assumed in the presence of carbon and liquid iron [13], where the produced oxygen and carbon can diffuse within the melt and nucleate carbon monoxide gas bubbles at impurities yielding a cleaning effect based on flotation [14, 15]. The formation of a dense layer of secondary corundum is proposed as the main reason for the termination of the reactive filtration phase [16], as it is suspected to suppress further carbon supply from the carbon-rich filter to the melt [17]. To capture this, a multi-component/multi-phase phase-field modeling approach is pursued to describe the diffusion, phase transitions, and chemical reactions in the system consisting of filter and melt material. In a first study following this concept, simulation results for the in situ layer formation are presented assuming a simplified filter-melt system.
16.2 Methods
The focus of this section is on the numerical methods that have been established to generate artificial foam structures, to model the effective thermo-mechanical behavior of the foam filters, and to analyze the formation of in situ layers during steel melt filtration with carbon-bonded alumina filters.
16.2.1 Representative Volume Elements of Foam Structures
There are several ways to generate RVEs of foam structures. One could use data from computed tomography to map realistic structures. However, in order to reproduce specific geometrical features and to study their influence on the macroscopic behavior, a generic approach is beneficial [1]. The process used here starts with generating a spatially periodic sphere packing. The spheres may be arranged in a lattice structure as shown in Fig. 16.3, or randomly, using algorithms proposed in the literature [18‐20]. From the sphere packing, a Laguerre tessellation is generated, which contains convex cells around each sphere. But these cells usually exhibit some very small facets, which is not observed in real foams. A relaxation procedure is therefore applied in order to minimize the total foam surface by topological changes and thereby the surface energy [21]. The edges of the relaxed foam structure are mapped into a 3D voxel grid, where each voxel contains the length of struts within itself as a numerical value. This strut density field is smoothed out by a Gaussian filter and finally the surface of the foam structure can be represented by an isocontour (level set). By changing the value of the standard deviation of the Gaussian filter and/or the value for the isocontour, the shape of the struts and the relative density of the generic foam can be adjusted [1].
A set of 4 3-D illustrations of foam structures with increasing pores from left to right. The first illustration has 2 pores, the second illustration has 8 pores, the third illustration has 16 pores, and the fourth illustration has 216 pores.
Fig. 16.3
Generic foam structure RVEs with a relative density of \(\rho _{\text {rel}}=20\%\) and 2, 8, 16 and 216 pores. The foam with two pores is also known as a Kelvin foam and the one with eight pores as the Wheire-Phelan foam [22]
×
Subsequently, the generated RVEs can be discretized in different ways with finite elements to perform numerical studies. To analyze mechanical properties such as elasticity, plasticity or creep behavior, the struts are meshed using hexahedral elements as was shown in [1, 23], where sharp edge cavities had to be taken into account for analyzing the fracture mechanical behavior. The same geometries are also used for fluid mechanical simulations, utilizing fine meshed voxel grids (see Fig. 16.4) and assigning either solid or fluid properties to the individual elements [11]. Furthermore, the closed surfaces of the individual RVEs can be exported as STL-files, which may be used to feed 3D printers [24, 25] or being analyzed by software tools like MeshLab to estimate the surface area or topological properties such as the number of cell windows [11].
16.2.2 Hybrid Homogenization Approach
In our notation, homogenized quantities are denoted by an overbar. Homogenized stresses \(\bar{\boldsymbol{\sigma }}\) and strains \(\bar{\boldsymbol{\varepsilon }}\) are then introduced as the volume-averaged microscopic stress \(\boldsymbol{\sigma }\) and strain \(\boldsymbol{\varepsilon }\) field, i.e.,
All foam RVEs considered in this research, see the examples in Fig. 16.3, are periodic in all three spatial dimensions. For the mechanical analyses periodic boundary conditions are applied, since this gives the optimum response [26]. The local displacements \(\boldsymbol{u}\) at a position \(\boldsymbol{x}\) on the RVE boundary \(\partial V\) depend on the macroscopic strain \(\bar{\boldsymbol{\varepsilon }}\) and a fluctuation \(\tilde{\boldsymbol{u}}\)
at homologous points \(\boldsymbol{x}^-\) and \(\boldsymbol{x}^+\) at the opposing boundaries of the RVE. The fluctuations can be eliminated from (16.2) yielding
which is implemented as multi point constraints for all sets of homologous points at the RVE boundary. In case of non-periodic meshes a technique developed by Storm et al. [27] may be applied.
Elastic Properties and Failure Limit Surfaces for Foam RVEs
To predict the effective elastic properties and failure limit surfaces local stress fields are determined for six linearly independent strain controlled load cases, denoted by an upper index k in parentheses
\(\bar{\;\!\varepsilon }_i^{(k)}\) denotes the effective strain state for load case k in Voigt notation, whereas \(\varepsilon _0\) is a small scalar strain value for which no inelastic response of the RVE is expected. Since \(\bar{\varepsilon }_i^{\;\!(k)}\) has only a single non-zero component \(\delta _i^{(k)}\), we can determine the \(k^\text {th}\) column of the effective stiffness tensor in Voigt notation using
For a failure limit surface the local stress fields \(\boldsymbol{\sigma }^{(k)}(\boldsymbol{x})\) for each load case k are computed. Here, \(\boldsymbol{x}\) denotes the position within the model. Then, for an arbitrary effective stress state \(\bar{\boldsymbol{\sigma }}\) the corresponding strain state is \(\bar{\boldsymbol{\varepsilon }}=\bar{\mathbb {C}}^{-1}: \bar{\boldsymbol{\sigma }}\). The resulting local stress field can be computed using the superposition principle
Different failure criteria may be applied, depending on the material. Just to name a few, there are the maximum principal stress criterion, where \(\sigma _\text {f}(\bar{\boldsymbol{\sigma }})=\max \left[ \sigma _1(\boldsymbol{x})\right] \ge \sigma _\text {c}\), or the maximum equivalent stress criterion with \(\sigma _\text {f}(\bar{\boldsymbol{\sigma }})=\max \left[ \sigma _{\text {eq}}(\boldsymbol{x})\right] \ge \sigma _\text {c}\), and the Weibull stress criterion for brittle ceramics, where \(\sigma _\text {f}(\bar{\boldsymbol{\sigma }})=\sigma _{\text {W}}\left( \boldsymbol{\sigma }(\boldsymbol{x})\right) \ge \sigma _\text {c}\). The equivalent stress is defined as
where the microscopic loading state depends on the macroscopic stresses \(\sigma _\text {f}(\bar{\boldsymbol{\sigma }})\). To obtain their shape for a specific failure criterion on the microscale, these functions need to be evaluated for many different macroscopic stress states. The following approach takes samples in stress space in a systematical way. A direction in stress space is defined by
with \(\alpha \in \left[ -\frac{1}{2}\pi ,\frac{1}{2}\pi \right] \) denoting the angle between the stress direction and \(\pi \)-plane. The deviatoric part is defined as
where \(\boldsymbol{n}\) denotes the eigendirections of the stress tensor, \(\theta \in \left[ -\frac{\pi }{6},\frac{\pi }{6} \right] \) the Lode angle and \(\boldsymbol{I}\) the second order identity tensor. The angles \(\alpha \) and \(\theta \) can be obtained from the corresponding stress tensor \(\bar{\boldsymbol{\sigma }}\) and its invariants
This concept of defining a direction \(\boldsymbol{N}\) can be used to any second order tensor and will be applied for different purposes in this chapter. Equation (16.11) a priori satisfies \(\Vert \boldsymbol{N}\Vert =1\) and a critical effective stress can be defined \(\bar{\boldsymbol{\sigma }}_\text {c} = \lambda _c \boldsymbol{N}\) with \(\lambda _\text {c}\) expressed as
Similar approaches have been applied by Zhang et al. [28] and Storm et al. [29] for different realizations of Kelvin foams.
A 3 D structural illustration of a foam with 216 pores. The stress ranges from 3.921 e minus 07 to 9.829 e minus 04, as indicated by a color gradient scale on the left.
Fig. 16.4
Stress distribution in a F216 foam RVE loaded by an effective shear load case. Here, the FE mesh is made up by cubic voxels for a specific strut geometry and relative density of \(\rho _{\text {rel}}=10\%\)
×
A suitable failure criterion for a fracture mechanical analysis is the coplanar equivalent stress intensity factor (SIF)
It takes the local stress-state at the crack tip into account caused by superposition of mixed mode loading \(K_I\), \(K_{II}\) and \(K_{III}\) [30]. For a homogeneous isotropic bulk material the three local SIFs \(K_i\) can be related to the J-integral using
Therein, E and \(\nu \) denote Young’s modulus and Poisson’s ratio, respectively. Furthermore, \(Y_{ij}\) is the Irwin matrix for a homogeneous isotropic material
The interaction integral \(J_j^{\text {int}}\) and the inverted Irwin matrix together with the unit value SIF \(k_0\) are used to separate the SIFs \(K_i\)
as the superposition of the actual fields of stress \(\sigma _{kl}\), strain \(\varepsilon _{kl}\) and displacement gradient \(u_{k,i}\) with their corresponding auxiliary fields \(\sigma _{kl}^m\), \(\varepsilon _{kl}^m\) and \(u_{k,i}^m\), respectively. The auxiliary fields are the a priori known near crack tip solutions for pure mode m loading causing a unit value SIF \(k_0\). To evaluate the interaction integrals a finite element model of a Kelvin foam as shown in Fig. 16.5 is utilized. Additional sub-models in form of tube like crack models (white mesh in Fig. 16.5 right) are used to compute the interaction integrals along the corresponding crack fronts. The sub-models get as boundary conditions the interpolated displacements of the global Kelvin cell model. More details about this particular modeling approach are given in Settgast et al. [30].
A set of 2 3 D models. The left model represents the Kelvin foam with multiple hexagonal rings interconnected in a complex arrangement. The right model is a wireframe 3 D model with a narrowing neck, a wider base, and a central trigonal cavity. Numerous lines define its contours.
Fig. 16.5
Finite element model of the Kelvin foam including the sharp edged strut cavities and a detail of the mesh (black grid) together with the submodels (white grid) for the calculations of stress intensity factors [30]
×
Hybrid Constitutive Model Incorporating Neural Networks
In this section a hybrid constitutive model for a macroscopic foam structure is derived. For the sake of simplicity the described model is restricted to elastic-plastic material behavior in a small strain setting. The strain rate tensor is additively split into an elastic and a plastic part
where the components of the effective stiffness tensor \(\bar{\mathbb {C}}\) can be determined by at most 6 linear independent load cases as described in detail in [31]. A yield function on the macroscale \(\bar{\phi }\) is defined, where the norm of the effective stress is compared with a value, which is predicted by a neural network\(\text {NN}^{\bar{\phi }}\)
The input values for the network are two angles describing the stress direction \(\boldsymbol{N}=\bar{\boldsymbol{\sigma }} / \Vert \bar{\boldsymbol{\sigma }} \Vert \), where \(\alpha \) denotes the angle between stress direction and \(\pi \)-plane and \(\theta \) the Lode angle, see (16.14). In other words for a given stress direction and accumulated plastic strain \(\bar{\varepsilon }_{\text {q}}^{\text {pl}}=\int _0^t \Vert \dot{\bar{\boldsymbol{\varepsilon }}}^{\text {pl}}\Vert \,\text {d}t\)internal variable, the neural network \(\text {NN}^{\bar{\phi }}\) predicts the stress amplitude \(\Vert \bar{\boldsymbol{\sigma }}\Vert \) required to achieve plastic flow.
In general, foams often show non-associated plastic flow, which would require a second dissipation potential for the constitutive description, from which the plastic flow direction is derived. Here, it is assumed that the eigendirections of stress and plastic flow are coaxial, which holds for isotropic foams. The plastic strain rate is defined as
and the direction of plastic flow is given using (16.11) via two different angles \(\alpha ^{\varepsilon }\) and \(\theta ^{\varepsilon }\), which are predicted by a second neural network
complete the model. Details of the implementation of such a model into FE codes are given in Malik et al. [32] and Settgast et al. [33].
Generation of Training Data and Neural Network Training
The process for the generation of sample data for the neural network training is independent of the choice of the RVE. Using values for \(\alpha ~\in ~[-\pi /2,\pi /2]\), \(\theta ~\in ~[-\pi /6,\pi /6]\) and \(\lambda ~\in ~[0,1]\) a unit direction in stress or strain space is defined according to Eqs. (16.11) and (16.12). For isotropic foams the eigendirections \(\boldsymbol{n}_k\) coincide with the three base vectors. The effective loads for the RVE are defined as
where \(\sigma _{\max }\) and \(\varepsilon _{\max }\) denote a maximum amplitude either in stress or strain space. For each loading direction a number of \(n_{\lambda }\) effective stress and strain tensors are computed using the FEM. The norm of the effective stress tensor is determined as target for the training of \(\text {NN}^{\bar{\phi }}\). The effective plastic strain is determined using \(\bar{\boldsymbol{\varepsilon }}^{\text {pl}}=\overline{\boldsymbol{\varepsilon }}-\bar{\mathbb {C}}^{-1}:\bar{\boldsymbol{\sigma }}\) and the effective strain rate is approximated by the finite difference between two subsequent increments indicated by the indices \(n-1\) and n
The angles \(\alpha ^\varepsilon \) and \(\theta ^\varepsilon \) are derived using (16.14) from the effective plastic strain rate, or \(\alpha \) and \(\theta \) are either given if a stress controlled RVE simulation is performed or derived from the simulated effective stress tensor using (16.14) for a strain controlled simulation. The data base of training samples used by Malik et al. [32] was generated using \(n_{\alpha }=39\), \(n_\theta =19\) and \(n_{\lambda }=20\) variations for \(\alpha \), \(\theta \) and \(\lambda \). The neural networks are fully connected feed forward neural networks containing a layer of input units (or neurons), one or more hidden layers and a layer of output units. The numbers of units for the input and output layer depend on the problem, whereas the number of hidden layers and units can be chosen freely. Each hidden neuron with index i represents a sigmoidal activation function depending on the input \(x_i\), which is the sum of weighted outputs of the neurons from the previous network layer plus a bias
The weights \(w_k\) and biases \(b_i\) are the free parameters of the network, which are determined during a training process. The necessary number of neurons within the hidden layers depend on the complexity of the function which the network is supposed to represent. A good approximation quality can be achieved if the number of free parameters of the neural network corresponds approximately to the square root of the number of training samples. To supervise the training process, the training data are split into two different batches. 90% of the data are used for the network training and the remaining 10% are solely used for the validation of the approximation quality, see Fig. 16.6. The numerical realization of the neural networks, the training and validation process is done with the Python package FFNET [34]. It allows especially the conversion of the neural network into Fortran code, which can be integrated in user material routines (UMAT) for the finite element program Abaqus (Fig. 16.6).
2 scatter plots. The left scatter plot plots target versus output. The data representing N N phi is concentrated on the top right and N N n is concentrated on the bottom left diagonally. The right scatter plot of e train m s e versus epochs plots 4 descending curves and a vertical line.
Fig. 16.6
a Accuracy of a trained network. b Evolution of training and validation mean square error during the learning process [32]
×
An alternative homogenization approach to model the inelastic behavior of foam structures is the FE2 method. The concept of FE2 is to run separate microscale FE simulations of the RVE on each macro material point. Recently, an efficient monolithic formulation was proposed by Lange et al. [10]. The FE2 method is supposed to produce most accurate homogenized results, however, it is still computational expensive compared to the presented hybrid approach, which resorts to training data that has been computed offline before.
In order to be able to capture size effects of foam structures [35] within the scale transition, higher gradient homogenization schemes have to be considered. Micromorphic theories [36, 37] allow to model these effects by extending the continuum model at the macro scale with additional degrees of freedom. Which of the available theories is the most dominant in the considered foam structures is part of ongoing research [38, 39].
16.2.3 Phase-Field Modeling of Multi-component/Multi-phase Systems
The phase-field method is a versatile tool to model the microstructure evolution of multi-phase systems, see [40] for an overview. Moreover, phase-field modeling approaches have been applied to various coupled problems such as electrochemistry [41], chemically reactive systems with phase separation [42] and, very recently, to chemo-mechanics focusing on various aspects [43‐48]. Available commercial multi-component/multi-phase simulation tools such as DICTRA employ a sharp interface description. The intentionally restricted applicability to 1D-problems could be justified so far by the number of practically relevant cases [49]. In addition, DICTRA has access to the mobility database of the Thermo-Calc software necessary to describe the diffusion processes. Furthermore, Thermo-Calc is directly called in order to determine the existing phases [49, 50]. However, the advantage of the phase-field method is its straightforward application to 2D and 3D-problems covering the evolution of complex phase-interface topologies in a diffusive sense.
Recently, a unifying framework for phase-field modeling of multi-component/multi-phase
chemo-mechanics has been proposed by Svendsen et al. [51]. Our phase-field approach is based on a mixture theory assuming \(i=[1,\dots ,M]\) components, where their individual masses \(m_i\) sum up to the total mass m. The constitution of the mixture is characterized by the mass fractions \(c_i\) of the components which obey additional sum relations. In addition, order parameters\(v_j\) are defined as mass fractions to distinguish \(j=[1,\dots ,K]\)phases, e.g., different aggregate states. The following relations hold
The evolution of the component mass fractions and the order parameters are described by coupled, generalized Cahn-Hilliard and Allen-Cahn type equations [52], respectively. The phase-field model considers the gradients \(\boldsymbol{\nabla }c_i\) and \(\boldsymbol{\nabla }v_j\) in order to capture the influence of phase interfaces in a diffusive sense. We derive the balance equations and necessary constitutive relations on the basis of an extended Coleman-Noll procedure for generalized continua as proposed by Hütter [53]. After introducing the conventional mass balance and diffusion equations, a gradient-extended energy balance is exploited in the framework of the thermodynamics of internal variables yielding the desired Cahn-Hilliard and Allen-Cahn equations. Thermodynamically consistent relations for the dissipative mechanisms are derived with help of a dissipation potential approach. Especially, we focus on the incorporation of chemical conversion processes based on the theory of equilibrium reactions. In Sect. 16.2.3, the phase-field model is re-cast into a mixed rate variational setting as proposed by Miehe [43, 54, 55] which enables a proper numerical treatment using finite elements.
Field Equations for Phase-Field Models
The mass balances of the single components can be cast into diffusion equations
for the evolution of the density, where \(\dot{{\boldsymbol{u}}}\) denotes the material (barycentric) velocity of the mixture. In order to ensure the balance of total mass we employ restrictions on the flux and source terms \({\boldsymbol{h}}_i\) and \(h_i\) [51], respectively,
The energy balance is discussed for the special case of mass transport and phase transition assuming isothermal conditions. The global balance for a spatial region \(\mathcal {P}_{\text {}}\) convecting with the material body \(\mathcal {B}_{\text {}}\) reads
with \(e\) being the mass specific internal energy and a flux of generalized power acting on the boundary \(\partial \mathcal {P}_{\text {}}\) with outward unit normal \({\boldsymbol{n}}\). Therein, \({\boldsymbol{\Pi }}_i\) and \({\boldsymbol{\Phi }}_j\) are power-conjugated microstresses to \(c_i\) and \(v_j\), respectively. Applying the divergence theorem to (16.36), we obtain the localized energy balance
The generalized flux terms introduce the gradients \(\boldsymbol{\nabla }\dot{c}_i\) and \(\boldsymbol{\nabla }\dot{v}_j\) into the energy balance. Next, the second law of thermodynamics is considered based on a balance of entropy and Gibbs’ fundamental equation [57, 58]. We assume a local entropy balance of the general form
where the restriction \(\pi \ge 0\) on the source term states the second law of thermodynamics [57]. The mass specific entropy is denoted as \(\eta \).
The entropy flux \({\boldsymbol{j}}\) and the source \(\pi \) are specified following the formalism of the thermodynamics of internal variables as comprehensively described by Lebon et al. [58]. To this end, the internal energy \(e\) is exchanged by the Helmholtz free energy \(\psi \) applying the established Legendre transformation \(e=\psi +\vartheta \,\eta \). The time derivative of the latter yields the localized Gibbs fundamental equation
which can be exploited to specify the terms of the entropy balance (16.38) [57, 58]. The free energy depends on the temperature \(\vartheta \), and additionally on the deformation gradient which is convenient for a chemo-mechanical extension [43, 51]. However, we introduce the following set of constitutive state variables
for the mass transport and phase transition processes. Note that this set of variables can be extended by kinematic deformation measures in case of chemo-thermo-mechanical processes. Here, we finally specify the Gibbs equation (16.39) inserting the time derivative of the free energy (16.41) together with the local energy balance (16.37) and the diffusion equations (16.33) as
Following the Coleman-Noll argumentation [53], the terms in brackets in (16.42) having rates of state variables as pre-factors could be chosen so that their dissipative contribution vanishes completely. On the other hand, to discuss the remaining terms in view of (16.38) some dissipative mechanisms could be additionally introduced as discussed by Hütter [53]. We define energetic constitutive relations for the entropy and the microstresses
where \(\mu _i\) and \(\Phi _j\) are generalized chemical potentials and dissipative microforces of the phase transition, respectively. The chemical potential yields the classical definition \(\mu _i=\partial _{c_i}{\hat{\psi }}\,\) for vanishing gradient terms in the free energy. As mentioned by Hütter [53], the additional balance equations (16.46)–(16.47) appear naturally without an a priori introduction of microforce balances as preferred by Gurtin [52]. Furthermore, boundary conditions based on microtractions can be prescribed which includes as a special case the generalized (trivial) boundary conditions utilized by Svendsen et al. [51].
Applying the sum relations (16.32) and (16.35) to the remaining terms in (16.42) yields
We interpret \(\mathcal {D}\) as the dissipation. To ensure \(\mathcal {D}\ge 0\), the microforces \(\tilde{\Phi }_j\), the mass fluxes \({\boldsymbol{h}}_i\) and the mass sources/sinks \(h_i\) must be chosen accordingly. The latter are related to chemical conversion processes which is addressed in the following.
The source terms \(h_i\) are specified for \(P\) parallel chemical reactions written in the generalized format
Therein, \(\textrm{X}_i\) denotes the chemical formula of the component i. According to the stoichiometry, a prefactor \(\nu _{i}^{\alpha }\) for every component i appears depending on the reaction \(\alpha \). Here, we follow the convention
Due to the stoichiometric relations and the mass conservation, the mass of components taking part in a chemical reaction \(\alpha \) can only change proportional to a single conversion rate density \(r^{\alpha }\). According to [56] we decided for linear relations
In order to ensure positive dissipation, linear Onsager-relations could be chosen between the conjugated variables in (16.55) as established by the thermodynamics of irreversible processes, see [56, 57]. In view of the desired variational setting, however, we proceed with a dissipation potential approach. For details we refer to the paper of Miehe [54] and the literature therein. In this spirit, a mass specific dissipation potential \(\phi \) with the rate \(\dot{{\underline{\boldsymbol{z}}}}\) as variables and state variables \({\underline{\boldsymbol{z}}}\) as additional parameters is introduced
The potential \(\phi \) is convex w.r.t. its arguments and zero at \(\dot{{\underline{\boldsymbol{z}}}}={\underline{0}}\). A set of generalized dual dissipative forces \({\underline{\boldsymbol{Z}}}\) is derived as
The mentioned properties of the dissipation potential \(\phi \) ensure \(\mathcal {D}_{\dot{z}} \ge 0\). Moreover, a dual dissipation potential\(\phi ^*\) can be defined
The viscosity parameters \(\beta _j\) and the parameters \(k^{\alpha }\) of the conversion rate densities are positive numbers. The Onsager-mobility tensors\(\boldsymbol{M}_{ab}\) obey \(\boldsymbol{M}_{ab}=\boldsymbol{M}_{ba}\). The microforces \(\tilde{\Phi }_j\) are derived from \(\phi \), whereas the fluxes \({\boldsymbol{h}}_i\) and the conversion rates \(r^{\alpha }\) stem from the dual potentials
At this point, the balance equations (16.33), (16.46), (16.47) and the constitutive relations (16.44), (16.45), (16.63)–(16.65) are combined to generate the field equations. These form the basis for the variational formulation. For the specific problems to be analyzed we assume that \(\rho = \rho _0\), i.e., the density of the mixture keeps its initial value and is additionally constant over the whole domain. Furthermore, we set \(\dot{{\boldsymbol{u}}}={\boldsymbol{0}}\). The field equations are written in terms of variational derivatives which are defined as
Note that \(2\left[ M-1\right] \) and \(K-1\) equations need to be solved due to the sum relations of the mass fractions. Moreover, the equations (16.70) have the structure of generalized Allen-Cahn equations . Combining (16.69) with (16.68) yields generalized Cahn-Hilliard equations . In the previous set of equations, relative quantities w.r.t. the last component \(M\) or last phase \(K\) are denoted by a tilde, e.g.,
with prescribed values at the respective boundaries marked with a bar.
Mixed Rate-Type Variational Setting
Following Miehe [54, 55], the field equations (16.68)–(16.70) appear as the Euler-Lagrange equations of a
variational problem based on the mixed rate potential
The necessary condition for stationary points reads \(\delta \,{\Pi ^*\left( \dot{c}_i,\dot{v}_j,\tilde{\mu }_i\right) }=0\), where \(\delta \,{\Pi ^*\left( \dot{c}_i,\tilde{\mu }_i,\dot{v}_j\right) }\) denotes the first variation of the rate potential. We allow arbitrary variations of the defined field variables, except at boundaries where the field variables are fixed. Note that variations of \(\delta \,{\!}{\dot{c}_i}\) and \(\delta \,{\!}{\dot{v}_j}\) are restricted by the sum relations (16.32). Therefore, we finally find (16.68)–(16.70) as the Euler-Lagrange equations.
For the numerical treatment via the finite element method, a time incremental counterpart of the variational principle is derived as explained in detail by Miehe et al. [43, 55]. For discretization in time, an Euler-backward scheme is utilized. The spatial discretization is performed by iso-parametric finite elements with low-order, C0-continuous shape functions for all considered primary field variables, which are the mass fractions \(c_i\), the chemical potentials \(\tilde{\mu }_i\) and the order parameters \(v_j\). The final non-linear system of algebraic equations to be solved for any time point t is handled with a monolithic Newton-scheme. Due to the variational structure, the considered problem provides an inherently symmetric tangent matrix as pointed out in [43, 55].
16.3 Results and Applications
16.3.1 Application of Foam Models
The generic foam structures described in Sect. 16.2.1 were widely used within the CRC 920 either as geometric models in numerical simulations or as 3D printed structures to conduct specific experiments. Foam RVEs were used to simulate the influence of foam morphology on effective properties (hydraulic turtuosity, viscous and inertial permeability, filtration coefficient) related to metal melt filtration [11]. Asad et al. [59] investigated the immersion process of a ceramic filter in a steel melt. Lehmann et al. [12] used the RVEs as base structures and investigated the influence of specific geometric modifications like additional struts, closed foam windows or streamlined strut cross sections on effective hydraulic properties and filtration performance. 3D printed filter structures based on the RVEs have been used by Wetzig et al. [60] for real world filtration experiments, whereas Bock-Seefeld et al. [25] and Herdering et al. [24, 61] used the RVEs to produce polymer filter templates for customized ceramic foam structures to estimate filtration efficiency, structural filter strength, and integrity. A comparison of mechanical properties between generic and real foam structures performed by Settgast et al. [62] proved the high accuracy of the generic foam models.
16.3.2 Thermo-mechanical Behavior of Foams and Filter Structures
The generation procedure of the foam model allows modifications of the strut shape as explained in Sect. 16.2.1. The influence of the strut shapes on the elastic properties was investigated by Storm et al. [63] on spatial periodic Kelvin cells. The corresponding FE-models were either build from beam or 3D elements. We found that beam models systematically underestimate the stiffness of foams, since the additional reinforcement of the foam nodes cannot be represented by beam elements. However, volumetric models follow almost perfectly the Gibson-Ashby relation [64, 65] for the relative elastic modulus depending on the relative density as \(\bar{E}/E=\left[ \bar{\rho }/\rho \right] ^2=\rho _\text {rel}^2\).
A set of 2 3-D structures. The left structure is a generic model with a smooth skeleton and well-defined pores. The right structure is a C T scan, highlighting multiple depressions on the skeleton and irregularly shaped pores.
Fig. 16.7
Comparison of a generic foam structure (left) with a relative density of \(\varrho = 20\%\), and a CT-scan image (right) of a 10ppi foam sample. For better comparison, all elements of the generated model lying outside the bounding cylinder (\(d = 20\,\textrm{mm}\), \(h = 20\,\textrm{mm}\)), are removed [1]
×
Furthermore, Storm et al. [66] investigated the influence of the strut shape on the strength of a foam. The cross section radius of the struts was varied according to a hyperbolic equation along their longitudinal axis, where \(r_{\text {e}}/r_{\text {m}}\) defines the ratio of the strut radii at the end and the middle of the strut. A ratio \(r_{\text {e}}/r_{\text {m}}=1\) characterizes a strut with constant cross section along its axis. Struts with \(r_{\text {e}}/r_{\text {m}}>1\) have thicker foam nodes and the cross sections become smaller towards the strut centers, which is the typical form in real struts as shown in Fig. 16.7. Interestingly, an optimal ratio can be found for \(r_{\text {e}}/r_{\text {m}} \approx 1.4\) maximizing the bending strength of the struts, illustrated in the corresponding Fig. 16.8a. The curvature of struts measured in terms of the ratio between length and curvature radius \(l_0/r_{\text {c}}\) has only a minor influence on the strength of foams, which is shown in Fig. 16.8b.
Further results from Storm et al. [66] include the comparison of continuum and hybrid models, where only the foam nodes where meshed using continuum elements, whereas large parts of the struts are modeled using beam elements. The hybrid models show an elastic behavior comparable to the continuum models, but the numerical effort required to compute them is considerably reduced.
2 scattered line graphs. The left graph plots 6 curves with an ascending region, a peak, and a descending region. The right graph plots 3 curves and 3 sets of scattered data in descending trend.
Fig. 16.8
Influence of the a material distribution along the strut axis (hyperboloid end to middle radius \(r_e/r_m\)) and b strut curvature (strut length to curvature radius \(l_0/r_c\)) on the norm of critical deviatoric stress in \(\pi \)-plane for continuum and beam models of the Kelvin cell [66]
×
Small sharp edged cavities remain along the center line of the struts due to the production process. It was found that the influence of the cavities does not have a significant effect on the elastic properties of the foams as long as the relative density remains constant [31].
Failure limit surfaces allow the assessment of strength of filter structures utilizing different failure criteria. Firstly, we determined yield surfaces analytically for Kelvin cells where the struts are assumed to follow the Euler-Bernoulli beam theory. Using this model, results for different failure criteria as von Mises, maximum absolute principal, maximum principal and principal stress criterion with tension compression asymmetry have been calculated [29]. Storm et al. [31] studied the influence of geometrical variations for a Kelvin cell foam topology. These investigations considered the change of cross-section shape of the struts, strut curvature and pore anisotropy. All these features influence the elastic properties significantly. For the generic and almost isotropic foam structure F216 as shown in Fig. 16.3, investigations tackling failure assessment have been conducted using the von Mises criterion, considering different relative densities and strut geometries [1]. Two examples of yield surfaces of the F216 foam RVE are displayed in Figs. 16.9 and 16.10 for a von Mises and Weibull failure criterion, respectively. The macroscopic failure surface for a local von Mises criterion, as shown in Fig. 16.9, is point symmetric with respect to the origin. The shape of the meridian cross section depends on the Lode angle \(\theta \). For a Lode angle \(\theta =0\) the maximum \(J_2\)-value is located where the hydrostatic stress vanishes (\(I_1=0\)). For negative or positive Lode angles the maximum is shifted towards the negative and positive hydrostatic stress values, respectively. The rightmost diagram shows cross sections for different angles \(\alpha \) (cf. (16.11)) projected onto the deviatoric plane. For \(\alpha =0\), which implies a pure deviatoric stress state, the cross section exhibits a hexagonal shape, whereas for positive or negative angles the cross section transforms towards a triangular shape but with opposite orientation. To the authors knowledge, for such shapes of failure surfaces no closed form analytical description is available in literature.
3 graphs. The left graph is a 3 D surface plot, with the 3 plot peaking when the 3 variables are 0. The center graph is a multi-line graph with the theta = 0 curve normal distribution, negative pi by 6 curve left skewed, and pi by 6 curve right skewed. The right graph is a radial graph of alpha.
Fig. 16.9
Failure limit surfaces for a F216 foam RVE with a relative density of 20% for a von Mises criterion. (Left) limit surfaces in principal stress space. (Center) meridian cut for different Lode angles \(\theta \). (Right) hydrostatic cuts at different angles \(\alpha \)
A set of 3 graphs. The left graph is a 3 D surface plot, with the 3 plot peaking when the 3 variables are 0 and forming closed surfaces. The center graph is a multi-line graph with the theta = negative pi by 6 curve left skewed and pi by 6 curve right skewed. The right graph is a radial graph.
Fig. 16.10
Failure limit surfaces for a F216 foam RVE with a relative density of 20% for a Weibull criterion. (Left) limit surfaces in principal stress space. (Center) meridian cut for different Lode angles \(\theta \). (Right) hydrostatic cuts at different angles \(\alpha \)
×
×
The failure surface for the Weibull criterion (cf. Fig. 16.10) given at a failure probability of \(63.2\%\) shows a strong asymmetry with respect to the hydrostatic stress axis. The strength in hydrostratic compression is much higher than in hydrostatic tension. The ratio between compressive and tensile strength depends on the relative density of the foam. The meridian cross sections show a similar contour as for the von Mises criterion, but the maximum value for negative Lode angles exceeds the one for positive Lode angles. In general the maxima of the meridian cross sections are observed for negative hydrostatic stress values (\(I_1<0\)) and exceed the ones of the von Mises criteria. The cross section shapes for constant angles \(\alpha \) vary from an almost circular shape for \(\alpha =0.425\) to triangular shapes down to \(\alpha \approx -0.125\). Around \(\alpha \approx -0.25\) an almost hexagonal shape is observed, before for further decreasing \(\alpha \) values the shape morphs again towards a triangle, but with opposite orientation compared to those for positive angles \(\alpha \).
The spatial orientation of non-isotropic foam RVEs with respect to the principal directions of the effective stress tensor has an influence on the shape of the corresponding failure surfaces and has been investigated by Zhang and Storm for the special case of orthotropic Kelvin cells [28, 29]. For RVEs containing an increasing number of randomly arranged pores the effective behavior tends towards isotropy and the orientation of the RVE becomes negligible.
The sharp edged cavities remaining from the production process required a fracture mechanical analysis, conducted by Settgast et al. [30]. With help of the interaction integral (16.19), the SIFs are computed for each point along the sharp edges of a Kelvin cell as illustrated in Fig. 16.5. Since the FE analysis is linear elastic, any load case and therewith the corresponding SIFs can be constructed by a scaled superposition of these six independent base load cases. It was found that the fracture limit surfaces always enclose the failure limit surfaces for a von Mises or a maximum principal stress criterion using credible values obtained by experiments for the considered material. From these observations it can be concluded that the local failure is preferably triggered by a critical principal stress on the outer strut surface rather than by the stress concentration at the crack front inside the strut cavities. Figure 16.11 shows exemplary the values of the three SIFs along the local parametric coordinate \(l_{\text {pc}}\) of a sharp edge around a square shaped foam window for an effective shear loading.
A set of 2 line graphs plot K 1, K 2, and K 3 versus l p c curves. Left graph. K 2 and K 3 curves (0, 0), (2, 0), and (4, 0). The K 1 curve (0, 0), (1.5, 1.1), (2.5, 1.2), (4, 0). Right graph. K 1 (0, 0), (1, 0), (4, 0), K 2 (0.4, 0.5), (2.6, 0.4), (4, 0), K 3 (0, 1.2), (2.5, 2.5), (4, 1.2).
Fig. 16.11
SIFs along the sharp edged cavity within the four-sided strut loop of a Kelvin cell having the normal (001) under effective a uniaxial tension \(\bar{\sigma }_{11}=1/16\,\textrm{MPa}\) and b shear loading \(\bar{\sigma }_{13}=1/16\,\textrm{MPa}\) [30]
3 scatterplots. A plots J positive versus h data as follows. J positive values increase with the increase in h. B plots J positive versus rho relative percent. J increases with an increase in rho relative. C plots J positive versus l and J increases with an increase in l.
Fig. 16.12
a\(J^+(h)\) for \(\rho _\text {rel}=10\%\), \(l=4\,\textrm{mm}\) and different \(E(\vartheta )\), b\(J^+(\rho _\text {rel})\) for \(h~=~2000\,\textrm{Wm}^{-2}\textrm{K}^{-1}\) and \(l=4\,\textrm{mm}\) and c\(J^+(l)\) for \(h=2000\,\textrm{Wm}^{-2}\textrm{K}^{-1}\) and \(\rho _\text {rel}=10\%\) taken from [23]
×
×
A more critical scenario is thermal shock loading, which occurs at the beginning of the casting process when the melt enters the filter, which may cause the filter to fail even before the melt filtration process starts. Here, the sharp edged cavities are potentially the most critical locations. The occuring stress intensity factors or J-integral values depend strongly on the material properties of the foams bulk material. The dependencies of the maximal J-integral value \(J^+\) along all cavity edges on thermal conductivity h, relative density \(\rho _\text {rel}\) and RVE length l are presented in Fig. 16.12.
\(J^+\) depends strongly on the thermal conductivity h with an almost quadratic dependence \(J^+ \propto h^2\) , because the increasing heat flux leads to higher temperature gradients and therefore higher thermal stresses as well as a higher fracture mechanical loading at the crack front. The influence of the different elastic moduli of the three considered filter materials slightly increases with an increasing heat transfer coefficient. If the relative density \(\rho _\text {rel}\) is increased the size of the sharp-edged cavities does not change. However, the outer surface area \(A_{\text {o}}\) of the foam increases and therewith the heat flux and the fracture mechanical loading at the cracks, with \(J^+ \propto \rho _\text {rel}^{1/\kappa }\) with the value of the denominator of the exponent \(\kappa \gtrsim 1\). The largest impact on \(J^+\) is observed for a changing size of the RVE. It is a combination of two effects. First the crack size changes linearly with the cell size and \(J \propto l\). The second effect is that the foam surface area \(A_{\text {o}} \propto l^2\). Both effects result in an almost cubic dependence \(J^+ \propto l^3\) . The detailed material parameters and the temperature dependent elastic moduli as well as the loading conditions can be found in reference [23].
In some applications, e.g., continuous casting, the ceramic foam filter is subjected to a permanent loading caused by the molten metal for more than a couple of minutes. For that reason, it is important to estimate the long time behavior of the material at high temperatures. In Settgast et al. [67] the creep properties of foams are compared with the creep properties of the bulk material. The following stress-time dependent creep law was assumed
for the macro scale, wherein the effective equivalent creep strain rate \(\dot{\bar{\varepsilon }}^{\text {cr}}_\text {eq}\) depends on the effective equivalent stress \(\bar{\sigma }_\text {eq}\), the structure specific effective stress factor \(\bar{A}\), with exponent n and time t with its corresponding exponent m. When applying the same creep law for the bulk material at the micro scale, a good agreement between experimental and simulated creep curves is achieved for uniaxial loading at elevated temperature (cf. Fig. 16.13). Interestingly, the exponents n and m are the same for the local and homogenized constitutive material laws. The ratio between effective and local stress factor
depends on the structure via exponent a and the relative density \(\rho _{\text {rel}}\) of the foam. The exponent a is used to describe the proportionality
which is fitted by multiple FE analyses where the relative density is varied. Assuming a constant relative density the creep rate can be minimized if the strut cross sections are constant along the strut axis [62, 67].
A line graph plots epsilon c r versus t simulation and experiment curves. As the A overbar and A ratio increase and the rho relative increases, the slope of the curves decreases. The curves are right-opening parabolas with an origin at (0, 0) and a decreasing trend.
Fig. 16.13
Comparison of creep curves of three tested and simulated foam samples for uniaxial loading (\(\bar{\sigma }=0.1\,\textrm{MPa}\) and \(\vartheta =1350\,^\circ \)C [62]
×
Utilizing an FE2 approach in a monolithic setting, proposed by Lange et al. [10], for which a separate microscale FE simulation of the RVE is run on each macro material point, the influence of foam morphology on the strength and mechanical creep behavior can be investigated during realistic filtration scenarios. In [22] a flow-through filter application (cf. Fig. 16.14) was investigated. As foam RVE the Wheire-Phelan cell is utilized, because of less computational costs compared to the F216 foam RVE. The filter causes a pressure drop and therewith mechanical forces inside the filter. For a constant volumetric flow rate but increasing relative density the pressure drop increases, which can be described by the Darcy-Forchheimer law. The developed model allows the prediction of elastic and creep deformations for realistic filtration scenarios. The shape of the foam struts has almost no influence on the elastic deflection for a constant relative density. However, with an increasing relative density \(\rho _{\text {rel}}\) the deflection decreases. An increasing Reynolds number of the flow leads to an increasing deflection, due to the larger pressure drop [11, 22]. In Fig. 16.14 the deformation including creep after one hour of a flow through filtration scenario is illustrated with the underlying FE2-micro model.
A deformation u magnitude profile on the left and a zoomed-in view of a section of the ultrastructure of the foam on the right. The left profile has a central focus of high u magnitude which decreases radially outward. The right profile represents the sigma and E ratio along the foam structure.
Fig. 16.14
Deformation of a filter during a flow through scenario utilizing the FE2 homogenization approach [22]
×
The hybrid approach described in Sect. 16.2.2 employing neural networks was developed to further reduce the computational costs. First, the feasibility of this approach was examined in 2D on a Kelvin cell by Settgast et al. [68]. The model was then extended by a damage formulation [33, 69]. In Fig. 16.15 stress-strain curves are compared from simulations using a fully discretized RVE, with those predicted by the proposed hybrid approach. The hybrid approach can reproduce the RVE simulations almost exactly, even if elastic unloading is taken into account. Also the material degradation due to damage can be reproduced. For these simulations a speed-up factor of 4000 compared to the fully resolved RVE was achieved.
3 graphs. A plots sigma bar 22 over sigma Y versus epsilon bar 22 E over sigma Y. B plots sigma bar 11 over sigma Y versus epsilon bar 11 E over sigma Y. C plots sigma bar 22 over sigma Y versus epsilon bar 22 E over sigma Y. All graphs plot a set of parallel lines joined by an outer line.
Fig. 16.15
Macroscopic stress-strain curves for a uniaxial loading \(\bar{\varepsilon }_{22}\), and b biaxial loading \(\bar{\varepsilon }_{11} = \bar{\varepsilon }_{22}\), where damage is considered. E and \(\sigma _\text {Y}\) denote the elastic modulus and yield stress for a local von Mises plasticity model, respectively [33]
A set of 5 illustrations of yield surfaces in principal stress space. A illustrates a smooth, elliptical 3 D shape. B includes meridian cuts. C, D, and E provide zoomed-in views of the meridian cuts from B with 3 axes each labeled sigma 1, sigma 2, and sigma 3.
Fig. 16.16
Determined yield surfaces
of the investigated Wheire-Phelan RVE for different values of \(\bar{\varepsilon }_{\text {q}}^{\text {pl}}=0.005\), 0.05, 0.1, 0.15, 0.2, 0.25 utilizing the approximated yield function \(\text {NN}^{\Phi }\) in principal stress space. The thin lines indicate the yield surfaces for different values of \(\bar{\varepsilon }_{\text {q}}^{\text {pl}}\). The symbols \(\sigma _1\), \(\sigma _2\) and \(\sigma _3\) denote the three principal stresses. a represents the 3D yield surface, b shows the meridian cut with marked deviatoric cutting lines c, d, and e along the hydrostatic axis in \(I_1^{\sigma }\)-direction without, negative and positive hydrostatic part, respectively taken from [32]
×
×
Malik et al. [32] extended the hybrid approach to model also 3D material behavior. As RVE a Wheire-Phelan foam as shown in Fig. 16.3 is used. The corresponding yield surface
was identified and its evolution due to strain hardening could be depicted, which is displayed in Fig. 16.16. The comparison of simulations using a fully discretized RVE with the results of the hybrid approach was conducted as well. Good agreement was achieved (cf. Fig. 16.17) for different load cases and even for cyclic loading, despite the fact that information on cyclic load cases have not been part of the training data set for the neural networks. It should be emphasized that the flexibility of the neural networks allows the exact reproduction of the short \(\bar{\sigma }_{11}\)-drop after reaching the maximum load in the elastic regime before hardening, as noticed in the simulations using fully resolved RVEs.
A set of 4 line graphs plot sigma bar slash sigma 0 versus epsilon bar E slash sigma 0 curves for sigma bar 11 N N, sigma bar 22 N N, sigma bar 33 N N, sigma bar 11 R V E, sigma bar 22 R V E, and sigma bar 33 R V E in different trends.
Fig. 16.17
Stress-strain curves of proportional a uniaxial tension \(\Delta \overline{\varepsilon }_{11}=0.1,\, \Delta \overline{\varepsilon }_{22}=\Delta \overline{\varepsilon }_{33}=0\), b uniaxial compression \(\Delta \overline{\varepsilon }_{11}=-\,0.1,\, \Delta \overline{\varepsilon }_{22}=\Delta \overline{\varepsilon }_{33}=0\), c deviatoric shear \(\Delta \overline{\varepsilon }_{11}=-\,2\,\Delta \overline{\varepsilon }_{22}=-\,2\,\Delta \overline{\varepsilon }_{33}\) and d cyclic deviatoric loading \(\Delta \overline{\varepsilon }_{11}=-\,2\,\Delta \overline{\varepsilon }_{22}=-\,2\,\Delta \overline{\varepsilon }_{33}\) [32]
×
16.3.3 Verification of the Phase-Field Model
In order to verify the FE-implementation we consider a reactive two phase/two component system, where a dissociation reaction
between the fictitious chemical components \(\textrm{A}_{2}\) (diatomic, \(i=1\)) and \(\textrm{A}\) (monoatomic, \(i=M=2\)) can occur. The free energy of the mixture is chosen according to a regular solution model with binary interactions
Therein, \(\psi _{\text {bul,id}}\) and \(\psi _{\text {bul,ex}}\) denote the phase dependent energy parts of the ideal solution and excess model, respectively. Moreover, \(\psi _{\text {int,c}}\) and \(\psi _{\text {int,v}}\) comprise the interface energies between phases characterized either by different component mass fractions \(c_i\) or order parameters \(v_j\). The order parameter \(v_2\) is directly substituted by \(v_2=1-v_1\) in the given energy. Moreover, \(v_1\) appears via the interpolation function \(p\left( v_1\right) \) which mixes the chemical potentials of the pure substances in the respective phases. Furthermore, the double well function \(g(v_1)\) ensures extrema of the bulk energy at \(v_1=0\) and \(v_1=1\)
appears in the logarithmic term of ideal solution and in the excess energy in order to be consistent with a description in molar fractions, as used, e.g., by Bai et al. [47]. The mobility tensors are simplified as \(\boldsymbol{M}_{ab}=M_{ab}\boldsymbol{I}\). A constant conversion rate \(k^1\) for reaction (16.84) is considered. The studied model parameters are summarized in Table 16.1.
Table 16.1
Model parameters of the verification example (dissociation reaction & phase separation)
are plotted for \(v_1=1\) (phase 1) and \(v_1=0\) (phase 2). By arguments of equilibrium thermodynamics, possible final states of the model system are known for different initial component/phase compositions which can be illustrated with help of Fig. 16.18a. We firstly assume that the mass fraction of \(A_2\) is homogeneous across the domain with \(c_1=0.92\) at the beginning. The bulk energies can be read off at point 0 in Fig. 16.18a. For a fixed composition, i.e., no chemical reaction, the energy of the system can be minimized, if the system is allowed to decompose into sub-regions with compositions given by the common tangent construction (dotted line in Fig. 16.18a) and mass fractions according to the lever rule [70]. The energetic state of the system with decomposed phases can be found on the common tangent, i.e., point 2 in Fig. 16.18a. During a simulation, this behavior is only observable if the initial order parameter is not constant \(v_1=0\) or \(v_1=1\) throughout the system, i.e., an initial perturbation is necessary. For chemically reactive systems, the mass fractions can evolve, e.g., along the energetic states 1-2-3, until the global minimum of the bulk energy which is point 4 in Fig. 16.18a, with a homogeneous composition and stable phase 2 everywhere. At the intermediate point 3, phase 1 vanishes completely.
2 line graphs. A plots normalized bulk energy versus mass fraction psi superscript j bul c 1 of component A 2 curves for phase 1 and phase 2. Phase 1 and phase 2 curves have an ascending trend and coincide at (0.92, 0.4). B plots 4 descending parabolic normalized energy psi k versus normalized time curves.
Fig. 16.18
Energy landscapes of the verification example: a Bulk energies of the considered phases, and b time evolution of system averaged energy contributions. Numbers refer to bulk energies at particular time points
×
For the 2D-simulation, a square domain is considered with edge length \(l_0=10\) m. As mentioned, the initial mass fraction \(c_1\) of \(A_2\) is homogeneously prescribed:
Trivial Neumann boundary conditions are assumed for the mass flux as well as for the microtractions at all boundaries. The domain is discretized by 100 \(\times \) 100 = 10\(^\text {4}\) quadrilateral finite elements of equal size.
The process to reach energetic minimum highly depends on the chosen model parameters. For the considered case, phase transition and diffusion proceed much faster than the chemical reaction. In Fig. 16.18b, the evolution of different contributions to the system’s whole energy are plotted over time. The specific energies are illustrated as normalized quantities \(\bar{\psi }_{k}\) given by
where the index k refers to the terms defined in Eqs. (16.85)–(16.89). The bulk energy \(\bar{\psi }_{\text {bul}}\) is introduced as the sum of the ideal solution and excess parts, \(\bar{\psi }_{\text {bul,id}}\) and \(\bar{\psi }_{\text {bul,ex}}\), respectively. Furthermore, \(\bar{\psi }\) is the normalized total free energy.
The numbers 0–4 in Fig. 16.18b belong to specific states of the system’s bulk energy and time points which are correspondingly highlighted in the energy landscape Fig. 16.18a. Additionally, the spatial distribution of the order parameter \(v_1\) and the mass fraction \(c_1\) of \(\textrm{A}_{2}\) are plotted in Fig. 16.20 at these time points. The plots belong to cuts along the x-axis at \(y=0\). These reduced representations contain the full field information of the considered primary variables since the fields just vary along x for the specific problem as shown in Fig. 16.19 for a particular time point.
A set of 2 spatial distribution charts. A plots y versus x data and the distribution of v 1 is as follows. Up to x = 3, v 1 = 0 to 0.5. From 3.2 onwards, v 1 = 0.8 to 1. B plots y versus x data and the distribution of c 1 is as follows. Up to x = 3, c 1 = 0.65 to 0.7. From 3 onwards, c 1 = 0.8 to 1.
Fig. 16.19
Spatial distribution of a the order parameter and b the mass fraction \(c_1\) of component \(\textrm{A}_{2}\) at time point 2 marked in Fig. 16.18
2 line graphs. A plots order of parameter v 1 versus x curves for 0, 1, 2, 3, and 4. The 0 curve is a step curve. The curves representing 1, 2, 3, and 4 are sigmoid curves. B plots mass fraction c 1 of A 2 versus x curves for 0 and 4 which are lines parallel to x, and sigmoid curves for 1, 2, and 3.
Fig. 16.20
Spatial distribution of a the order parameter and b the mass fraction \(c_1\) of component \(\textrm{A}_{2}\) at time points 0–4 marked in Fig. 16.18
×
×
According to Fig. 16.18b, the total free energy of the system is monotonously minimized during time evolution until the final energetic state 4 which matches the expectation from equilibrium consideration as highlighted in Fig. 16.18a. The following intermediate states of interest should be discussed. From point 0 to 1, the interface energy \(\bar{\psi }_{\text {int,v}}\) is quickly decreased, see Fig. 16.18b, in order to relax the initially sharp gradient of \(v_1\) to a typical \(\tanh \)-shape of the interface [71], compare Fig. 16.20a. The second gradient energy \(\bar{\psi }_{\text {int,c}}\) participates, too, but the change is not visible in the diagrams since the chosen gradient parameter seems to have a lower energetic effect. The considered domain consists of two distinct phases at point 1, but the energetic state and the mass fraction \(c_1\) in the phases are still above the common tangent and far away from the intermediate equilibrium constitution according to Fig. 16.20b. Due to the high mobility of diffusion and the slow conversion rate of the chemical reaction, the common tangent is reached at point 2 while the overall composition does not change during this short time span. The interface has moved to \(x\approx 3\) mm yielding a larger region of \(v_1=1\) according to the lever rule. The mass fractions \(c_1\) in the phases located at the left and the right from the interface exhibit the expected equilibrium values at the common tangent. From point 2 to 3 in Fig. 16.18b, the chemical reaction becomes dominant which decreases the overall mass fraction \(c_1\) of \(\textrm{A}_{2}\) shifting the bulk energy along the common tangent as additionally illustrated in Fig. 16.18a. Correspondingly, the domain of phase 1 shrinks, see Fig. 16.20a, but the mass fractions far from the interface keep the equilibrium values, see Fig. 16.20b. At point 3, the interface starts to decay which is visible as the jump-like change in the interface energy \(\bar{\psi }_{\text {int,v}}\) in Fig. 16.18b. The common tangent is left and phase 2 becomes stable during the rest of the purely reactive process. The mass fraction \(c_1\) takes a constant value across the spatial domain and the reaction decreases \(c_1\) to 0.13 as shown in Fig. 16.20b. This state is found at the global energetic minimum marked as point 4 in Fig. 16.18a. Here, the affinity \(A^1\propto \tilde{\mu }_1\) as driving force of the reaction consequently vanishes which defines the (dynamic) chemical equilibrium.
The obtained final state of the simulation matches the equilibrium thermodynamics solution, where intermediate evolution steps towards the chemical equilibrium are in qualitative accordance with our expectations explained above. The proposed two phase/two component problem is a versatile benchmark example to verify the numerical treatment relying at least on knowledge and solutions of equilibrium thermodynamics. For the same purpose similar simplified problems have been discussed in recent literature especially for Cahn-Hilliard type problems [47, 55]. Such additional cases have been successfully tested, too, e.g., chemical reaction during spinodal decomposition (Cahn-Hilliard type) for fixed phase 1. However, in this paper we focused on simultaneously active chemical reaction and phase transition for the sake of brevity. We conclude that our FE-implementation of the multi-component/multi-phase model yields trustworthy results and can be utilized for more complex problems.
16.4 Conclusions
Along the road tools have been developed to virtually design foam structures. These models can be used in simulations to predict and analyze their thermo-mechanical properties. Using additive manufacturing methods the models can be physically realized and used in experimental investigations regarding filtration phenomena. The conducted numerical investigations provide insights into structure-property relations of filter structures, but can now also contribute statements about the behavior of the entire foam structure in a homogenized manner. Different homogenization approaches can be applied to efficiently predict effective elastic, plastic, creep, fracture, and damage properties of foam structures under loading conditions in metal melt filtration applications. Geometric and topological variations of foams are discussed to improve the thermo-mechanical integrity and the filtration efficieny. Further, thermo-chemical phenomena of filtration processes, as the formation of in-situ layers can be modeled with developed phase-field models. These new models are able to describe diffusion, phase transition, and chemical reactions in the multi-component/multi-phase systems. The in-situ layer formation affects the reactive filtration phase. The knowledge about the kinetics of this effects is very valuable. The developed modeling tools and approaches are not restricted to applications regarding metal melt filtration. They can be applied in many different fields of thermo-chemo-mechanics and research on porous media or meta materials.
Acknowledgements
The authors gratefully acknowledge the work and efforts of Dr. -Ing. J. Storm, Dr. -Ing. C. Settgast and Dr. -Ing. Dongshuan Zhang done in the previous funding periods of the CRC 920. This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)—Project-ID 169148856—CRC 920: Multi-Functional Filters for Metal Melt Filtration—A Contribution towards Zero Defect Materials, subproject B05. Furthermore, the authors acknowledge computing time on the compute cluster of the Faculty of Mathematics and Computer Science of Technische Universität Bergakademie Freiberg, operated by the computing center (URZ) and funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)—Project-ID 397252409.
Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.
The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.
Die im Laufe eines Jahres in der „adhäsion“ veröffentlichten Marktübersichten helfen Anwendern verschiedenster Branchen, sich einen gezielten Überblick über Lieferantenangebote zu verschaffen.
