Scipedia

Polycarpou, P.

Papaloizou, L.

Kyriakides, M.

An Experimental Study for the Assessment of Numerical Impact Models for Simulating Earthquake Induced Pounding

Cerbe, F.

Hühne, C.

Sinapius, M.

Junge, A.

The 9th European Congress on Computational Methods in Applied Sciences and Engineering 
ECCOMAS Congress 2024 
3 – 7 June 2024, Lisbon, Portugal 
 
 
 
Simulation of SMPs for the bonding process of a gyroid-core and a CFRP 
aerodynamic profile 
Ferdinand Cerbe*¹, Alexander Junge¹, Christian Hühne² and Michael Sinapius³ 
1 Institute of Mechanics and Adaptronics, TU Braunschweig, f.cerbe@tu-braunschweig.de  
2 German Aerospace Center (DLR), Christian.Huehne@dlr.de 
3 Institute of Mechanics and Adaptronics, TU Braunschweig, m.sinapius@tu-braunschweig.de 
 
Key Words: CFRP bonding, Shape-memory polymers 
 
 
Additive manufacturing (AM) allows the manufacturing of highly complex structures, such as 
lattices or structures derived from triply periodic minimal surface (TPMS), such as the Gyroid. 
The TPMS can be used as a core to increase the specific stiffness and strength of a sandwich 
structure, which are relevant in aerospace engineering. A significant challenge is the 
manufacturing process, where the TPMS core must be inserted and bonded to an aerodynamic 
profile made of carbon fibre-reinforced plastic (CFRP). This integration can be achieved by 
using the shape memory effect of a polymeric AM material. The idea is to program the TPMS 
core to a compact shape, insert it into the airfoil, and eventually expand it during the bonding 
process to apply pressure on the bonding area. 
 
In our contribution, we present a thermo-mechanical Abaqus simulation of a shape memory 
polymer (SMP) based on a viscoelastic material model with temperature dependence [1]. In our 
presentation, we describe the straightforward process from material characterization by 
dynamic mechanical analysis to the simulation of the shape memory effect. In conclusion, we 
show the applicability of our approach in the bonding process of a gyroid core to a CFRP 
aerodynamic profile. 
 
Acknowledgments. This research was funded by the German Federal Ministry for Economic 
Affairs and Climate Action (BMWK) under project number 20A2103D (MuStHaF). 
 
REFERENCES 
[1] W. A. Azzawi, J. A. Epaarachchi, M. Islam, J. Leng.  Implementation of a finite element analysis procedure for structural 
analysis of shape memory behaviour of fibre reinforced shape memory polymer composites. Smart Materials and Structures, 
26: 125002, 2017. 
 
[2] H. Traub, M. Sprengholz, D. Teufel, C. Hühne. Structural-mechanical characterisation of triply periodic minimal surface 
sheet networks: simulation and experiment. AIAA SCITECH 2023 Forum. https://doi.org/10.2514/6.2023-2076 
 
[3] M. L. Williams, R. F. Landel, J. D. Ferry. The Temperature Dependence of Relaxation Mechanisms in Amorphous 
Polymers and Other Glass-forming Liquids. J. Am. Chem. Soc. 1955, 77, 14, 3701–3707. 
https://doi.org/10.1021/ja01619a008 
 
[4] N. W. Tschoegle. The Phenomenological Theory of Linear Viscoelastic Behavior: An Introduction. ISBN: 978-3-642-
73602-5. 2012. 
 
[5] W. Wagermaier, K. Kratz, M. Heuchel, and A. Lendlein. Characterization Methods for Shape-Memory Polymers. In: A. 
Lendlein (eds) Shape-Memory Polymers 2010 https://doi.org/10.1007/12_2009_25 


Simulation of SMPs for the Bonding Process of a Gyroid-Core with a CFRP Aerodynamic Profile

Bartolini, G.

Falco, A.

Landi, F.

Sequential Sensitivity Analyses on Masonry Mechanical Parameters in the Substructure Analysis of a Monumental Construction

Costa, S.

Bastos, L.

Rietter, L.

Oliveira, R.

Freitas, R.

Rocha, A.

Vale, B.

Ribeiro, C.

Serrão, D.

Silva, J.

Gonçalves, N.

Carneiro, F.

Silva, S.

Portinha, A.

Bernardo, P.

Dias, G.

Numerical simulation and experimental evaluation for the development of a light electronic enclosure used in automotive industry

Pingaro, N.

Buzzetti, M.

Milani, G.

On the nonlinear behaviour of masonry cloister vaults: A case study in Piacenza, Emilia Romagna, Italy

Lavore, E.

Ruano, J.

Oliva, A.

Trias, X.

The 9th European Congress on Computational Methods in Applied Sciences and Engineering
ECCOMAS Congress 2024
3–7 June 2024, Lisboa, Portugal
INVESTIGATION OF LENGTH SCALE DEFINITION
INFLUENCE IN LES MODELS
E.Di Lavore , J.Ruano, A.Oliva and F.X.Trias
Heat and Mass Transfer Technological Center, Technical University of Catalonia, C/Colom 11,
08222 Terrassa (Barcelona), enrico.di.lavore@upc.edu
Key words: LES, length-scale, mesh anisotropy, OpenFOAM
Summary.
Large Eddy Simulations (LES) have been increasing popularity due to the decrease of com-
putational power cost. Indeed, engineering applications which have previously involved mainly
RANS, are now shifting to LES. Many eddy-viscosity models have been developed during the
last decades. They all share a general structure, which derive from a dimensional analysis:
νSGS = (Cm∆)
2D(Ū), where in order appear the model constant, the length scale and the dif-
ferential operator underlying the model. The aim of this paper is to investigate the influence of
length scale ∆ definition on highly anisotropic grids, because most of the research has focused
mainly on the model constant and differential operator roles. The main length scale definitions
that will be compared are: i) the most popular is ∆vol = V
1
3 ii) ∆ω iii) ∆lsq. In order to
do so, we first calibrate the model constant, then we carry out a set of simulations in Open-
FOAM to assess the ∆ definition influence. In order to assess the resilience of the models for
highly anisotropic meshes, these simulations will be carried out also on meshes having control
volumes with high aspect ratios. Results for standard test cases such as a decaying Homoge-
neous Isotropic Turbulence (HIT) and a turbulent periodic plane channel will be presented and
compared to reference cases.
1 INTRODUCTION
LES equations are obtained by applying a spatial filter to the Navier-Stokes equations. Usu-
ally a simple ”box filter” is used, though more complex kernels can be used. This results in:{
∂tū+ (ū · ∇)ū = ν∇2ū−∇p̄−∇ · τ
∇ · ū = 0
In order to model the effects of turbulence acting on scales smaller than the filter length ∆ ,
a eddy-viscosity model (e.g. WALE [7], Vreman’s[9], QR[12], Sigma[11], S3PQR [16]) is usu-
ally applied, where the sub-grid scale viscosity is defined as νSGS = (Cm∆)
2D(Ū) . This
is possible because the smallest resolved-scales motions provide information that can be used
to model the largest SGS motions [5] . The aim of this paper is to verify the behavior of an
eddy-viscosity turbulence model while different length-scale definitions are used on progressively
higherly anisotropic meshes.
Di Lavore et al.
1.1 Length-scale definitions
In the LES community, the far more used mesh characteristic length-scale definition is the
cubic root of the volume ∆vol , introduced from Deardorff in 1970 [2]. Though, its performance
is affected from a high aspect-ratio of the volumes, as in unstructured or highly anisotropic
structured meshes. Indeed, this cause the model to be deactivated while the volume shrinks
along one direction. In this situation, defining an optimal ∆ it’s not a trivial task. Nevertheless,
in the DES community, a wider range of ∆ ‘s definitions have been settled, to overcome different
kinds of issues, for example taking into account flow and mesh anisotropy interaction for ∆ω .
Two ∆ definitions are identified for bounding others ∆ behavior: ∆max and ∆min . Indeed,
while refining the mesh in one direction, ∆max is not affected at all, and ∆min is progressively
deactivating the model.
Among the others, ∆ω [10] is meant to progressively deactivate the model while the mesh is
refined in either direction of each plane normal to the vorticity versor. This because the in-plane
velocity components associated to their vorticity component are acting in a progressively finer
mesh, hence making the LES model progressively less necessary in that plane. Instead, ∆lsq [17]
it’s based on a Taylor series approximation of the gradient model. The latter is characterized
by keeping the influence of ∆ as three separate components, in the form of diagonal of a tensor.
In the next section, the HIT case setup will be presented. Then the results in the form of
energy density spectra will compared for different models and degrees of anisotropies. It will
follow the Periodic Plane Channel case, where average velocities and Reynolds stresses will be
plotted. Finally, some concluding remarks to summarise the ∆ definitions behaviors.
2 TEST CASES
2.1 HIT
The first test case is the Homogeneous Isotropic Turbulence (HIT) decay in a periodic cube.
For this case we have used as reference the results from the historical paper by Comte-Bellot
and Corrsin (CBC) [3].
The comparison between the time-evolving HIT and the spatial decay of stationary turbulence
in a wind-tunnel, it’s relying on Taylor’s hypothesis. This latter leads to the assumption that
the spatial correlations of the turbulent velocity components correspond well enough to their
self-correlations in time i.e. (using CBC notation) R11(r, 0, 0; t0, 0) ≈ R11(0, 0, 0; t0, r/Ū) , which
have been verified from Favre et alia [1].
The non-dimensionalization procedure preserves the Taylor micro-scale Reynolds number
Reλ =
urms·λ
ν = 71.6, with a cubic box side length Lbox = 0.09 · 2π, a urms = 0.222 and a
viscosity ν = 10−5 .
The mesh anisotropy is increased by progressively refining the mesh along z-axis. We started
with a uniform mesh of 643 volumes, with the following mesh refinement: 64 × 64 × Nz with
Nz ∈ {64, 128, 256, 512, 1024}.
Regarding the turbulence model, we have chosen Smagorinsky for its simplicity, due to the ab-
sence of walls in this case; whereas for the time integration, a backward 2ndorder explicit Euler
is used , with a constant ∆t = 10−3.
2
Di Lavore et al.
HIT Results
The resolved kinetic energy density spectrum are computed at dimensionless times tU0/M =
42 ; 98 ; 171 and compared with the results from CBC [3]. In Figure 1 are summarized the
outcomes of the HIT case. In each plot is presented the temporal decay of the energy density
spectra, where time correspondence is done by keeping the dimensionless similarity with the
reference.
(a) (b)
(c) (d)
Figure 1: Energy density spectra
In the plot 1(a) is illustrated the spectrum of the main case run in OpenFOAM i.e. on a
isotropic mesh using Smagorinsky model with ∆vol . A consistent energy loss in the energy-
containing range can be noticed, but overall the Kolmogorov 53 law is fulfilled in the inertial
subrange. Anyway, the model produces results (red dash-dotted) in agreements with the refer-
ence (always represented as thinner black dotted) [3].
3
Di Lavore et al.
In the plot 1(b) , keeping the reference in the background, we still show the spectrum at
the last of the considered time-steps, run on isotropic mesh with ∆vol (red dash-dotted). Then
we compare the latter, at first with the same case run without LES model (blue continuous)
to verify that the model is acting properly by modelling the SGS scales and dissipating the
non resolved modes. Indeed, νsgs is properly damping the turbulent kinetic energy pile-up that
can be seen instead in the model-less case. We also show the ∆max (thick green dashed) and
∆min (thin green dashed) spectra from the most anisotropic case (64x64x1024). Intuitively,
∆min is showing a dissipating behavior similar to the model-less case, because the z-axis grid
spacing it’s so small that is almost deactivating the model. Instead, all ∆ ’s behaves equally
to the ∆vol on isotropic mesh, because they are returning the same length-scale, hence same νsgs.
In the plot 1(c) we show the pencil of spectra generated by increasing the degree of anisotropy
in the case of ∆vol . Here the isotropic case spectrum (blue dash-dotted) and the highly
anisotropic case spectrum with Nz = 1024 (green dashed) upper/lower bound the results: the
isotropic case upper-bounds in the inertial subrange and lower-bounds in the dissipation range
; increasing anisotropy leads to an over dissipation in the inertial sub-range and an under-
dissipation in the dissipation range due to the non realistic turbulence seen in the plot 1(b).
This leads graphically to the pivoted pencil visible in plot 1(c).
Finally in plot 1(d) we do the most relevant comparison: are plotted the spectra of all the
used models in the higly anisotropic case. We have already seen the behaviors of ∆min and
∆max , and here we see that they upper and lower bound the results of all the other length-
scale definitions. ∆vol , even if in a minor extent, still suffers model deactivation due to high
anisotropy, because the volume is simply reduced by the refinement. The models that show to
perform the best are ∆ω and ∆lsq , which substantially agree with the reference [3].
2.2 Periodic Plane Channel
The second test case is the periodic plane channel where the flow’s turbulence developement
is guaranteed by applying periodic boundary condition at stream-normal faces while the station-
arity is enforced by a volume force in the stream-wise (x) direction. The reference for this case
is the DNS from Kim, Moin and Moser [4] . We have used a channel size (20π, 2, π) . Viscosity
is set in order to keep the viscous Reynolds number Reτ = 180 as in the reference. Due to
the wall-normal (y-axis) mesh inhomogeneity, we have refined the mesh in either stream-wise or
span-wise directions. Wall-normal spacing have been kept the same as the initial mesh, which
is a homogeneous 323 except along wall-normal direction, where a hyperbolic tangent mapping
is applied to resolve the boundary layers. The meshes are respectively Nxx32x32 and 32x32xNz
with Nx, Nz ∈ {64, 128, 256, 512, 1024} .
For this case we have used the turbulence model WALE [7] that properly deactivates the model’s
effect close to the walls. For time integration, we have used a CFL-adaptive time-stepping, with
Cmax = 0.6 on a 3
rd order Runge-Kutta scheme implemented in the solver RKsymFoam, which
has implemented a symmetry-preserving spatial discretization. [8, 13, 18]
4
Di Lavore et al.
Periodic Plane Channel Results
In order to compare the run cases with the reference [4], we plot in both the homogeneous-
ensemble average velocity and the Reynolds stresses.
(a) (b)
Figure 2: Ensemble-average stream-wise velocity and trace of Reynolds stress tensor in 323 mesh
In figure 2 are shown the behaviors of different ∆ definitions on the same 323 mesh (which is
not isotropic due to the channel aspect ratio). Again, ∆max and ∆min are bounding the others
∆ definitions. Anyway, it can be seen that even with ∆min , the model is over-dissipating
turbulence, which leads to a higher bulk velocity due to a lack of transport of momentum from
the walls towards the central plane of the channel.
In figure 3 is shown how the mesh refinement along the x-axis have reduced this over-
dissipation effect by reducing the pencil-shape volume’s length such that the volume become
a x-axis pancake-shape one. This is considered a natural behavior since in finer meshes, turbu-
lence models become progressively less necessary, up to a resolution that allows for DNS.
Also, this refinement led all the ∆ definitions to converge to the same average velocity profile,
which still present an over-dissipation of turbulence, hence making this discrepancy independent
of LES model.
5
Di Lavore et al.
(a) (b)
Figure 3: Ensemble-average stream-wise velocity and trace of Reynolds stress tensor in
512x32x32 mesh
As it can be seen in figure 4, the mesh refinement impact consistently the behavior of each
∆ definition. We can see in plot 4(b) that the refinement makes ∆ω to converge closer to the
actual result, even if the limit curve it’s still not the reference from [4]. Anyway, it shows
an asymptotical behavior to a steady value. Instead, ∆vol in plot 4(a) shows a less significant
dependency on the mesh. Indeed the mesh refinement don’t show a monotonic approach to
a limit curve, but the curves’ pencil shows a consistent behavior compared to the asymptotic
value of ∆ω . The ∆ definition which is showing the most resilience to mesh anisotropy is the
∆lsq [17], which results overlap while refining the mesh along x-axis.
6
Di Lavore et al.
(a) ∆vol (b) ∆ω
(c) ∆lsq
Figure 4: Ensemble-average stream-wise velocity: length-scale definition influence on x-
refinement
3 CONCLUSIONS
In this work we have tested different length-scale definitions on a range of mesh anisotropies
in two cases, namely the HIT and the plane channel. We have seen how different ∆ definitions
have shown quite different resilience properties with respect to mesh refinement. In particular,
∆lsq [17] have shown to be the least affected from the anisotropies among the tested length-scales,
which is in line with the fact that it has been designed with this purpose. Further research can
be done with other lenght-scale definitions as ∆SLA [15], which is a modification of ∆̃ω [14]
suited to trigger Kelvin-Helmoltz instability in the shear layer, or ∆Sco [6] which is one of the
corrections of Deardorff’s ∆vol [2] required for highly anisotropic meshes as for near wall turbu-
lence and shear layers.
7
Di Lavore et al.
Acknowledgments. This work is supported by SIMEX project (PID2022-142174OBI00) of
Ministerio de Ciencia e Innovacion and the RETOtwin project (PDC2021-120970-I00) ofMinis-
erio de Economia y Competitividad, Spain. E.D.L is supported by Programa Investigo (2023
INV-2 00044) of the Catalan Agency for Management of University and Research Grants (AGAUR).
Calculations were carried out on JFF2 computer cluster at CTTC. The authors thankfully ac-
knowledge these institutions.
REFERENCES
[1] A.Favre, J.Gaviglio, R.Dumas, Appareils de mesures de la correlation dans le temps et
l’espace. Quelques mesures de correlation dans le temps et l’espace en soderie. Proc. 8th
Int. Gong. for Appl. Mech., Istanbul, 304-314, 314-324, 1952
[2] J. W. Deardorff, Numerical study of three-dimensional turbulent channel flow at large
Reynolds numbers, J. Fluid Mech. 41, 453–480, 1970
[3] G. Comte-Bellot and S. Corrsin, Simple Eulerian time correlation of full- and narrow-band
velocity signals in grid-generated, isotropic turbulence, J. Fluid Mech. 48, 273–337, 1971
[4] J. Kim, P. Moin, and R. D. Moser, Turbulence statistics in fully developed channel flow at
low Reynolds number, J. Fluid Mech. 177, 133 , 1987
[5] Germano, M., Piomelli, U., Moin, P. and Cabot, W.H. , A dynamic subgrid-scale eddy
viscosity model. Physics of Fluids A: Fluid Dynamics, 3(7), pp.1760-1765, 1991
[6] A. Scotti, C. Meneveau, and D. K. Lilly, Generalized Smagorinsky model for anisotropic
grids, Phys. Fluids A 5(9), 2306–2308 , 1993
[7] Nicoud, F., Ducros, F. Subgrid-Scale Stress Modelling Based on the Square of the Velocity
Gradient Tensor. Flow, Turbulence and Combustion 62, 183–200, 1999
[8] Verstappen, R.W.C.P. and Veldman, A.E.P., Symmetry-preserving discretization of turbu-
lent flow, Journal of Computational Physics, 187(1), pp.343-368, 2003
[9] W. Vreman, An eddy-viscosity subgrid-scale model for turbulent shear flow: Algebraic
theory and applications, Phys. Fluids 16(10), 3670–3681, 2004
[10] Deck, S. Recent improvements in the Zonal Detached Eddy Simulation (ZDES) formulation.
Theor. Comput. Fluid Dyn. 26, 523–550, 2012
[11] Nicoud, Franck, Hubert Baya Toda, Olivier Cabrit, Sanjeeb Bose, and Jungil Lee, Using
singular values to build a subgrid-scale model for large eddy simulations, Physics of fluids
23, no. 8, 2011
[12] R. Verstappen, When does eddy viscosity damp subfilter scales sufficiently?, J. Sci. Comput.
49(1), 94–110, 2011
8
Di Lavore et al.
[13] F. X. Trias, O. Lehmkuhl, A. Oliva, C. D. Perez-Segarra, R. W. Verstappen, Symmetry-
preserving discretization of Navier–Stokes equations on collocated unstructured grids, Jour-
nal of Computational Physics 258 246–267, 2014
[14] C. Mockett, M. Fuchs, A. Garbaruk, M. Shur, P. Spalart, M. Strelets, F. Thiele, and A.
Travin, Two non-zonal approaches to accelerate RANS to LES transition of free shear layers
in DES, in Progress in Hybrid RANS-LES Modelling, Notes on Numerical Fluid Mechanics
and Multidisciplinary Design Vol. 130, edited by S. Girimaji, W. Haase, S.H. Peng, and D.
Schwamborn, pp. 187–201, 2015
[15] Shur, M.L., Spalart, P.R., Strelets, M.K. et al. An Enhanced Version of DES with Rapid
Transition from RANS to LES in Separated Flows. Flow Turbulence Combust 95, 709–737,
2015
[16] Trias, F. X., David Folch, Andrei Gorobets, and A. Oliva, Building proper invariants for
eddy-viscosity subgrid-scale models, Physics of Fluids 27, no. 6, 2015
[17] F.X.Trias, A.Gorobets, M.H.Silvis, R.W.C.P.Verstappen, A.Oliva. A new subgrid charac-
teristic length for turbulence simulations on anisotropic grids. Physics of Fluids 29, 115109,
2017.
[18] E.M.J. Komen, J.A. Hopman, E.M.A. Frederix, F.X. Trias, R.W.C.P. Verstappen, A
symmetry-preserving second-order time-accurate piso-based method Computers & Fluids,
225: 104979, 2021
9


Investigation of length scale definition influence in LES models

Hua, Y.

A Novel Deformable Limit Analysis Modeling for the Masonry Arch or Vault Structures

Chiozzi, A.

Nale, A.

Senhora, F.

Paulino, G.

A multi-material topology optimization approach for large-scale additive manufacturing of non-periodic architected structures

SIDHARTH, P.

Rao, N.

The 9th European Congress on Computational Methods in Applied Sciences and Engineering
ECCOMAS Congress 2024
3–7 June 2024, Lisboa, Portugal
AN OPEN-SOURCE PHASE FIELD MODEL
FOR FRACTURE IN ORTHOTROPIC FGM
WITH ADAPTIVE MESH REFINEMENT
P.C.SIDHARTH1 AND B.N.RAO1
1 Department of Civil Engineering
Indian Institute of Technology Madras (IITM)
Chennai, Tamilnadu, 600036, India
e-mail: sidharthpccalicut@gmail.com
Key words: Adaptive mesh refinement, fracture, Functionally graded material, Phase-field,
MOOSE
Summary. Fracture prediction in Orthotropic Functionally Graded Materials (FGM) poses
a formidable challenge due to the intricate interplay of spatially varying material properties.
Orthotropic FGMs, characterized by distinct material behaviors along different axes, introduce
heightened complexity in modeling, particularly when addressing fracture phenomena [1]. The
need for accurate and computationally efficient fracture models becomes crucial for understand-
ing and predicting the structural response of such materials. In this context, phase field modeling
emerges as a highly feasible approach, offering an effective means to simulate crack initiation,
propagation, and their intricate interactions in a numerically stable manner. This paper details
implementing a phase field model for fracture prediction in Orthotropic Functionally Graded
Materials (FGM) within the Multiphysics Object-Oriented Simulation Environment (MOOSE),
an open-source finite element framework [2]. Leveraging MOOSE’s flexibility, the study incor-
porates a hybrid split of strain energy proposed by Ambati and introduces an adaptive mesh
refinement scheme to enhance computational efficiency. The primary emphasis is on the metic-
ulous implementation aspects, showcasing how MOOSE accommodates the complexities of or-
thotropic materials. Through reproducing crack paths from literature using in-house codes, the
implemented model’s reliability is demonstrated, validating its accuracy. The innovative use of
adaptive mesh refinement not only reduces computational costs significantly but also highlights
the efficiency gains achieved in simulating fracture in orthotropic FGMs. The paper underscores
the advantages of utilizing MOOSE as an open-source finite element framework for transparency,
collaboration, and accessibility within the research community, showcasing its role in advancing
accurate and efficient fracture predictions.
1 INTRODUCTION
Functionally Graded Materials (FGMs) are advanced engineered materials characterized by
a gradual variation in composition and properties, introduced by Japanese researchers in the
1980s to overcome the limitations of traditional homogeneous materials. This graded structure
allows for a tailored distribution of mechanical, thermal, and electrical properties, optimizing
performance in specific applications by reducing stress concentrations and enhancing structural
integrity. FGMs are used across various industries, including aerospace, biomedical, and en-
P.C.Sidharth and B.N.Rao
ergy sectors, with fabrication methods like powder metallurgy, sol-gel processes, and additive
manufacturing offering unique advantages. Their use in biomedical implants, particularly in
orthopedic applications, is notable for mimicking bone properties. However, FGMs present
challenges in fracture behavior, particularly when anisotropy is introduced through manufac-
turing processes. The phase-field modeling (PFM) approach has emerged as a powerful tool
for simulating fracture phenomena in FGMs, addressing the complexities of anisotropic fracture
behavior. Recent advancements in PFM, including the use of exponential finite element (EFE)
shape functions, have improved computational efficiency and accuracy in modeling crack prop-
agation in FGMs. This paper explores these advancements, comparing numerical simulations
with experimental data to enhance the application of FGMs in engineering.
2 FUNCTIONALLY GRADED MATERIAL (FGM) PLATES
This study focuses on two-dimensional FGM plates with property gradation along the plate’s
length. The Voigt rule of mixtures is applied to homogenize the graded properties across the
FGM plate, with property variations given by:
P (x) = P1 + (P2 − P1)Vx/y (1)
𝑥
𝑦
𝑙𝑥
𝑙𝑦
𝑃1
𝑃2
Figure 1: Model geometry and coordinate system used in Equation 1.
Here, P denotes a material property such as Young’s modulus (E), Poisson’s ratio (ν), or
the critical energy release rate (Gc). Vx/y is the volume fraction of the initial constituent, with
gradation directions defined by:
Vx =
(
1
2
− x
lx
)p
, Vy =
(
1
2
− y
ly
)p
(2)
2
P.C.Sidharth and B.N.Rao
2.1 Orthotropic Functionally Graded Materials
Orthotropic materials display distinct mechanical properties along three mutually perpendic-
ular axes, making them suitable for high-performance structures in aerospace, automotive, and
civil engineering applications. Understanding their fracture mechanics is crucial for ensuring
structural integrity.
3 PHASE-FIELD FRACTURE MODEL FOR ORTHOTROPIC FGMs
Figure 2: Schematic representation of an orthotropic FGM domain with a geometric disconti-
nuity in the PFM framework.
In phase-field models, cracks are represented by a phase-field parameter, which transitions
from 0 (intact) to 1 (cracked). The total energy, Π, in an isotropic system is:
Π =
∫
Ω
ψ(ε) dV +
∫
Ω
Gcγ(ϕ,∇ϕ) dV (3)
Here, Gc is the critical energy release rate, ψ is the bulk energy, and ϕ is the phase-field
order parameter. The crack surface regularization, based on Bourdin’s variational formulation,
is expressed as:
γ(ϕ,∇ϕ) = Gc
(
1
2ℓ
ϕ2 +
ℓ
2
|∇ϕ|2
)
(4)
The bulk strain energy, incorporating a degradation term for stiffness reduction, is:
Πb =
∫
Ω
[(1− ϕ)2 + η]ψ0(ε(u)) (5)
ψ0 =
1
2
εTCε (6)
For orthotropic FGMs, the Cauchy stress tensor σ is defined as:
σ = [(1− ϕ)2 + kp]D(x, θ)ε (7)
where kp ensures numerical stability, and D(x, θ) incorporates material orientation. The
linear system is then solved using the Bubnov-Galerkin method, resulting in:
3
P.C.Sidharth and B.N.Rao
Kuuu
h = fuu, Kϕϕϕ
h = fϕ (8)
This study considers a linearly elastic solid with spatially varying toughness Gc(x) undergoing
small strains. The strong form of the governing equations, without inertia forces, is described
by the following system:
∇e · σ = 0 in Ω,
−Gc(x)ℓ∇ϕA∇ϕ+
[
Gc(x)
ℓ
+ 2H+
]
ϕ = 2H+ in Γ,
(9)
where A = I+ β [I− n⊗ n] accounts for crack orientation based on material orientation. Here,
n = {cos θ, sin θ}T , and β serves as a penalty factor, with β = 20 used to constrain crack
propagation perpendicular to the cleavage plane.
The Cauchy stress tensor for an orthotropic functionally graded material is given by:
σ =
[
(1− ϕ)2 + kp
]
D(x, θ)ε, (10)
where kp is a small positive number introduced for numerical stability. The material stiffness
matrix D(x, θ) is derived as:
D(x, θ) = TTQ(x)T, (11)
with the transformation matrix T and tensor Q(x) defined as:
T =
 c2 s2 2css2 c2 −2cs
−cs cs c2 − s2
 , (12)
and
Q(x) =
Q11 Q12 0Q21 Q22 0
0 0 Q66
 , (13)
where the components of Q(x) are computed as:
Q11 =
E1(x)
1− ν12ν21
,
Q22 =
E2(x)
1− ν12ν21
,
Q12 =
ν12E2(x)
1− ν12ν21
=
ν21E1(x)
1− ν12ν21
,
Q66 = G12(x),
ν21 =
E2(x)
E1(x)
ν12.
(14)
4
P.C.Sidharth and B.N.Rao
The system of linear equations is then derived using the standard Bubnov-Galerkin method:
Kuuu
h = fuu,
Kϕϕϕ
h = fϕ,
(15)
where
Kuu = Σ
∫
Ω
[
(1− ϕ)2 + kp
]
BTD(x)BdΩ,
Kϕϕ =
∫
Ω
{[
2H+ +
Gc(x)
ℓ
]
NTN+Gc(x)ℓB
T
ϕBϕ
}
dΩ,
fuu = Σ
∫
Ω
NT t̄ dΩ,
fϕ = Σ
∫
Ω
NT 2H+ dΩ.
(16)
4 FINITE ELEMENT IMPLEMENTATION
The domain is divided into 4-noded quadrilateral elements as part of a discretization strategy
used in the finite element method to solve the phase-field fracture issue. Weak versions of the
equations are obtained from this division, leading to the following construction of residual and
stiffness matrices.
4.1 Discretization scheme
The domain is partitioned into 4-noded quadrilateral elements for the purpose of applying
the finite element technique to phase-field fracture issues. Stiffness and residual matrices are
obtained from the weak versions of the equations. Standard linear finite element (LFE) shape
functions are utilized to approximate displacement, and exponential finite element (EFE) shape
functions are employed to approximate the phase-field parameter. The Newton-Raphson tech-
nique is used to solve the nonlinear equations iteratively. The coupled equations are solved using
Idaho National Laboratories’ state-of-the-art, open-source MOOSE finite element framework.
For increased efficiency and accuracy, MOOSE makes use of the Jacobian-Free Newton-Krylov
(JFNK) solver. It also supports a number of physics modules for intricate simulations, such
as phase-field modeling and fracture analysis. The solution is obtained through a staggered
approach, calculating displacement and phase degrees of freedom alternately in load increments.
The solution is obtained through a staggered approach, calculating displacement and phase
degrees of freedom alternately in load increments.
4.2 Exponential finite element shape functions
Kuhn and Müller [3] gave the 1D analytical solution for a homogeneous bar exposed to
tension, which was found by solving the coupled equations in (16). The phase-field variable ϕ
may be expressed as follows:
exp
(
−|x|
ℓ
)
(17)
5
P.C.Sidharth and B.N.Rao
The exponential nature of the analytical solution serves as an inspiration for the creation of
exponential finite element (EFE) shape functions. To properly capture this exponential behavior,
the authors recommended using EFE shape functions rather than more conventional linear finite
element (LFE) form functions. The 1D exponential finite element (EFE) shape functions are
given by: ξ′ ∈ [−1, 1], which is a function of the natural coordinate in 1D.
N̄ e1 (ξ
′, δ) = 1−
exp
(
−δ(1 + ξ
′)
4
)
− 1
exp
(
−δ
2
)
− 1
N̄ e2 (ξ
′, δ) =
exp
(
−δ(1 + ξ
′)
4
)
− 1
exp
(
−δ
2
)
− 1
(18)
Unlike linear finite element (LFE) shape functions, which are symmetric and independent of
ϕ, proper orientation of EFE form functions with respect to nodal phase-field values is critical
to accuracy. See [4] for further information on the orientation scheme, derivatives of EFE shape
functions, and their two-dimensional representation.
4.3 Automatic orientation of EFE shape functions
The asymmetry of exponential finite element (EFE) shape functions requires alignment with
the evolving crack path, which is unknown until a critical load increment is reached. The
study uses an automatic orientation scheme that starts with linear finite element (LFE) shape
functions. Elements are marked for correction with EFE shape functions based on the sum of the
phase-field variable at their nodes. The orientation of EFE shape functions is then determined
for these elements, and the analysis is restarted from a lower load increment.
5 NUMERICAL EXAMPLES AND DISCUSSION
5.1 Single Edge Notched Orthotropic Rectangular FGM Plate (SEN-FGM)
Traditional fracture analysis has mainly focused on homogeneous orthotropic materials. With
the introduction of functionally graded materials (FGMs), new challenges and opportunities
have emerged. FGMs exhibit a gradual variation in properties, which enhances their ability
to withstand thermal and mechanical stresses. This study investigates the fracture behavior
of orthotropic FGM plates, focusing on unique crack propagation mechanisms due to material
gradients. The study simplifies the complex variation of material properties by assuming an
exponential gradation:
E1(x) = E1 · eαx
E2(x) = E2 · eβx
G12(x) = G12 · eγx
(19)
Material properties vary proportionally with indices set to 0.2. The finite element mesh
includes approximately 24,000 elements with selective refinement where cracks are anticipated.
The chosen length scale parameter is 0.0015 mm, and an average mesh size of over 0.003 mm
ensures the applicability of EFE shape functions [3].
6
P.C.Sidharth and B.N.Rao
2L
1
2
u
a
x
y
L
Figure 3: Geometry and boundary conditions of SEN-FGM specimen 5.1 showing functionally
graded properties along x-axis
Property Value
Longitudinal stiffness E11 114.8 GPa
Transverse stiffness E22, E33 11.7 GPa
Shear stiffness G12 9.66 GPa
Major Poisson’s ratio ν12 0.21
Critical energy release rate Gf 2.7 N/mm
Table 1: Material properties of SEN-FGM specimen (5.1)
5.1.1 Direction of FGM Gradation
The gradation direction significantly influences crack nucleation and propagation. Examining
the plate with a crack perpendicular to the gradation direction (Y-direction), we find that crack
propagation paths are similar regardless of whether gradation is along the X or Y direction.
However, Figure 5 shows that load responses differ, with increased stiffness for X-direction
gradation compared to Y-direction.
5.1.2 Orientation of the Orthotropic Material Axis
This study assesses the effectiveness of EFE shape functions for predicting fractures in or-
thotropic FGMs. Crack paths for various material axis orientations (30°, 45°, and 60°) are
7
P.C.Sidharth and B.N.Rao
(a) (b) (c)
Figure 4: Crack trajectory predicted by EFE shape functions for SEN-FGM specimen: (a)
θ=30◦, (b) θ=45◦, (c) θ=60◦
(a)
0 . 0 0 0 . 0 1 0 . 0 2 0 . 0 3 0 . 0 4 0 . 0 5 0 . 0 6 0 . 0 7
0
2 0
4 0
6 0
8 0
1 0 0
1 2 0
1 4 0
1 6 0
1 8 0
2 0 0
Pe
ak
 lo
ad
 (N
)
D i s p l a c e m e n t  ( m m )
 P r o p o r t i o n a l  X
 P r o p o r t i o n a l  Y
 N o n - P r o p o r t i o n a l  Y
 N o n - P r o p o r t i o n a l  X
(b)
Figure 5: (a) Comparison of load-displacement response for the orthotropic FGM sample with
different material orientations; gradation of FGM properties in X and Y directions. (b) Com-
parison of load-displacement response for the orthotropic material with different material orien-
tations, SEN-FGM sample.
compared with experimental data from Cahill et al. [5]. Table 2 compares crack propagation
angles with experimental and phase-field model results.
5.1.3 Gradation of FGM Properties
The impact of material property gradation on fracture responses is examined. Two types
of exponential gradation are considered: proportional (with α, β, and γ set to 0.2) and non-
8
P.C.Sidharth and B.N.Rao
Orientation of fibres (θ) 30◦ 45◦ 60◦
Cahill et al. (Experimental [5]) 30 45 60
Hirshikesh et al. (Phase-field [?]) 29.1 43 55
Present study (EFE-Phase-field) 30 43 57
Table 2: Comparison of crack propagation angles with experimental data and phase-field models,
SEN-FGM specimen
proportional (with α=0.5, β=0.4, and γ=0.3). The material properties remain consistent with
previous parameters. Force-displacement responses indicate that non-proportional grading yields
stiffer responses compared to proportional grading, with slight increases in peak load for X-
direction gradation.
6 CONCLUSIONS
In this study, the use of Exponential Finite Element (EFE) shape functions to assess the
mechanical response of an orthotropic functionally graded plate has provided significant insights.
By examining benchmark examples, various factors affecting crack propagation were explored,
including the orientation of material orthotropy, the direction of FGM gradation, and the type of
gradation. The EFE shape functions effectively predicted crack propagation paths in alignment
with material directions, demonstrating their capability to capture the complex behavior of
orthotropic functionally graded materials (FGMs).
REFERENCES
[1] Hirshikesh and Mart́ınez-Pañeda, Emilio and Natarajan, Sundararajan, Defence Technol-
ogy, Vol. 17, pp. 185–195, 2021.
[2] D. Schwen and L.K. Aagesen and J.W. Peterson and M.R. Tonks, Rapid multiphase-field
model development using a modular free energy based approach with automatic differenti-
ation in MOOSE/MARMOT Computational Materials Science, Vol. 132, pp. 36-45, 2017.
[3] C. Kuhn, R. Müller, Journal of Theoretical and Applied Mechanics 49 (2011).
http://www.ptmts.org.pl/jtam/index.php/jtam/article/view/v49n4p1115/227
[4] P.C. Sidharth, B.N. Rao, Engineering Fracture Mechanics 291, 109576 (2023).
https://doi.org/10.1016/j.engfracmech.2023.109576
[5] L.M.A. Cahill, S. Natarajan, S.P.A. Bordas, R.M. O’Higgins, C.T. McCarthy, Compos.
Struct. 107, 119 (2014). https://doi.org/10.1016/j.compstruct.2013.05.039
9


An Open-Source Phase Field Model for Fracture in Orthotropic FGM With Adaptive Mesh Refinement

Mueller, J.

Chiba, K.

Oba, Y.

Identification of losses from engine-airframe interaction for a passenger aircraft through integrated simulation

Not a member yet