SIMULIA LS Med Devices Ebook

2015-10-19

: Ensight Simulia-Ls-Meddevices-Ebook SIMULIA-LS-MedDevices-ebook AbaqusRUM_2015 CEI Houston

Open the PDF directly: View PDF PDF.
Page Count: 34

DownloadSIMULIA-LS-Med Devices-ebook
Open PDF In BrowserView PDF
Optimize Medical Device Reliability
through Multiphysics Simulation

Life Sciences

Industrial Applications of Multiphysics Simulation to
Improve Device Performance, Reduce Time & Costs,
and Support Regulatory Approval
Content
Simulation Now Recognized by FDA as Essential to Medical Device Evaluation

3

SIMULIA Community News

The Living Heart Project: A robust and integrative simulator for human heart function

6

Modeling Hemodynamics with Abaqus/CFD Steady State Solver:
FDA Benchmark Nozzle Model

7

Brian Baillargeon, Nuno Rebelo, and David D. Fox (Dassault Systèmes)
Robert L. Taylor (University of California at Berkeley, Berkeley, California)
Ellen Kuhl (Stanford University, Stanford, California)
European Journal of Mechanics A/Solids

Santanu Chandra, PhD and Richard Swift, PhD, PE (MED Institute Inc.)
Ramji Kamakoti, PhD (Dassault Systèmes)
SIMULIA Community Conference, Providence, Rhode Island

Novo Nordisk Makes Designing Injection Pens a Snap (Fit)

16

Optimization of Surgical Positioning in Total Hip Replacement

18

Computational Study of Cortical Bone Screw Pullout using the
eXtended Finite Element Method (XFEM)

19

SIMULIA Spearheads The Living Heart Project

25

Abdominal Aorta Blood Flow Analysis with Abaqus/CFD

26

Fluid-Structure Interaction (FSI) of a Prosthetic Aortic Valve using Abaqus/Explicit Smoothed
Particle Hydrodynamics

31

Leading innovators in diabetes care use realistic simulation to improve product integrity from design to manufacture
SIMULIA Community News
Jacob Elkins, John Callaghan, Douglas Pedersen, and Thomas Brown
(University of Iowa, Iowa City, Iowa)
SIMULIA Community News

Emer M. Feerick and Patrick McGarry (National University of Ireland, Galway, Ireland)
SIMULIA Community Conference, Providence, Rhode Island
Bringing the Medical Community Together for Improved Patient Care
SIMULIA Community News
Sharat Prasad (Dassault Systèmes)
Abaqus Technology Brief

Cheryl Liu (Dassault Systèmes)
Abaqus Technology Brief

www.3ds.com/simulia

2

Compliance

Simulation Now Recognized
by FDA as Essential to
Medical Device Evaluation

FDA sees increasing numbers of applications that
include simulation
As Life Sciences engineers embrace simulation, they are
achieving increasingly accurate levels of precision when
evaluating device function, including the ability to evaluate
aspects of device performance not possible with bench tests
alone. As a result, the Food and Drug Administration’s (FDA)
Center for Devices and Radiological Health (CDRH) is seeing
a growing number of submissions for medical devices that
include a simulation-data component.

Cheryl Liu, Senior Technical Marketing Specialist,
Life Sciences Industry, SIMULIA
One of the toughest design engineering challenges is making
a medical device that works flawlessly with the human body.
The unique anatomy and physiology of every patient create
physical complexities, and ever-shifting functional parameters,
that must be thoroughly accounted for when producing a
therapeutic product that may need to last a lifetime.

The CDRH is responsible for regulating firms that manufacture,
repackage, relabel, and/or import medical devices sold in the U.S.
The submissions for these therapeutic devices typically contain
data from four types of evaluation models—animal, bench,
computational, and human—to demonstrate a reasonable
assurance of safety and effectiveness. When a company
submits simulation metrics that supplement bench testing, this
can help promote approval by demonstrating both the integrity
of the proposed device and the required realistic device failure
analysis. As the ultimate safety-and-effectiveness regulatory
body between medical device manufacturers and patients, the
FDA recognizes the value of such advancing technologies—and
its own need to stay abreast of them—and has now begun
actively encouraging the use of simulation in device evaluation.

Domestic inpatient procedures involving medical devices—
stents, heart valves, dental implants, spine and joint implants,
surgical tools, blood pumps, endovascular grafts, drug-eluting
devices, and more—totaled 46 million in the U.S. alone in 2006,
according to the Centers for Disease Control and Prevention
(CDC). It’s a global market that is growing along with aging
populations everywhere.
Computer simulation, already widely accepted in many
industries, is increasingly being viewed as an important tool
by medical device companies and their designers. It helps them
visualize what they cannot see, explore the design space more
fully, refine their ideas faster and more accurately—and reduce
expensive prototyping and testing.

However, the FDA has also put the industry on notice that
verification and validation (V&V) must go hand-in-hand with
the use of simulation in applications. The CDRH is looking to
quantify when a computational model is credible enough, and
whether its intended purpose is appropriate for a regulatory
submission. Unclear reporting standards, insufficient data
about geometries and boundary conditions, lack of validation
metrics, incomplete understanding of physiological loads in
the body, and variations in patient populations—any and all
of these uncertainties can impact the relevance of simulation
outputs.

Solid mechanics simulations can help determine proper implant
size, evaluate manufacturing tolerances, compare design
geometries, or consider next-gen devices. Fluid dynamics
can be employed to identify high-shear stresses on blood
vessels, regions of low flow, and potential for blood damage.
And simulation-based product development processes can be
linked in automated workflows, optimizing huge quantities of
design data to provide exquisitely fine-tuned results that are of
particular value for creating patient-specific medical devices.

a)

b)

c)

Figure 1. This example of modeling and simulation of a medical device shows an aortic valve geometry (a), a model of the effect of blood
flow on the valve in a blood vessel (b), and an Abaqus finite element analysis (FEA) of the stress on the valve leaflets during the diastolic
phase (c). This work was performed by SIMULIA in conjunction with the FDA’s Center for Devices and Radiological Health (CDRH).

www.3ds.com/simulia

Previous

3

Next

Contents

Compliance
SIMULIA contributes to advancement of knowledge

The Virtual Physiological Patient

Noticing that a significant proportion of the applications they
have seen in recent years have included simulations with
Abaqus finite element analysis (FEA), the CDRH reached
out to us in 2010 for support in developing their own
internal framework, and in-house expertise, for validating and
regulating industry-submitted simulations.

DBS
Massachusetts
General Hospital

ODEC

SIMULIA presented at the FDA’s 3rd workshop on Computational
Modeling of Medical Devices the same year. We continue to
deliver on-site training courses to FDA reviewers about best
practices in modeling and simulation and to partner with the
FDA on aortic valve model development (see Figure 1). The FDA
has also presented at our SIMULIA Community Conference
and Regional User Meetings. Realizing the importance of
model V&V in early 2011, ASME and FDA launched the V&V 40
subcommittee to develop V&V guidelines for the medical device
industry specifically; we are actively participating, along with
others in the industry and software communities.

Johns Hopkins
University

CABAS

IMCAFS
Georgia Institute
of Technology
Penn Medicine

As one outcome of these efforts, the FDA will publish a
guidance document titled “Reporting Computational
Modeling Studies in Medical Device Regulatory Submissions”
in 2013. Appendices will cover fluid and mass transport,
solid mechanics, electromagnetism, control loops, thermal
transport, and ultrasound. Publication date updates can be
found on the CDRH website at www.fda.gov/MedicalDevices/
DeviceRegulationandGuidance/GuidanceDocuments/
ucm371016.htm

BioEVAL
University of
Pittsburg

CHAP
Dartmouth
Hitchcock
M2S
SENTARA
Penn Medicine

HIPS
Cleveland Clinic
Massachusetts
General Hospital
BARD
Boston Scientific
Abbott
ndc
SIMULIA

As knowledge about the importance of simulation grows,
another priority for the FDA is the creation of a publicly available
‘Virtual Physiological Patient’ of human body computer models
in different disease states (see Figure 2). This is not intended to
be a single model encompassing every function and disease at
once. Rather, the project will comprise a library of verified and
validated submodels and data based on the combined expertise
of those groups in the relevant disciplines, i.e., cardiology,
orthopedics, software, and so forth.

VFAM
IT’IS Foundation

The goal of the Virtual Physiological Patient project is a
shared point of reference that will improve understanding of
model attributes and limitations, and provide discrete models,
data, and simulations validated for regulatory evaluation. Peer
review by experts in academia, government, and industry will
ensure robust V&V and provide periodic assessment. SIMULIA
is contributing expertise to a group that is developing a
computational model for the evaluation of a diseased femoral
artery for stent evaluation.

Previous

University of
Colorado School
of Medicine
The University
of Iowa

ASPECT

The ‘Virtual Patient’ idea is born

www.3ds.com/simulia

CPI/CFD
Cleveland Clinic
Mississippi State
University
Penn State
RIT
University of
Pittsburg

Figure 2. Broad cross-industry collaboration between medical device
manufacturers, academia, and software companies is being harnessed
for the FDA’s Virtual Physiological Patient project.

4

Next

Contents

Compliance
Newly launched public-private partnership
benefits all parties

the function of the living body is critical to every medical-device
developer, and sharing the data that lie at the core of that
understanding can be accomplished without infringing on any
one company’s patents.

Concurrent with the development of the Virtual Physiological
Patient concept, the FDA is reaching outward to device
manufacturers, software providers, and medical professionals
to form a Regulatory Science Public-Private Partnership.
Launched in December of 2012, the partnership is called the
Medical Device Innovation Consortium (MDIC).

The FDA views modeling and simulation as incentives to
innovation that can reduce the time and cost of device design,
assessment, and manufacturing. It is in all our interests—the
medical device industry, the regulatory agency, and software
companies—to collaborate to ensure that the power of simulation
is increasingly utilized to solve the wide range of challenges in
medical device development. We can all agree that the ultimate
goal is the safety and effectiveness of medical devices for every
physician who uses them, and every patient who needs them.

The idea is to create an opportunity for information gathering
in a pre-competitive state, i.e., not device-specific, but diseasespecific. For example, if the heart valve community were
interested in a comprehensive evaluation of the structure and
function of heart valves, costs could be minimized through
nonprofit group funding and participation in the development
of tools and resources for modeling and simulating of a range
of valves. All results would be shared. End-stage renal disease
is another area recently identified by the FDA as a priority.
Industry forums on this topic are already underway.

Special thanks to Dr. Tina Morrison, Dr. Nandini Duraiswamy,
and Dr. Donna Lochner of the FDA for their assistance in
preparing this article.
Read more about how the FDA is promoting innovation in
“High Stakes Balancing Act” in Compass magazine—
www.compassmag.3ds.com.

The medical device industry can only benefit from such
endeavors. Individual device design copyrights certainly need
to be protected, but the tradition of publishing evidencebased research results in order to move the entire body of
medical knowledge forward has resonated in the life sciences
throughout the history of medicine. A deep understanding of

www.3ds.com/simulia

Previous

For More Information
www.fda.gov
www.deviceconsortium.org

5

Next

Contents

Cardiovascular

The Living Heart Project: A robust and integrative simulator
for human heart function
Brian Baillargeon a, Nuno Rebelo a, David D. Fox b, Robert L. Taylor c, Ellen Kuhl d,*
Dassault Systèmes Simulia Corporation, Fremont, CA 94538, USA
Dassault Systèmes Simulia Corporation, Providence, RI 02909, USA
c
Department of Civil and Environmental Engineering, University of California at Berkeley, Berkeley, CA 94720, USA
d
Departments of Mechanical Engineering, Bioengineering, and Cardiothoracic Surgery, Stanford University, Stanford,
CA 94305, USA
a

b

ABSTRACT
The heart is not only our most vital, but also our most complex organ: Precisely controlled by the interplay of electrical and
mechanical fields, it consists of four chambers and four valves, which act in concert to regulate its filling, ejection, and overall pump
function. While numerous computational models exist to study either the electrical or the mechanical response of its individual
chambers, the integrative electro-mechanical response of the whole heart remains poorly understood. Here we present a proofof-concept simulator for a four-chamber human heart model created from computer topography and magnetic resonance images.
We illustrate the governing equations of excitation-contraction coupling and discretize them using a single, unified finite element
environment. To illustrate the basic features of our model, we visualize the electrical potential and the mechanical deformation
across the human heart throughout its cardiac cycle. To compare our simulation against common metrics of cardiac function, we
extract the pressure-volume relationship and show that it agrees well with clinical observations. Our prototype model allows us to
explore and understand the key features, physics, and technologies to create an integrative, predictive model of the living human
heart. Ultimately, our simulator will open opportunities to probe landscapes of clinical parameters, and guide device design and
treatment planning in cardiac diseases such as stenosis, regurgitation, or prolapse of the aortic, pulmonary, tricuspid, or mitral valve.
Read the entire paper at http://dx.doi.org/10.1016/j.euromechsol.2014.04.001.
Mechanical displacement

© 2014 The Authors. Published by Elsevier Masson SAS. This is an open access article under the CC BY-NC-ND license
(http://creativecommons.org/licenses/by-nc-nd/3.0/).

* Corresponding author.
E-mail address: ekuhl@stanford.edu (E. Kuhl).

www.3ds.com/simulia

Previous

6

Next

Contents

Cardiovascular

Modeling Hemodynamics with Abaqus/CFD Steady State
Solver: FDA Benchmark Nozzle Model
Santanu Chandra, PhD and Richard Swift, PhD, PE (MED Institute Inc.)
Ramji Kamakoti, PhD (Dassault Systèmes)
Abstract: Understanding the blood flow dynamics (hemodynamics) and the fluid forces exerted on the blood by implantable medical
devices and predicting blood damage is an intricate part of interventional medical device design. Computational techniques such as
Computational Fluid Dynamics (CFD) are increasingly being used as a tool for describing complex hemodynamics and calculating
the fluid forces such as pressure and shear stress. However, this technique is challenged by the lack of standardized methods for
validation and verification of the results.
In an effort to standardize the process, the FDA made an initiative in 2008 to engage the Medical Device/CFD Community
worldwide to participate in an open challenge that involves CFD computation of blood flow and its experimental validation with
Particle Image Velocimetry (PIV) on a benchmark nozzle model. The goal of the CFD phase of the study was to understand the
limitations of CFD, understand the variability that arises due to choice of software or solver, and the influence of user expertise level
and diverse modeling and meshing techniques. This information would then be employed to identify best practices and define the
critical techniques necessary to achieve credible results. Though a variety of software packages have been used by 28 different
groups from 6 countries, Abaqus/CFD has not been publicly applied in this challenge.
In this study we aim to assess the performance of Abaqus/CFD in modeling hemodynamics using the FDA Benchmark Nozzle model.
The Nozzle model consists of a tube with a straight section, followed by a conical section, a section with reduced tube diameter
and a section with a sudden expansion in tube diameter. The device was designed to include accelerating flow, decelerating flow,
variations in shear stress and velocities, and recirculating flow, all of which may be present in a medical device and relate to blood
damage. Five different flow rates corresponding to fully laminar, transitional and fully turbulent flow were modeled with laminar
and turbulent solvers. In this study we have used pure hex mesh as well as the steady state solver introduced recently to solve
both laminar and turbulent flow problems and compared the results with published experimental data. The hex-dominant mesher
available in the 3DExperience platform was also tested with this model (results not presented here). For solving the turbulent flow at
higher flow rates we tested three models, namely Spalart-Allmaras, k-ε and k-ω models. Results of interest were the axial velocity,
pressure and fluid shear stress along the centerline, axial velocity along a radius of cross section at selected planes, and wall shear
stress.
Overall, Abaqus/CFD results were found to match very well with experiment results and faired competently with other software
results. The k-ω model was found to perform best among the three turbulent models. This detailed study provides valuable insight
into effective strategies for modeling hemodynamics using Abaqus/CFD.
Keywords: CFD, hemodynamics, FDA Nozzle, Sudden expansion, steady state solver, turbulent models, Spalart-Allmaras, RNG k-ε
model, SST k-ω model

1. Introduction
Computational Fluid Dynamics (CFD) techniques are intensely used in the design and analysis of components in heavy industries,
such as the aerospace and automotive industries. Following the same trend this technology is being used in the medical device
industry to design, develop and analyze devices like ventricular assist devices, prosthetic heart valves, stents, blood filters and
hemodialysis catheters. However this technique can not only help in the design process of a device but also can predict the
change in the hemodynamic environment (i.e., the change blood flow dynamics and/or the fluid forces exerted on the device). This
information is crucial as blood flowing through medical devices is subject to hemolysis or thrombosis. Understanding the effect of
hemodynamic forces like shear stress and exposure time on hemolysis has been an area of active research focus for the past decade.
Though crucial understanding has been achieved and several power law based empirical formulas have been modeled, the use of
CFD and blood damage models in predicting medical device safety has not been adequately proven yet. To successfully implement
these techniques for medical device design and/or for regulatory reviews a better understanding of the computational model and
experimental validation must be achieved and steps towards “standardized” practices should be taken.

www.3ds.com/simulia

Previous

7

Next

Contents

Cardiovascular
With an aim to utilize this technology to its fullest and standardize it, FDA initiated a critical path initiative named as” Computational
Fluid Dynamics (CFD)/Blood Damage Project”. The initiative consist of two projects, the first one is CFD/PIV analysis of a Nozzle
model that represents an idealized medical device geometry and the second one is a CFD/PIV and blood damage experiments on a
heart pump model. The goal of the first project, which has been evaluated in this study, was to assess the current state of the art in
CFD analysis in modeling hemodynamics. The experimental PIV analysis was performed in 4 labs, 3 from University research groups
and one from the FDA. The averaged data was used for validation purposes. Twenty eight groups from 8 countries participated
in the CFD challenge and submitted their simulation results. This project was initiated in 2008 and the experimental results were
published in 2010 (Hariharan, 2010). The computational results with experimental validation were published in 2011 (Stewart,
2011).
Abaqus/CFD, being a relatively recent addition to the SIMULIA brand has not been documented in this challenge. Therefore the
need to validate the standard hemodynamics process and assess the performance of Abaqus/CFD was the prime motivation for
this study. The easy accessibility of the experimental results from the FDA repository is an additional reason for choosing this
benchmark model as a validation standard.
The FDA Nozzle model consists of a straight tube, a conical section that connects to a straight throat region of smaller diameter,
and a sudden expansion region where the diameter suddenly increases back to the original tube diameter. The device was designed
to include accelerating flow, decelerating flow, sudden expansion/gradual diffusion, variations in shear stress and velocities, and
recirculating flow, all of which may be related to blood damage in medical devices. For this geometry, flow in one direction is a
sudden expansion problem, whereas flow in the opposite direction is a conical diffusor problem. In this study we concentrated on
the sudden expansion problem. The model was simulated with 5 different flow rates that cover the range of fully laminar to fully
turbulent flow.
Through this study we present the capabilities of Abaqus/CFD as a comprehensive tool offered by SIMULIA, capable of solving for
laminar, transitional and turbulent flow with accuracy. Abaqus/CFD can be used within the Abaqus/CAE environment in all three
stages of analysis - modeling/meshing, solving and post-processing. In version 6.13, Abaqus/CFD has introduced a steady state
solver that is robust and uses a 2nd order accurate SIMPLE based algorithm. Abaqus/CFD has also introduced three turbulent flow
models, the first one is a one equation Spalart-Allmaras model, and the other two are two equation RNG k-ε and SST k-ω models.
The first two models are available in CAE and the last one is available through insertion in the input deck. Details of the modeling
technique followed are discussed in the method section, and the simulation results are presented in the Results and Discussion
section.

2. Methods
Nozzle Geometry
The dimensions of the nozzle geometry were obtained from the published literature (Stewart, 2011). The FDA Nozzle model
consists of a straight tube of diameter (D) 12 mm, a conical section that reduces the tube diameter to 4 mm in the throat region,
and a small straight throat section of 40 mm length which is followed by a sudden increase in tube diameter (D) to 12 mm. The
inlet length and the outlet length were left for the analyst to choose. We have chosen the inlet length to be 10 times D. We have
used a parabolic inlet velocity profile therefore larger inlet length was unnecessary. For assessing the outlet length, simulations
were performed with different outlet length and the velocity profile at the outlet observed. Our objective was to obtain a parabolic
velocity profile at outlet. Results with shorter outlet length resulted in a velocity profile where the peak velocity is eccentric. As
the length was increased, the outlet velocity was fully developed and peak velocity was concentric. Depending on this criterion we
chose the outlet length to be 60D.

www.3ds.com/simulia

Previous

8

Next

Contents

Cardiovascular
The 3D model of the CFD domain as shown in Figure 1a was developed using the solid modeling features in Abaqus/CAE .

Figure 1. (a) FDA Nozzle model geometry with cross sectional surfaces identified by numbers b) final hex mesh of the
model c) view of the radial section of the mesh at inlet

Meshing
Linear Hex elements were used to mesh nozzle model in this study. The model volume was subdivided along the axial direction
and also along the radial direction as shown in in Figure 1b. There were 13 cross sectional surfaces of interest identified for
post processing of results that were used for volume subdivision along the axis. The subdivision in the radial direction helped in
prescribing adequate seeds near the walls. A mesh refinement was performed by changing the global seed from 0.001 to 0.004
along with adequate seeding of radial and axial edges. The final mesh consists of 1.39 million elements. The final mesh pattern is
shown in Figures 1b and 1c.
Another meshing technique was also tested and used for this model. It is the hex-dominant mesher available as part of the R2014
3DExperience platform. It is capable of handling complex geometries in a robust manner, and creating high-quality Hex-dominant
3D meshes. The Hex-dominant 3D mesher includes surface wrapping technology to automatically extract the fluid domain by
creating a clean closed ‘wetted’ surface for subsequent volume meshing. A snapshot of the mesh created for the nozzle geometry is
shown in Figure 2 (a). A closer look of the mesh, shown in Figure 2 (b) shows the smooth transition from a finer mesh (near the wall)
to a coarser mesh away from the wall. The robustness and the time it takes to create high quality meshes, makes this mesher useful
for creating CFD meshes for both simple and highly complex geometries. The results obtained by using this mesh are not included
in this paper, but the results were very comparable to the more traditional hexahedral element mesh created using Abaqus/CAE.

www.3ds.com/simulia

Previous

9

Next

Contents

Cardiovascular

(a)

(b)
Figure 2: Nozzle mesh using R2014 3DExperience Mesher

Materials
Blood is modeled as incompressible fluid with a density of 1056 kg/m3. It is assumed be a Newtonian fluid with dynamic viscosity
of 0.0035 N.s/m2 as per the guidelines of the challenge (Stewart, 2011)

Boundary Conditions
Simulations were performed with 5 different inlet velocities ranging from inlet Reynolds number 167 to 2167 as documented in
(Stewart, 2011).The flow rates were labeled as Q1, Q2, Q3, Q4 and Q5.The average inlet velocity, flow rate, average velocity at
throat and throat Reynolds number are presented in Table 1. Though the inlet Reynolds number indicates the flow to be in a fully
laminar zone at Q1 to a transitional zone at Q5, the throat Reynolds number indicates that the flow type ranged from fully laminar
at Q1 to fully turbulent at Q5 in the sudden expansion region. In this study we have chosen to represent the flow rate Q1 as fully
laminar flow, flow rate Q3 as representative of transitional flow and flow rate Q5 as fully turbulent flow. A parabolic inlet velocity
profile Vz was used to model the inlet flow distribution following Equation 1
											

1.

where Vmax is 2 times the average inlet velocity, r is the radial coordinate , R is the inlet radius and a = Vmax and b = Vmax/R2.
The outlet boundary condition was set to zero pressure and a no-slip boundary condition was applied to the wall surface of the
model
Table 1. Inlet flow rate, velocity, Reynolds number used in the simulation
Inlet
Re

Average
Inlet
Velocity

Flow
Rate

Peak
Velocity
(Parabolic)

Rei

m/s

m3/s

m/s

a

167

0.05

5.22e-6

0.09

Q2

667

0.18

2.08e-5

Q3

1167

0.32

3.65e-5

Q4

1667

0.46

Q5

2167

0.60

Q1

www.3ds.com/simulia

Velocity
at throat

Throat
Re

b

m/s

Ret

0.09

2563.52

0.42

500

0.37

0.37

10234.73

1.66

2000

0.64

0.64

17906.93

2.90

3500

5.21e-5

0.92

0.92

25579.14

4.14

5000

6.77e-5

1.20

1.20

33251.35

5.39

6500

Previous

10

Parabolic Profile

Next

Contents

Cardiovascular
Analysis Settings
A steady state incompressible fluid flow analysis was performed in Abaqus/CFD. The steady state method is a new addition to
Abaqus/CFD version 6.13 and is based on a second order accurate SIMPLE (Semi-implicit Method for Pressure Linked Equation)
based algorithm. Depending on the problem, a proper choice of an under-relaxation factor can significantly increase the convergence
rate of the solution. For this study we chose 0.7 for the Momentum equation and 0.3 for the pressure equations and obtained a
converged solution in 1000-1500 iterations. For the lowest flow rate Q1 the laminar flow model was used. For all other higher flow
rates, we have used the laminar flow model along with three turbulent flow models. Abaqus/CFD offers a Spallart-Allmaras model
and RNG k-ε model in CAE and the SST k-ω model can be accessed through the input deck. The Spallart-Allmaras model is a one
equation Reynolds Averaged Navier Stokes (RANS) turbulent model that solves for the Reynolds’s stress when the Spalart-Allmaras
variable kinematic turbulent viscosity is known. Abaqus/CFD recommends calculating the kinematic turbulent viscosity (ν∼) as 3
to 5 times of the kinematic viscosity (ν) of the fluid. We used a value of 1.3258e-5 Stoke, which is 4 times of kinematic viscosity.
RNG k-ε model and SST k-ω models are two-equation RANS models that solve for Reynold’s stress by including two transport
equations to represent the turbulent properties of the flow, like turbulent kinetic energy(k) and turbulent dissipation (ε) and /or
specific dissipation(ω). The RNG k-ε model is a variant of k-ε models where the Re-Normalization Group (RNG) method is used
to re-normalize the Navier stokes equation. The parameters k and ε can be calculated using Equation 2 if the inlet velocity (u0),
turbulent intensity (I) and turbulent length scale (l) is known.
												

2.

Cµ is an empirical constant whose value is 0.085. The SST k-ω model is a two equation eddy viscosity model where a Shear Stress
Transport (SST) formulation is combined with a standard k-ω model and represents the “best of both” of the models. The variables
for the transport equation are available in Abaqus/CFD whereas the parameters k, and ω can be calculated if the inlet velocity,
turbulent length scale and turbulent intensity are known using Equation 3.
													

3.

For this study the turbulent intensity and length scales were not available so we assumed the values to be 0.057 and 0.84 mm,
given that turbulence intensity (I ) at core of a fully developed pipe flow can be approximated by
and the
turbulent length scale at inlet can be approximated by 0.07*L, where L is the characteristic length, in this case it is the inlet diameter
of the model. Table 2 summarizes the parameter values used for 4 different flow rates (Q2, Q3, Q4 and Q5).
Table 2. Turbulence parameters used for RNG k-ε model and SST k-ω model, for flow rates corresponding to throat Reynolds number

 

u (inlet)

k

eps (ε)

omega(ω)

Turbulence
Intensity

length scale

Q2:Re = 2000

0.1839

1.64E-04

3.96E-04

28.31

0.057

8.40E-04

Q3:Re = 3500

0.3218

5.04E-04

2.12E-03

49.53

0.057

8.40E-04

Q4:Re = 5000

0.4606

1.03E-03

6.23E-03

70.89

0.057

8.40E-04

Q5:Re = 6500

0.5986

1.74E-03

1.36E-02

92.13

0.057

8.40E-04

0

Simulations were performed on a desktop computer equipped with an Intel CoreTM i7 CPU and 16 GB of RAM. Simulation times
were between 5 hours to 12 hours.

www.3ds.com/simulia

Previous

11

Next

Contents

Cardiovascular
3. Results & Discussion
3.1 Centerline Velocity and Pressure
Results were post-processed in the Abaqus/CFD visualization module. Snapshots of the velocity magnitude plot from three different
flow rates are shown in Figure 3. The snapshot shows the jet length in the sudden expansion zone for laminar flow model with
flowrate Q1 and turbulent flow model with flowrates Q3 and Q5. It should be noted here that there are some instabilities seen for
the Laminar case near the outlet. This is an outcome of the steady state plot shown at an instant prior to fully converged solution.

Figure 3 : Velocity magnitude plot obtained from laminar model for flow rate Q1 and turbulent model (SST k-ω model)
are plotted for flow rates Q3 and Q5.

A direct comparison of the Abaqus/CFD results with the experimental data obtained from FDA repository is presented in Figure 4.

Figure 4 : Centerline velocity in the axial direction, obtained from laminar and turbulent simulations, are plotted for flow rates (a) Q1 , (b) Q3 and (c) Q5.
Similarly, wall pressure, represented by the centerline pressure in the axial direction is plotted for flow rates (d) Q1 (e) Q3 and (f) Q5. Flow rates Q1, Q3
and Q5 are representative of fully laminar flow, transitional flow and fully turbulent flow.

www.3ds.com/simulia

Previous

12

Next

Contents

Cardiovascular
Centerline velocity in the axial direction and pressure data were recorded using the ‘path’ feature in the visualization module and
plotted in Figure 4. The data obtained were for three different flow rates that are representative of fully laminar (Q1), transitional
(Q3) and fully turbulent (Q5) flow types. Due to very little radial pressure gradient, wall pressure is represented by the centerline
pressure. The experimental data presented in Figure 4 are an average of 3 different experimental data sets obtained from the FDA
repository [3]. The experimental and computational pressure values are referenced with respect to the pressure at the location of
the sudden expansion (z=0.0). Our simulation results matched very closely with the experimental data for laminar, transitional
and turbulent flow. As the flow through the nozzle before the sudden expansion (z=0) region is truly laminar, the centerline axial
velocity characterized by the laminar flow model was in excellent agreement with the experimental data for all three flow rates.
Flow beyond the sudden expansion zone is complicated and none of the models were an exact match with the experimental data.
Similar behavior was reported by the other contributors to the challenge and this was comprehensively discussed in the FDA
publication (Stewart, 2011). At the sudden expansion zone the laminar flow model did well for lower flow rates and demonstrated
fluctuations for transitional and fully turbulent flow rates. Here, the SST k-ω model and Spallart-Allmaras model achieved the best
agreement with the experimental data.
Figure 4(d-f) represents the pressure plots. Clearly pressure from laminar flow model for flow rate Q1 (Figure 3d) did not match with
the experimental data. As reported (Stewart, 2011) results from other computational studies did not match with the experiment
either. The errors are attributed to experimental error in pressure measurement, and were partly due to wrong reference scaling and
were partly for not using differential pressure transducers and were discussed in the literature (Stewart, 2011). For higher flow rates
(Q2 & Q3), the largest discrepancy in pressure was observed for the RNG k-ε model. The discrepancy in the pressure characteristics
by the RNG k-ε model can be attributed to its inability to characterize turbulence in multiple different length scales as the eddy
viscosity is determined from a single selected turbulence length scale. It is also well known that the RNG k-ε does not perform well
for transitional Reynolds number, which we believe is the case in cases Q2 and Q3. In addition, the impact of wall mesh resolution
requirement (y+ near the wall) could also be a contributing factor to the deviation from experimental results. Future work will involve
using the newly introduced realizable k-ε turbulence model, which uses a hybrid wall functions approach thereby alleviating the
need for creating meshes with specific near wall resolution (any y+ value near the wall) as opposed to the RNG k-epsilon, which had
a stringent near wall mesh requirement (y+>20 near the wall). Similar behavior of the k-ε model was observed by users of other
software as well (Stewart, 2011). Though none of the turbulent flow models were an exact match with the experimental data at
each location, the Spalart-Allmaras and SST k-ω models did a superior job in characterizing the axial velocity and pressure before
and beyond the sudden expansion zone, whereas the laminar model results deviated slightly in the sudden expansion region. The
next section describes the axial velocity profile in the radial direction at specific cross sectional surfaces along the axial length of the
model.

3.1. Axial velocity profile at selected cross sections

Figure 5 : a) Axial velocity in radial direction at specific cross sectional surfaces (Surf-1,Surf-3,Suf-6,Surf-8 and Surf-13)
are obtained from simulation and plotted along each row for three different flow rates (b) Q1 (c) Q3 and (d) Q5.

www.3ds.com/simulia

Previous

13

Next

Contents

Cardiovascular
The nozzle model was subdivided in such a way that 13 cross sectional surfaces were created as shown in Figure 1a. Out of the
13 surfaces, results from only the following surfaces of interest were presented in this study: Surf-1 is closest to inlet, Surf-3 is at
the conical section, Surf-6 is at the throat region, Surf-8 is in the sudden expansion zone and Surf 13 is the farthest away from the
sudden expansion region. Axial velocities along the radial direction at these cross-sections were recorded for three strategic flow
rates (Q1,Q3 and Q5) and presented in Figure 5 (a, b and c) respectively. The experimental data presented here is from one set of
experiments obtained from the FDA repository [3] that best matches the average experimental data. As expected, the velocity profile
from the laminar flow model matched very well before the sudden expansion zone for all three flow rates. As flow rate increases,
the turbulent flow models did a better job in matching the experiment results at and beyond the sudden expansion region, as
demonstrated clearly by the Surf-13 results described in Figure 5b and Figure 5c.

3.2. Radial shear rate profile at selected cross sections

Figure 6 : Shear Rate in radial direction at specific cross sectional surfaces (Surf-1, Surf-3, Suf-6, Surf-8 and Surf-13), obtained
from simulation and plotted along each row for three different flow rates (b) Q1 (c) Q3 and (d) Q5.

We have further post processed hemodynamic parameters such as shear rate at selected cross sectional surfaces (Figure 6). Shear
Rate is the ratio of velocity and distance and fluid shear stress can be calculated by multiplying shear rate and dynamic viscosity.
Results from the laminar flow models matched the shear rate characteristics very well for lower flow rate Q1 and transitional flow
rate Q3. With flow rate Q3, a larger discrepancy in shear rate was observed at Surface 10 and none of the models were effective in
capturing the true behavior of the flow observed through experiment. However, as the flow rate increased to Q5, the turbulent flow
models were successful in capturing the flow behavior at this location. As flow rate Q5 is still in a transitional zone in the inlet tube,
the laminar flow model did a better job in capturing the flow behavior at this location. We have further post-processed the wall
shear stress (data not presented here) and found very good correlation with the experimental data. Overall, the simulation results
are in very close agreement with the experimental data and were satisfactory, particularly in light of the wide variability of data
presented in the challenge.

4. Conclusion
Abaqus/CFD was used to model, simulate and visualize the flow through a benchmark FDA nozzle model and the simulation results
were validated by experimental data obtained from the FDA repository and presented in this technology brief. Abaqus/CFD now
has the steady state feature, three turbulent models (Spalart-Allmaras model, RNG k-ε and SST k-ω) and a laminar flow model. All

www.3ds.com/simulia

Previous

14

Next

Contents

Cardiovascular
of these features were put on test for modeling hemodynamics. The results illustrate the challenging nature of this benchmark and
highlight areas where particular user attention is required. This study also revealed that it can be challenging to predict which model
type is appropriate for any general medical device geometry for a given flow rate, as laminar, transitional and turbulent flow regions
may co-exist in different locations at the same time. This also emphasizes the fact that CFD is not a “push button” technology, and
user skill is an important factor in achieving accurate results, while the necessary details in modeling, meshing and analysis set-up
play a crucial role in achieving higher levels of accuracy.

5. References
1.	 Abaqus User’s Manual, Version 6.10-1, Dassault Systèmes Simulia Corp., Providence, RI.
2.	 FDA Critical Path Initiative-Computational Round Robin #1 (https://fdacfd.nci.nih.gov/interlab_study_1_nozzle)
3.	 Hariharan, P. et al., Multi-laboratory Particle Image Velocimetry Analysis of the FDA Benchmark Nozzle Model to Support
Validation of Computational Fluid Dynamics Simulations, Journal of Biomechanical Engineering, April 2011, Vol. 133 / 041002-1
4.	 Stewart, S. et al., Assessment of CFD Performance in Simulations of an Idealized Medical Device: Results of FDA’s First
Computational Inter-laboratory Study, Cardiovascular Engineering and Technology, Vol. 3, No. 2, June 2012 ( 2012) pp. 139–160

www.3ds.com/simulia

Previous

15

Next

Contents

Drug Delivery

Novo Nordisk Makes
Designing Injection Pens
a Snap (Fit)

Injection typically involves twisting a short needle onto the
pen, turning a dial to the required dose, and pushing a button
to deliver the medication under the skin. After a given number
of doses are injected, the cartridge is exchanged for a new one
(with a durable device) or discarded (with disposable pens).
Audible clicks that occur at key stages of this procedure reassure
the patient that they are engaging the device correctly at each
step. It looks pretty easy. A one-minute video of a woman
checking her blood sugar and then using a NovoPen to inject
insulin is available here: http://www.novonordisk.com/press/
broadcastroom/default.asp. But every one of those reassuring
clicks represents a challenge that has been overcome by the
engineers who created the pens. So do the clicks the patient
never hears: those that occur as the pen parts are assembled in
the factory before use.

Leading innovators in diabetes care use
realistic simulation to improve product
integrity from design to manufacture

“Parts that click into place with ‘snap fit’ instead of screw
connectors are very efficient to assemble within mass
production,” says Torben Strøm Hansen, principal scientist in
the Device R&D division of Novo Nordisk, near Copenhagen,
Denmark. “Snap fit is the commonly used way to connect parts
in our device mechanisms, and it also signals reliability when the
internal components have optimal connections that don’t rattle.
It’s very efficient when designed correctly.”
Getting those designs correct from the start is the task that
Hansen and his Mechanical Analysis team focus on in close
collaboration with Novo Nordisk’s mechanical designers. “Even
though an injection pen is not that big, there are a lot of fine
details in its design,” he says. Whatever the configuration of
device, the plastic-polymer components must withstand the
rigors of both manufacturing and patient use, performing as
required at different temperatures and loads.

Many medical conditions can be treated with tablets, but others
require injections under the skin in order for therapeutic drugs
to reach the bloodstream. In the case of insulin administration
for diabetes treatment, patients need to self-inject the drug
daily.
Making those injections easy and safe is of prime importance
for Novo Nordisk, the Danish company that has been a world
leader in the production of insulin ever since it was discovered
by Canadian scientists in the 1920s. The company innovated
beyond standard syringe technology to produce the world’s
first patient-friendly self-injection system, the NovoPen, some
25 years ago.
With more than 350 million diabetics worldwide—8.3 percent
of the global population and growing, according to the
International Diabetes Federation—demand for insulin pens will
likely remain strong into the foreseeable future. Since effective
control of the disease is dependent on consistent use of the
drug, these delivery systems need to be portable, easy-to-use,
reliable, and even resistant to minor misuse by patients.

Small medical device, big design task
An insulin pen may be small, but it is a precision instrument
with a number of complex parts that must work in perfect
tandem. Some pens are durable, containing a replaceable
drug cartridge, while other disposable ones come pre-filled
with the drug.

www.3ds.com/simulia

Previous

(Top) CAD image of diabetes pen components. Grey and red
cylindrical parts are snapped onto the green. (Bottom) An injectionmolded ratchet component from a medical device used by Novo
Nordisk for a benchmark study.

16

Next

Contents

Drug Delivery
Teaming up to model polymer behavior
To ensure the integrity of their designs, Hansen and his team in
the Device Simulation department rely on computer simulation
with Abaqus finite element analysis (FEA).
“More than a decade ago, my colleagues and I explored a
number of commercial software codes,” he says. “We chose
Abaqus because it was a well-integrated solution that provided
both implicit and explicit capabilities, and could model the
nonlinear behavior of the fine details in our designs correctly,
including the high number of interfaces in contact.”

FEA analysis demonstrates the ultimate straining of the snap fit on
the device components during assembly.

Over time, the group’s device models have become more refined,
sophisticated, and computationally demanding. The SIMULIA
Polymer Customer Review Team (PCRT) has worked closely
with Novo Nordisk all along to provide updated enhancements
in Abaqus that enable the company to model and predict the
complexities of polymer behavior with increasing accuracy and
efficiency. (The PCRT includes members from the automotive,
high-tech, life sciences, and consumer goods industries.)

Final snap deformation model (left) in Abaqus FEA (using a parallel
framework to capture the changing behavior of the polymer material
under cyclic loading) shows a more accurate plastic deformation of
0.12mm. The earlier elastic-plastic model (right) over-predicted that
deformation would be 0.66mm.

The snap-fit challenge
A recent focus on snap fits in insulin pens demonstrates the
challenges the team has faced when modeling polymers.
“We concentrated on snap fits because they demonstrate
almost ideal cyclic loading, with parts repeatedly loading and
unloading from single to multiple cycles,” says Hansen.

framework makes use of an arbitrary number of viscoelastic
networks and an elastic equilibrium network to create a specific
nonlinear viscoelastic model that is used to predict and track
changes in the internal structural networks of a polymer as the
material responds to repeated cyclic loads during snap fit. The
material parameters in the FEA model are updated at each time
step to reflect the new, altered state of the polymer. Since every
type of polymer shows a different response to temperature,
load, etc., the team continues to explore ways to identify the
material characteristics of different polymer networks.

During such cycles, the viscous nature of the thermoplastic
material determines how the bouncing back to ‘normal’ occurs.
Prediction through analysis of such time-dependent behavior is
key to the device development process.
Since devices can be subjected to different environments,
including elevated temperatures, the function of the device
must be as unaffected as possible by such changes and
always comply with the specification even though the material
properties of the components vary. Even just sitting on a
pharmacy shelf or in the medicine cabinet, polymer materials
are prone to creep and relaxation over time at rates that can
vary with the temperature. Some polymers are also more
complex than others: those used in durable devices may
contain carbon or glass fillers that show anisotrophic behavior,
which can be hard to predict.

Not only are such advanced models useful to designers finetuning the latest pen configuration, the data can help inform
manufacturing processes in the factory. “We have a processsimulating capability, through Moldflow, for which Abaqus has
an interface. This allows us to input the stress fields that result
from the molding process right into our models,” says Hansen.
“As a result, we have greater insight into our manufacturing
process and are more able to design parts that have a very low
level of residual stresses in critical regions.
“SIMULIA is working closely with us to provide capabilities we
need,” says Hansen. “Having material models incorporating
time-dependent viscous behavior is very important for our work
and we’re now able to simulate both creep and relaxation with
Abaqus. We are investigating how well the model will adapt to
different kinds of thermoplastics, which may require different
networks. Calibration will be key going forward.”

“Modeling these diverse material characteristics as well as the
behavior of the polymer as the load induces larger strains closer
to yield is difficult,” says Hansen. “In order to predict such
viscoelasticity precisely, we needed a more refined model that
goes beyond a mainstream elastic-plastic approach.”

Modeling material that’s constantly changing

For More Information

The team is now using the ‘parallel rheological framework’
methodology available in Abaqus to model polymer nonlinear
viscoelasticity with greater accuracy than ever before. The

www.3ds.com/simulia

Previous

www.novonordisk.com
www.3ds.com/SCN-February2013

17

Next

Contents

Orthopedics

Optimization of
Surgical Positioning in
Total Hip Replacement

c)

d)

Jacob Elkins, John Callaghan, Douglas Pedersen, and
Thomas Brown (Department of Orthopaedics and
Rehabilitation, University of Iowa, Iowa City, Iowa)
a)

When joint pain and loss of mobility occur as a result of
end-stage osteoarthritis or other severe hip pathologies, over
250,000 people choose to have total hip replacement (THR)
surgery. Even though THRs are one of the most successful
surgical inventions in medical history, they do fail. THR failures
are often grouped as “early” or “late,” with early failure usually
due to dislocation of the head from the cup, and late failure
frequently due to adverse biologic reaction to wear debris
generated at the bearing surface. Despite nearly six decades of
investigation, the ideal surgical orientation of THR components
remains unclear. Positioning of total hip bearings involves
significant tradeoffs, as cup orientations most favorable in
terms of stability are not necessarily ideal in terms of reduction
of contact stress and wear potential. Previous studies and
models have not addressed these potentially competing
considerations for optimal THA function. Additionally, it is
currently unknown whether the ideal orientation varies on
implant parameters, such as variations in femoral head size.
We, therefore, investigated optimal surgical cup orientation
with a previously generated and physically validated finite
element (FE) model of metal-on-metal THR.

Figure 1. The FE model consisted of bony anatomy (a) and the hip soft tissues (b,
anterior region of capsule rendered transparent for clarity). Four values of femoral
anteversion were considered (c) as were five distinct femoral head sizes (d).

Fitted
ellipse

a)

f)

c)

g)

d)

Isosurface

h)

Figure 2. For every combination of femoral head size and femoral anteversion
(20 such combinations total), the Stability Score (a) and Wear Score (b) are
combined to determine the Performance Score (c). The optimal orientation is
determined as the center of an ellipse fitted to an isosurface of scores > 90 (d).
When considering all 20 combinations, regressions could be performed
demonstrating optimal surgical orientation (e-h).

Discussion
The “landing zone” of ideal cup orientation did not increase with
increased head size, challenging the presumption that larger
heads are more forgiving in terms of stability and durability.
Additionally, ideal cup positioning was considerably more
sensitive to cup anteversion than to inclination. Finally, the
current investigation is the first to quantitatively suggest that
ideal cup positioning varies with both femoral anteversion and
femoral head size.

The FE model consisted of bony anatomy and the hip soft
tissues (see Figure 1). Five dislocation-prone motions as well
as gait were considered, as were permutations of femoral
anteversion (0° to 30°), femoral head diameter (32 mm to
48 mm), cup inclination (25° to 75°), and cup anteversion (0°
to 50°), resulting in 4,320 distinct FE simulations. A novel
metric (“Performance Score”) was developed to delineate
optimized cup orientation by considering both surface wear
and component stability (see Figure 2 A-D).

Positioning THR bearings involves significant tradeoffs with
regard to stability and long-term bearing wear. The computational
analysis identified optimal orientations to balance these
considerations. These tradeoffs help explain the alarming rates
of adverse local tissue response reported for large head metalon-metal THR devices that have demonstrated an improvement
in joint stability. The conclusions from this study can readily be
translated to other hard bearing surfaces—including ceramics
and highly cross-linked polyethylene—suggesting careful
consideration of the choices and compromises in THA design
are required for all bearing couples.

All FE simulations were performed using Abaqus/Explicit.

Results
Ideal cup position was substantially more sensitive to cup
anteversion than to inclination. Regressions demonstrated
strong correlations between optimal cup inclination vs.
head diameter (Pearson’s r = -0.88), between optimal cup
inclination vs. femoral anteversion (r = 0.96), between optimal
cup anteversion vs. head diameter (r = 0.99) and between cup
anteversion and femoral anteversion (r = -0.98) (see Figure 2
E-H).

Previous

b)

e)

Method

www.3ds.com/simulia

b)

For More Information
www.uiortho.com

18

Next

Contents

Orthopedics

Computational Study of Cortical Bone Screw Pullout using
the eXtended Finite Element Method (XFEM)
Emer M. Feerick, Patrick McGarry (National University of Ireland, Galway, Ireland)
Abstract: A study of screw pullout from cortical bone has been conducted using the UDMGINI subroutine with the eXtended finite
element method (XFEM). XFEM alleviates mesh dependency when computationally modeling crack initiation and propagation.
Cortical bone is a naturally occurring composite with a distinctive aligned microstructure that leads to anisotropic material behavior
and damage. In this study, using a UDMGINI subroutine, stress components relative to a predefined osteon orientation were
computed and an anisotropic damage criterion was used to determine damage initiation and to predict crack propagation. A 2D
model of a single screw embedded in cortical bone was generated. A displacement boundary condition was applied to the top
surface of the screw. During pull-out, contact interactions were implemented between the newly formed surfaces during crack
propagation, eliminating unphysical over closure of elements. The model predicted the correct patterns of damage and crack
propagation for both longitudinal (osteons aligned parallel to screw axis) and transverse (osteons aligned perpendicular to screw
axis) pullout tests compared to experimental observations. In both cases crack propagation was predicted in the direction of osteon
alignment.
Keywords: Screw Pullout, Abaqus/Standard, Extended Finite Element Model (XFEM), Crack Propagation, Damage, Failure.

1. Introduction
1.1 The Extended Finite Element Method
The eXtended finite element method (XFEM) is a method developed for computationally modeling crack initiation and propagation
while alleviating mesh dependency. It can be implemented to simulate initiation and propagation of a discrete crack along an
arbitrary, solution dependent path. Additionally contact interactions can be applied to the newly formed surfaces exposed as a
result of crack propagation. Thus the effects of friction between newly exposed surfaces can be accounted for. Previously, work
has been conducted and microstructure models developed to model fracture of cortical bone. These models incorporated a damage
process that initiated damage within an entire element and when the properties of the element degraded such that they could no
longer carry load they were deleted from the mesh nucleating voids within the material that replicated experimentally observed
fracture mechanisms (Feerick et al. 2011). The possibility to alleviate mesh dependency within these materials offers an advance
on previous models as the density of the meshes becomes less influential. Also a material model that is capable of predicting the
correct crack patterns without incorporating a detailed geometry of the microstructure will also facilitate the development of larger
macro scale models requiring significantly lower computational power. A limitation of some current cortical bone failure models
is that they would be too computationally expensive to incorporate in macro scale models of whole bones. Previous studies that
generated models of cortical bone using XFEM incorporated detailed microstructures (Budyn et al. 2010; Abdel-Wahab et al. 2012).
If these previously developed models were applied to a macro scale simulation of whole bones the computational demand would
not be viable. However in the present study we investigate a new feature of Abaqus 6.11 which now facilitates the use of a user
subroutine (UDMGINI) for an anisotropic damage initiation criterion without the need for complex microstructure geometry.

1.2 Cortical Bone
Cortical Bone is a naturally occurring composite, consisting of several constituents that dictate the overall response of the material
during loading. In the present study we simplify the complex microstructure of bone to osteons (the fiber) embedded in an
interstitial matrix. Osteons consist of a series of concentric cylinders of lamellae with a central vascular canal. Osteons are typically
200 µm and are 1-2mm in length and run parallel to the long axis of the bone. The osteons are surrounded by an interstitial matrix
consisting of hydroxyapatite. (Rho et al. 1998; Cowin 2001).

1.3 Screw Pullout
Screws are used for orthopedic applications throughout the human body ranging from fracture fixation plates to spinal fusion
rods. Maximizing the screw insertion depth in cortical bone is regarded as a way of significantly increasing the pullout strength of
the screw (Pollard et al. 2010). A computational model capable of predicting the failure modes and loads that occur during screw
pullout would provide a design tool for the evaluation of future designs as well as design optimization. It is important to understand

www.3ds.com/simulia

Previous

19

Next

Contents

Orthopedics
the failure mechanisms that occur during screw pullout as this may determine the angle of insertion and geometry of screw threads
to maximize the holding power of the screw.

2. Materials & Methods
2.1 UDMGINI Failure Criteria
A new feature of Abaqus 6.11 is the ability to incorporate a user subroutine UDMGINI for user defined damage initiation criterion.
The user may define as many failure indexes as required. In the present study we define two failure indexes for damage initiation
in the x and y directions. Failure index one defines the fiber failure index for initiation. The UTS of the fiber as well as the ultimate
shear failure of the fiber is defined by the user. Once the value of σf =1 the damage is initiated.

Where:

σf 	
σ11 	
σ12 	
σf f 	
σfτ 	
f

is the fiber damage initiation criterion
is the current stress in the local x direction
is the current stress in the local x-y direction
is the UTS of the fiber
is the shear failure strength of the fiber

Failure index 2 defines the matrix failure index for initiation. The UTS of the matrix as well as the ultimate shear failure of the fiber
is defined by the user. Once the value of σm=1 damage is initiated.

Where:

σm 	
σ22	
σ12 	
σ mf 	
σ mτ 	
f

is the fiber damage initiation criterion
is the current stress in the local y direction
is the current stress in the local x-y direction
is the UTS of the matrix
is the shear failure strength of the matrix

As each of the failure indexes are calculated during an iteration of the simulation an array referred to as FNormal is compiled. This
array contains the normal direction to the fracture line for each failure mechanism. This means that as damage is initiated in a
particular failure index, the direction in which the crack will initiate is returned by the subroutine. The crack direction is defined in
terms of the local orientation assigned to the material. For the present study if failure index 1 (fiber failure) is initiated the crack
will propagate in the local y direction, perpendicular to the fiber direction. However, if failure index 2 (matrix failure) is initiated the
crack will propagate in the local x direction, parallel to the fiber direction.

2.2 Single Element Tests
2D single element tests were conducted to examine crack propagation based upon osteon alignment. A unit cell geometry was
created. The material properties assigned to the model are summarized in Table 1. An orientation was applied to the material so that
the local x direction represents the direction of the fibers. A single element was held fixed as shown in Figure 1. The block was held
fixed in the x and y directions at the bottom left corner. Also the top left and right hand corners were held fixed in the x direction. A
1mm displacement boundary condition was applied in the vertical direction. Longitudinal (fibers orientated parallel to the direction
of loading) and transverse (fibers orientated perpendicular to the direction of loading) simulations were conducted.

www.3ds.com/simulia

Previous

20

Next

Contents

Orthopedics
All simulations were conducted using Abaqus/Standard 6.11 with a CPS4 element. For simulation of crack growth a process of
damage evolution (DE) was applied to the model. Energy dissipation was used to determine crack growth. The value selected for
cortical bone energy release was based upon those reported in the literatrure (Abdel-Wahab et al. 2012).
Table 1. Cortical bone material properties.
E1

17100 MPa

σff

233 MPa

E2

10100 MPa

σ fτ f

85 MPa

E12

3300 MPa

σ mf

51 MPa

ν

0.3

σ mτ f

34 MPa

DEm

0.3 N/mm

DEf

0.3 N/mm

Figure 1. Single element boundary conditions.

2.3 Screw Pullout
A 2D model was developed for screw pullout to illustrate the potential application of XFEM for longitudinal and transverse pullout
simulations. A limitation of XFEM is that axisymmetric elements cannot be used with the enrichment feature required for XFEM.
Thus a plane stress model is presented here. The boundary conditions applied to the model are shown in Figure 2. A displacement
boundary condition of 5mm was applied to the top of the screw. The edge of the screw was held fixed in the global x direction.
The edge of the cortical bone was held fixed in both the global x and y directions. Local orientations were applied to assign fiber
directions for each simulation. For a longitudinal simulation the fiber direction (local x) was orientated parallel to the axis of the
screw. For a transverse simulation the fiber direction (local x) was orientated perpendicular to the axis of the screw. 35,000 CPS4
elements were used with element size of 0.05mm.

www.3ds.com/simulia

Previous

21

Next

Contents

Orthopedics

Figure 2. Screw pullout boundary conditions for a longitudinal and transverse simulation.

3. Results
3.1 Single Element Tests
The results for single element simulations are shown in Figure 3. A longitudinal simulation predicts vertical crack formation. A
transverse simulation predicts horizontal crack formation. In both a longitudinal and transverse simulation failure index two reaches
the value of 1 first. Thus cracks will propagate parallel to the direction of the fiber. Hence, cracks grow in the vertical direction for a
longitudinal simulation, while cracks grow in the horizontal direction for a transverse simulation.

Figure 3. Single element results for longitudinal and transverse simulations.

3.2 2D Screw Pullout Simulation
The results for longitudinal and transverse screw pullout simulations are summarized in Figure 4. A longitudinal pullout simulation
predicts localized deformation at the tips of the screw threads. Crack formation is predicted vertically upwards from the screw
thread tips with the material between the screw threads removed with the screw leaving no thread definition on the fracture

www.3ds.com/simulia

Previous

22

Next

Contents

Orthopedics
surface. A transverse pullout simulation predicts horizontal crack formation with deformation extending much further into the
surrounding bone. In both longitudinal and transverse simulations failure index 2 reaches the value 1 first. For both simulations
the cracks propagate parallel to the fiber direction results to significantly different failure modes as the screw is removed from the
cortical bone.

Figure 4. Crack patterns for longitudinal and transverse screw pullout simulations.

4. Discussion
Results highlight the ability of the models to predict significantly different failure mechanisms for cortical bone. This is achieved by
simply altering the orientation assigned to the material. It has previously been shown that fracture patterns during screw pullout
from cortical bone are dependent upon osteon alignment (Feerick et al. 2011). In the past complex geometries of cortical bone with
XFEM have been used to model the alternate fracture patterns of bone depending on osteon alignment (Budyn et al. 2010; AbdelWahab et al. 2012). However, in the present study, we have shown that without generating complex microstructure geometry and
using the UDMGINI sub routine it is possible to define crack directions based upon the fiber (osteon) orientation.

5. Conclusion
Screw pullout from cortical bone is one example of an application for the UDMGINI subroutine with XFEM. However the same
method could be applied to cortical bone simulations for a wide range of orthopedic applications. These applications could range
from orthopedic device design evaluation to modeling fracture incidence in whole bones. The use of the model is not limited to
applications that consider cortical bone but rather any composite material containing a known fiber alignment.

6. References
1.	 Abdel-Wahab, A. A., Maligno, A. R. and Silberschmidt, V. V., 2012. “Micro-scale modelling of bovine cortical bone fracture:
Analysis of crack propagation and microstructure using X-FEM”. Computational Materials Science 52, 128-135.
2.	 Budyn, É. and Hoc, T., 2010. “Analysis of micro fracture in human Haversian cortical bone under transverse tension using
extended physical imaging”. International Journal for Numerical Methods in Engineering 82, 940-965.
3.	 Cowin, S. C., 2001. Bone mechanics handbook, CRC press USA.
4.	 Feerick, E.M.; Mullett, H; FitzPatrick, D; McGarry, J.P. “Computational Investigation of Cortical Bone Failure Mechanisms During
Screw Pullout”, 19th Annual Symposium on Computational Methods in Orthopaedic Biomechanics (Pre-ORS), Long Beach,
California, January 12, 2011.

www.3ds.com/simulia

Previous

23

Next

Contents

Orthopedics
5.	 Pollard, J. D., Deyhim, A., Rigby, R. B., Dau, N., King, C., Fallat, L. M. and Bir, C., 2010. “Comparison of Pullout Strength between
3.5-mm Fully Threaded, Bicortical Screws and 4.0-mm Partially Threaded, Cancellous Screws in the Fixation of Medial Malleolar
Fractures”. The Journal of foot and ankle surgery : official publication of the American College of Foot and Ankle Surgeons 49,
248-252.
6.	 Rho, J.-Y., Kuhn-Spearing, L. and Zioupos, P., 1998. “Mechanical properties and the hierarchical structure of bone”. Medical
Engineering & Physics 20, 92-102.

7. Acknowledgements
The authors would like to acknowledge the assistance of Cheryl Liu at Dassault Systemes Simulia Corp for her assistance in the
development of the failure criteria for the UDMGINI subroutine for the applications presented in this study.

www.3ds.com/simulia

Previous

24

Next

Contents

Cardiovascular

SIMULIA Spearheads The Living Heart Project
Bringing the medical community together for improved patient care
The unveiling of Dassault Systèmes’ Living Heart Project at
this year’s SIMULIA Community Conference (SCC) is certainly
fitting. It represents the culmination of many months of
collaboration, planning, and software development that were
launched into action after the same event in 2013.
While reviewing presentations in the medical special-interest
group at the SCC in Vienna last year, Steven Levine, Senior
Director, SIMULIA Portfolio Management, observed that
simulation technology was approaching a tipping point in
replicating the physics of living systems. Aware of the impact
of cardiovascular disease, currently the #1 cause of death
worldwide, he and his team had already performed introductory
studies to explore the possibility of simulating a human heart.
However, to make the translation from proof-of-concept to
clinical practice, it would require an effort beyond what
SIMULIA could provide alone.

watching the project with great interest, as part of a Regulatory
Science initiative. Its aim is to inspire leadership from the
community—and SIMULIA and collaborators have stepped up
through The Living Heart Project.”

The plan was simple, yet visionary: use SIMULIA’s technology
and The 3DEXPERIENCE® platform to challenge the scientific
community to create a realistic 3D simulation of a beating,
human heart—one that could be employed as a foundation
to build personalized models based on a patient’s measured
data to gain insight into potential cardiovascular diseases and
methods of treatment. Levine approached Bernard Charlès,
Dassault Systèmes CEO, with the idea—and The Living Heart
Project was born.

The project has high aspirations. Simulating an entire human
heart that beats realistically is a highly coupled, multiscale, multi-physics problem layered on top of a complex
materials engineering problem. Initially, the focus will be
electromechanical applications for device design, such as
pacemaker leads, stents, and artificial valves. However, project
collaborators are already exploring its applicability to study
treatment for heart disease or personalized ventricular assist
devices. Future versions will include blood flow and thrombosisrelated applications, with increasingly detailed models of
electrical pathways eventually reaching the cellular level. The
knowledge developed under this project could one day pave the
way for fully 3D-printed bioficial hearts.

While simulation has gained broad adoption in many industries,
healthcare remains the last frontier. Strides have certainly
been made in medical device development, but the extreme
complexity of the human body continues to hamper the quest
for lifelike accuracy. What’s more, without such virtual models,
regulators often have no choice but to accept laboratory
testing as validation for new medical devices, with failures and
recalls grabbing news headlines with disappointing regularity.
Without question, simulation has great potential to reduce
testing risks and costs while further benefiting patient safety
and health.

“This project began with a challenge to realistically simulate
the physics of the heart,” says Levine. “As it evolves, we
already see it beginning to transform how people think about
what is possible and how this could serve as a platform to
translate science into meaningful medical practice. We can
imagine that this will lead to a new paradigm for data delivered
directly into the hands of physicians; not simply 2D patient
scan data, but fully analyzed 3D models where abnormal
behavior is clearly identified and treatment options evaluated.
The 3DEXPERIENCE platform could connect an ecosystem of
providers to allow any physician access to these state-of-theart diagnostics.”

In launching the project, SIMULIA reached out to researchers,
device developers, and cardiac physicians, emphasizing the
validity—and potential—of the effort. It was quickly recognized
that, if successful, they would have an important platform for
medical innovation.
“There are millions of patients out there who really need this
technology today, including my own daughter,” Levine says.
“Participation has gained momentum beyond our expectations.
Before we release it, the first model will be tested by the
community in their respective disciplines. In fact, the FDA is

www.3ds.com/simulia

Previous

For More Information
www.3ds.com/heart

25

Next

Contents

Cardiovascular

Abdominal Aorta Blood Flow
Analysis with Abaqus/CFD
Dassault Systèmes SIMULIA would like to thank Simpleware
for providing the geometric model of the abdominal aorta.
The progression of disease in arteries is affected by local
blood flow characteristics. Arterial structural features, such
as branches, bifurcations, and irregular shape and curvature
changes, introduce complexities to the blood flow field. Further
insight on how such mechanical factors affect arterial health
can be gained by employing computational fluid dynamics
tools. With this approach, accurate and detailed quantitative
data on hemodynamic factors can be obtained.
In this Technology Brief, we demonstrate the use of Abaqus/
CFD to simulate complex blood flow in an abdominal aorta and
its branches. Blood velocity profiles and wall shear stresses
are computed, and flow recirculation is visualized. It will be
shown that blood flow analyses using Abaqus/CFD can provide
comprehensive data that would be difficult to obtain from
experimentation.

Figure 1: Abdominal aorta

incompressible viscous flow solver based on an advanced
second-order projection method that uses a node-centered
finite-element discretization for the pressure. This hybrid
approach guarantees accurate solutions and preserves the
local conservation properties associated with traditional finite
volume methods. Abaqus/CFD can be used with unstructured
grids and with various element types such as hexahedral,
tetrahedral, wedge, and pyramid. The Abaqus/CFD solution
methodology incorporates iterative Krylov solvers with algebraic
multi-grid (AMG) preconditioning for solving the pressurePoisson equation.

Background
Many research studies have highlighted the importance of
hemodynamic factors in vascular disease progression. The
prevalence of disease in specific arterial locations is often
correlated with flow features such as wall shear stress,
stagnation, and recirculation zones.
For example, atherosclerotic lesions are more likely to be found
along the outer wall of the carotid sinus region of the carotid
artery where the wall shear stress is low, while they are less
likely to occur along the inner wall of carotid sinus where the
wall shear stress is high [1]. Similarly, atherosclerotic disease
is more likely to develop in the abdominal aorta below the
diaphragm [2]. Dilatation of the aortic wall in the abdomen as
seen in aortic aneurysms is also affected by blood flow behavior;
the growth and rupture of such aneurysms are affected by the
hemodynamics [3]. Renal artery stenosis may also alter flow
characteristics; narrowing of the artery can cause the blood
flow regime to transition from primarily laminar to turbulent.
Changes in the shear stress distribution affect the blood
pressure as well as the evolution of atherosclerosis [4]. Surgical
procedures such as coronary artery bypass, where grafts
reroute blood flow across a blocked artery, can cause changes
in flow characteristics that initiate a biological response such as
local growth and remodeling of surrounding tissues.

Key Abaqus/CFD Features and Benefits
•	 Transient incompressible viscous fluid flow analysis
•	 Non-Newtonian viscosity models
•	 Luminal time dependent pressure waveform boundary
conditions
•	 Pulsatile flow velocity boundary condition, derived from
Womersley theory, implemented in a velocity user subroutine
•	 Third-party boundary layer mesh import in Abaqus/CAE

Geometry and Model
Patient-specific abdominal aorta geometry is obtained using the
ScanIP software from Simpleware. ScanIP reads CT scan files
and creates a high quality triangulation for STL export. The STL
file of the abdominal aorta is imported into CATIA V6R2013. A
surface reconstruction is performed to obtain the CAD model
of the abdominal aorta. The geometry is then trimmed at the
boundaries to make it suitable for CFD analysis. Figure 1 shows
the final CAD model of the abdominal aorta and branches.

In this Technology Brief, we show that Abaqus/CFD can be
used to model the pulsatile blood flow in the abdominal aorta
and its various branches that supply blood to the organs in
the abdomen. Abaqus/CFD uses a time accurate transient

www.3ds.com/simulia

Previous

26

Next

Contents

Cardiovascular
Mesh

Table 1: Non-Newtonian fluid properties for the Carreau-Yasuda model

The CAD model is imported into ANSA v13.1.4 and meshed
with a boundary layer meshing technique. A total of 6 wedge
layers are generated along the wall. The mesh consists of
192375 nodes and 568769 elements. Of these, 314873 are
tetrahedral elements, 249660 are wedge elements and 4236
are hexahedral elements. The mesh is exported in the Abaqus
input file format and then imported into Abaqus/CAE where
the CFD analysis attributes are defined. A close-up view of the
boundary layer mesh is shown in Figure 2.

Shear viscosity at low shear rate

0.025

Shear viscosity at high shear rates

0.0035

Time constant

25

Flow behavior index

.025

Material constant

2.0

Figure 3: Viscosity v. shear rate for non-Newtonian blood properties
Figure 2: Boundary layer mesh of the abdominal aorta

the major branches of the abdominal aorta: the superior
mesenteric artery, the left and right renal arteries, the inferior
mesenteric artery, and the left and right iliac arteries.

Material
Blood is modeled as an incompressible fluid with a density
r = 1050 kg/m3 and a non-Newtonian viscosity described by
the Carreau-Yasuda model. The model properties are listed in
Table 1. The Carreau-Yasuda model describes the shear thinning
behavior of blood. It is often a reasonable approximation
to treat blood as a Newtonian fluid in blood vessels greater
than 0.5 mm in diameter. In vessels of such large diameter,
the viscosity is relatively constant due to high rates of shear.
However, non-Newtonian effects need to be accounted for in
smaller vessels. The dependence of viscosity on the shear rate,
on a log-log scale, is shown in Figure 3.

A no-slip/no-penetration wall boundary condition is applied
on the wall surface of the abdominal aorta. A time-dependent
luminal pressure waveform as shown in Figure 4 is applied at
the left and right iliac artery outlets. The pressure waveform is a
triphasic pulse appropriate for normal hemodynamic conditions
in the abdominal aorta below the renal artery ([7], [8]). The
volumetric flow rates at the abdominal aorta inlet and superior
mesenteric, left and right renal, and inferior mesenteric artery
outlets are shown in Figure 5 ([5], [6]).
These flow rate waveforms are represented as complex Fourier
series and their coefficients are utilized for evaluating the
velocity boundary condition at the inlet and outlets (except for
the left and right iliac artery):

Analysis Procedure
A transient incompressible laminar fluid flow analysis is
performed in Abaqus/CFD. Abaqus/CFD uses a projection
method that enables segregation of pressure and velocity
fields for efficient solution of the incompressible Navier-Stokes
equations. The advective and diffusive fluxes are both treated
implicitly. The solution method uses a second-order accurate
least-squares gradient estimation. Simulation is performed for
4 cardiac cycles, with each cycle lasting 0.9 sec.

Boundary Conditions
Solving the transient Navier-Stokes and continuity equations
requires specification of appropriate initial and boundary
conditions. Blood enters the abdominal aorta at a level below
the celiac artery. The model has six outlet sections representing

www.3ds.com/simulia

Previous

27

Next

Contents

Cardiovascular
Womersley theory [9] evaluates the velocity profile at any crosssection for the unsteady, laminar flow of an incompressible fluid
through a pipe of constant radius, when a time-varying pressure
gradient is applied. The profile is defined as

Here Qk represents the complex Fourier coefficients of
the flow waveforms. The velocity boundary condition is
specified at the inlet and outlets through user subroutine
SMACfdUserVelocityBC. The user subroutine calculates a
spatially- and time-varying velocity boundary condition based
on Womersley theory and the volume flow rates shown in
Figure 5.

In this equation, αn is the non-dimensional Womersley number,
where R denotes the radius, ρ the density, μ the viscosity and ω
the circular frequency. Real(*) denotes the real part of a complex
number and J0 and J1 are Bessel functions of the first kind of
order 0 and 1, respectively.
The choice of the Womersley velocity profile is a significant
improvement over spatially-constant or parabolic velocity
profiles since it captures the flow reversal at artery outlets
and provides a more realistic boundary condition for pulsatile
flows. Womersley theory, however, is limited to circular crosssections and hence, a mapping method ([10]) is used to map the
evaluated velocity profile to non-circular cross-sections.

Figure 4: Luminal pulsatile pressure waveform

Velocity User Subroutine Implementation
User subroutine SMACfdUserVelocityBC is written in C and
provides access to the necessary model and solution quantities.
Access is also provided to the MPI communicator and MPI
routines for parallel implementation. The following steps are
performed in the user subroutine to calculate velocities on the
inlet and outlets:
1.	 Find the facets of elements lying on the edge of the surface.
This is accomplished using a convex hull algorithm for
finding edge points amongst a point cloud.
2.	 Transform all surface points to a local coordinate system
specified by the user and with its origin lying at the surface
centroid.
3.	 Express the points in a polar coordinate system with its origin
located at the surface centroid.
4.	 Evaluate a normalized radius [10]:
a. For each internal point (r, θ), find the two edge nodes i and
j with θ coordinates such that θi ≤ θ ≤ θj.

Figure 5: Flow rates for abdominal aorta, superior mesenteric
artery, left and right renal artery, and inferior mesenteric artery

www.3ds.com/simulia

Previous

b. Evaluate a normalized radius R by linearly interpolating
the radius of edge nodes, Ri and Rj. The normalized radius
serves as the radius for all points lying on the line connecting
an edge point to the centroid.

28

Next

Contents

Cardiovascular
5.	 Evaluate the Womersley velocity using the radius r,
normalized radius R, current time value, and Fourier
coefficients for the flow waveform.
6.	 Transform the velocity to the global system for boundary
condition specification.
The above procedure can be skipped for circular surfaces. For
these cases, the radius of the circular surface serves as the
normalized radius for all internal points.

Results and Discussion
The kinetic energy of the fluid domain is plotted in Figure
6. It can be seen that the analysis reaches a steady periodic
condition after 1 cardiac cycle.
Figure 7 shows the velocity vectors along the midsagittal plane
in the infrarenal segment of the abdominal aorta at various
times during the cardiac cycle. A mildly recirculating flow
pattern can be seen near the anterior wall; this vortex is formed
during late systole and remains until peak systole when the
effect of strong systolic acceleration results in a completely
attached flow along the anterior wall. The particular location
of the vortex is due to the curvature of the infrarenal segment
which forces the flow to expand suddenly about the corner. The
anterior portion of the abdominal aorta displays lower velocities
than the core and posterior portion.

Figure 6: Total kinetic energy of the flow

Wall shear stresses during the cardiac cycle are contoured in
Figure 8. The shear stress magnitude can be used to estimate
the possibility of rupture of aneurysms.
Figure 9 shows the surface traction vectors on the infrarenal
segment of the abdominal aortic walls. Relatively low tractions
are found compared to the branches where the velocities are
higher. The primary movement of surface traction vectors
is from the top of the vessel to the bottom, following the
movement of the vortex during the cardiac cycle.

Figure 7: Velocity vector field (m/sec) along the midsagittal plane in the
infrarenal segment of abdominal aorta at (a) the beginning of the cycle,
(b) peak systole, (c) end systole, and (d) mid diastole. Midsagittal plane
definition (far right) and schematic representation of various time points
during cardiac pressure cycle

In Figure 10, streamlines are shown for the end systole. The
streamlines show the formation of the vortex in the infrarenal
segment and significant mixing of the blood. Figure 11 plots
the percentage error in the mass balance, defined as the ratio
of the sum of the inlet and outlet volume flow rates to the inlet
volume flow rate.

Summary
In this technology brief, we demonstrate that Abaqus/CFD can
be used to model pulsatile blood flow in the abdominal aorta
and the various branches that supply blood to the organs in
the abdomen. Newtonian and non-Newtonian blood properties
can be modeled and specialized boundary conditions can
be implemented with a user subroutine. This study can be
easily extended to perform CFD analyses of diseased arteries.
Parametric and quantitative studies on various hemodynamic

www.3ds.com/simulia

Previous

Figure 8: Wall shear stress (N/m2) at (a) the beginning of the cycle, (b)
peak systole, (c) end systole, and (d) mid diastole

29

Next

Contents

Cardiovascular
3.	 David A. Vorp, “Biomechanics of abdominal aortic
aneurysm,” Journal of Biomechanics, Vol. 40, pp. 18871902, 2007.
4.	 F. Liang, Ryuhei Yamaguchi, H. Liu, “Fluid dynamics in
normal and stenosed human renal arteries: an experimental
and computational study,” Journal of Biomechanical Science
and Engineering, Vol. 1, pp. 171-182, 2006.
5.	 D. N. Ku, S. Glagov, J. E. Moore, and C. K. Zarins, “Flow
patterns in the abdominal aorta under simulated
postprandial and exercise conditions: an experimental
study,” Journal of Vascular Surgery, Vol. 9, pp. 309–316,
1989.
6.	 Charles A. Taylor, Thomas J.R. Hughes, and Christopher K.
Zarins, “Finite Element Modeling of Three-Dimensional
Pulsatile Flow in the Abdominal Aorta: Relevance to
Atherosclerosis,” Annals of Biomedical Engineering, Vol. 26,
pp. 975-987, 1998.

Figure 9: Surface traction (N/m2) vector plot in the infrarenal segment at
(a) the beginning of the cycle, (b) peak systole, (c) end systole, and (d)
mid diastole

7.	 C. J. Mills, I. T. Gabe, J. H. Gault, D. T. Mason, J. Ross Jr., E.
Braunwald, and J. P. Shillingford, “Pressure-flow
relationships and vascular impedance in man,”
Cardiovascular Research, Vol. 4, pp. 405-417, 1970.
8.	 C. M. Scotti, E. A. Finol, “Complaint biomechanics of
abdominal aortic aneurysms: A fluid-structure interaction
study,” Computers and Structures, Vol. 85, pp. 1097-1113,
2007.
9.	 J. R. Womersley, “Method for the calculation of velocity, rate
of flow and viscous drag in arteries when the pressure
gradient is known,” Journal of Physiology, Vol. 127, pp.
553-563, 1955.
10.	 J. P. Mynard, P. Nithiarasu, “A 1D arterial blood flow model
incorporating ventricular pressure, aortic valve and regional
coronary flow using the locally conservative Galerkin (LCG)
method,” Communications in Numerical Methods in
Engineering, Vol. 24, pp. 367-417, 2008.
11.	 Abaqus 6.13 Analysis User’s Manual, Dassault Systèmes
Simulia Corp, 2013.

Figures 10 & 11: End systole streamlines (left), mass balance
percentage error (right)

factors could be performed as well. Such studies can offer
insight into the temporal and spatial variations of velocity and
pressure fields as well as variations of wall shear stresses in
vascular geometries.

References
1.	 C. K. Zarins, D. P. Giddens, B. K. Bharadvaj, V. S. Sottiurai, R.
F. Mabon, and S. Glagov, “Carotid bifurcation
atherosclerosis: Quantitative correlation of plaque
localization with flow velocity profiles and wall shear
stress,” Circulation Research, Vol. 53, pp. 502–514, 1983.
2.	 J. C. Roberts, C. Moses, and R. H. Wilkins, “Autopsy studies
in atherosclerosis: I. Distribution and severity of
atherosclerosis in patients dying without morphologic
evidence of atherosclerotic catastrophe,” Circulation, Vol.
20, pp.511–519, 1959.

www.3ds.com/simulia

Previous

30

Next

Contents

Cardiovascular

Fluid-Structure Interaction
Analysis of a Prosthetic
Aortic Valve using
Abaqus/Explicit Smoothed
Particle Hydrodynamics
Dassault Systèmes SIMULIA would like to thank Dr. Nandini
Duraiswamy of the U.S. Department of Health and Human
Services FDA Center for Devices and Radiological Health.
Durability is a key measurement of prosthetic heart valve
function. Assessment of fatigue life requires accurate estimates
of the stresses induced during the cardiac cycle. Finite element
(FE) studies have been used to estimate peak stresses in valves
[1], and computational fluid dynamics (CFD) studies have been
used to model blood flow around valves [2]. Fluid-structure
interaction (FSI) studies are less common, in part because the
closure of the valve creates CFD domain pinching.

CFD can model the behavior of the blood, and a coupled fluidstructure interaction (FSI) analysis can capture the effect of the
blood on the valve during the cardiac cycle.
There is a final condition during the cycle that presents a
challenge to coupled FSI: the fluid domain pinches during
valve closure, which is a condition most CFD packages cannot
handle. The Smoothed Particle Hydrodynamics (SPH) analysis
technique, available in Abaqus/Explicit 6.11-1, addresses this
challenge and makes modeling heart valves for the entire cardiac
cycle possible, thus increasing the accuracy of prosthetic valve
stress analysis.

The smoothed particle hydrodynamic (SPH) analysis method
in Abaqus/Explicit overcomes this difficulty. In this Technology
Brief, the SPH technique will be used to determine the FSI
response of a generic prosthetic heart valve.

Background

Key Abaqus Features and Benefits

There are two principal modes of aortic valve disease: aortic
stenosis, in which the valve no longer fully opens, and aortic
regurgitation, in which the valve no longer fully closes. Either
condition can eventually require the implantation of a prosthetic
valve to replace the underperforming original.

•	 Abaqus/Explicit Smoothed Particle Hydrodynamics capability
for analyses involving extreme deformation
•	 Robust hyperelastic material modeling
•	 General contact capability for simplified definition of contact
interactions

Surgically implanted or transcatheter-delivered bioprosthetic
aortic valve leaflets undergo dynamic cyclic loading and
large deformation during the cardiac cycle. This can cause
fatigue failure of the leaflets, compromising valve function and
potentially affecting the patient. Accurate stress analysis of
the valve during operation is therefore essential for designing
durable aortic valves and improving patient outcomes.

Analysis Approach
Smoothed Particle Hydrodynamics
SPH offers several advantages over CFD and coupled EulerianLagrangian methods in tracking free surface boundaries,
handling small material-to-void ratios, and modeling extreme
deformation with fragmentation. The latter capability makes it
ideal for simulating the behavior of blood during valve closure
and pressure changes.

The operating conditions of the aortic valve are complex. The
pressure on the aorta side of the leaflets is lower than that on
the ventricular side when the ventricle is pumping oxygenated
blood into the aorta, and the pressure on both sides varies
depending on the stage of the cardiac cycle. This can be
modeled by applying dynamic pressure loads (corresponding
to loads measured in the aorta and left ventricle) directly onto
the leaflets, which is an improvement in accuracy compared to
previous analyses that used only static load conditions.

SPH is part of a larger family of meshless numerical methods that
define a body by a collection of points, instead of using nodes
and elements. The SPH method implemented in Abaqus 6.111 uses a cubic spline kernel for interpolation, applying either a
fixed or a variable “smoothing” length to particles. Internally,
particle connectivity is determined based on smoothing length.
The particles can contact Lagrangian bodies (in this case, the

Even this method, however, does not account for the inertial
and viscous effects of blood contacting the leaflets during flow.

www.3ds.com/simulia

Previous

31

Next

Contents

Cardiovascular
valve leaflets) through the Abaqus/Explicit general contact
feature. In addition, particles can be “glued” to Lagrangian
bodies through *TIE constraints. SPH supports an extensive
library of solid and fluid materials, including user materials.
For this particular simulation, a finite volume of blood near
the aortic valve was modeled with one-node PC3D elements.
All particles had the same volume initially. There were 4956
particles, each with a radius of 1 mm.

Material Modeling
A generic aortic valve was meshed with shell (S4) elements.
The valve had a diameter of 26mm and a thickness of 0.5mm.
The junction between the aorta and the left ventricle was
represented with a rigid tube, and two rigid plates were used
to apply pressure on either side of the fluid particles (Figure 1).
Figure 2: Leaflet material test data and Marlow model representation

Figure 1: Aortic valve model with SPH particles

The material for the valve was modeled with the Marlow
isotropic hyperelastic representation, the general first-invariant
hyperelastic material model in Abaqus. This model can exactly
duplicate physical test data from one of several standard modes
of loading (uniaxial, biaxial, or planar). It works well in situations
where extensive data for one of the test modes is available. For
the present analysis, uniaxial tensile test was used. (Figure 2).
Figure 3: Pressure load profiles [1]

Boundary and Loading Conditions

Results and Conclusions

Translational degrees of freedom were fixed for the valve edges.
Left ventricle and aorta pressure profiles were applied to the
end plates (Figure 3) [1]. The pressure profiles start from the
point at which the pressure inside the left ventricle and the
aorta are the same since the initial condition of the valve was
stress-free. The same pressure was applied to the fluid as an
initial condition. The end plates were not allowed to rotate, and
because the finite volume of fluid is incompressible, the two
rigid plates were constrained to have the same displacement
along the axial direction using an equation constraint.

Peak stress in the valve leaflets occurs during the diastolic phase,
when the valve leaflets are closed. Higher stresses are observed
in the FSI analysis using SPH than the reference model (Figure
4). In addition, the distribution of stresses is also different.
Stress hot spots are observed in the middle of the leaflets as well
as near the corners where two leaflets meet. This shows that, in
addition to the pressure loads, the inertia effect of the fluid also
influences the stress analysis results.
The present SPH simulation capability is an important step
toward providing prosthetic valve designers with increased
simulation accuracy and the data needed to design more durable
valves.

As a reference model, a second analysis was run with the same
(uniform) pressure profiles directly applied on the valve leaflets
without the fluid. All other conditions were the same as the FSI
model.

www.3ds.com/simulia

Previous

32

Next

Contents

Cardiovascular

Figure 4: Mises stress on the leaflets during diastolic phase, SPH model (left) and reference model (right)

References
1.	 H. Kim, J. Lu, M. S. Sacks and K. B. Chandran, “Dynamic
Simulation of Bioprosthetic Heart Valves Using a Stress
Resultant Shell Model”, Annals of Biomedical Engineering,
Vol. 36, No. 2, 2008
2.	 W. Sun, K. Li, and E. Sirois, ”Simulated Elliptical
Bioprosthetic Valve Deformation: Implications for
Asymmetric Transcatheter Valve Deployment”, Journal of
Biomechanics 43(16), 2010,
3.	 Abaqus References
4.	 For additional information on the Abaqus capabilities
referred to in this brief, please see the following Abaqus
6.13 documentation reference:
5.	 Analysis User’s Guide:
hydrodynamic analyses

www.3ds.com/simulia

15.2.1

Smoothed

particle

Previous

33

Contents

Dassault Systèmes, the 3DEXPERIENCE Company, provides business and people with virtual universes to imagine sustainable innovations. Its world-leading
solutions transform the way products are designed, produced, and supported. Dassault Systèmes’ collaborative solutions foster social innovation, expanding
possibilities for the virtual world to improve the real world. The group brings value to over 170,000 customers of all sizes in all industries in more than 140
countries. For more information, visit www.3ds.com.

Americas
Dassault Systèmes
175 Wyman Street
Waltham, Massachusetts
02451-1223
USA

Europe/Middle East/Africa
Dassault Systèmes
10, rue Marcel Dassault
CS 40501
78946 Vélizy-Villacoublay Cedex
France

Asia-Pacific
Dassault Systèmes K.K.
ThinkPark Tower
2-1-1 Osaki, Shinagawa-ku,
Tokyo 141-6020
Japan

©2014 Dassault Systèmes. All rights reserved. 3DEXPERIENCE, the Compass icon and the 3DS logo, CATIA, SOLIDWORKS, ENOVIA, DELMIA, SIMULIA, GEOVIA, EXALEAD, 3D VIA, BIOVIA, NETVIBES, and 3DXCITE are commercial trademarks
or registered trademarks of Dassault Systèmes or its subsidiaries in the U.S. and/or other countries. All other trademarks are owned by their respective owners. Use of any Dassault Systèmes or its subsidiaries trademarks is subject to their express written approval.

Our 3DEXPERIENCE Platform powers our brand applications, serving 12 industries, and provides a rich
portfolio of industry solution experiences.



Source Exif Data:
File Type                       : PDF
File Type Extension             : pdf
MIME Type                       : application/pdf
PDF Version                     : 1.7
Linearized                      : No
Create Date                     : 2014:07:02 16:00:53-04:00
Creator                         : Adobe InDesign CS6 (Windows)
Modify Date                     : 2014:07:14 15:11:07-04:00
Has XFA                         : No
XMP Toolkit                     : Adobe XMP Core 5.2-c001 63.139439, 2010/09/27-13:37:26
Metadata Date                   : 2014:07:14 15:11:07-04:00
Creator Tool                    : Adobe InDesign CS6 (Windows)
Instance ID                     : uuid:1ee0c640-d6d0-4942-a0b6-f0ce3a6d04ef
Original Document ID            : xmp.did:40F58A5B04DAE311A105EF713449498F
Document ID                     : xmp.id:40E42C262302E411B948A0887374B344
Rendition Class                 : proof:pdf
History Action                  : converted
History Parameters              : from application/x-indesign to application/pdf
History Software Agent          : Adobe InDesign CS6 (Windows)
History Changed                 : /
History When                    : 2014:07:02 16:00:53-04:00
Derived From Instance ID        : xmp.iid:3FE42C262302E411B948A0887374B344
Derived From Document ID        : xmp.did:40F58A5B04DAE311A105EF713449498F
Derived From Original Document ID: xmp.did:40F58A5B04DAE311A105EF713449498F
Derived From Rendition Class    : default
Format                          : application/pdf
Producer                        : Adobe PDF Library 10.0.1
Trapped                         : False
Page Layout                     : SinglePage
Page Count                      : 34
EXIF Metadata provided by EXIF.tools

Navigation menu