Chapter 3

Simulation-Based Engineering

Melih Cakmakci, Gullu Kiziltas Sendur and Umut Durak

Abstract Engineers, mathematicians, and scientists were always interested in

numerical solutions of real-world problems. The ultimate objective within nearly all
engineering projects is to reach a functional design without violating any of the
performance, cost, time, and safety constraints while optimizing the design with
respect to one of these metrics. A good mathematical model is at the heart of each
powerful engineering simulation being a key component in the design process. In
this chapter, we review role of simulation in the engineering process, the historical
developments of different approaches, in particular simulation of machinery and
continuum problems which refers basically to the numerical solution of a set of
differential equations with different initial/boundary conditions. Then, an overview
of well-known methods to conduct continuum based simulations within solid
mechanics, fluid mechanics and electromagnetic is given. These methods include
FEM, FDM, FVM, BEM, and meshless methods. Also, a summary of multi-scale
and multi-physics-based approaches are given with various examples. With con-
stantly increasing demands of the modern age challenging the engineering devel-
opment process, the future of simulations in the field hold great promise possibly
with the inclusion of topics from other emerging fields. As technology matures and
the quest for multi-functional systems with much higher performance increases, the
complexity of problems that demand numerical methods also increases. As a result,
large-scale effective computing continues to evolve allowing for efficient and
practical performance evaluation and novel designs, hence the enhancement of our
thorough understanding of the physics within highly complex systems.

M. Cakmakci (&)
Department of Mechanical Engineering, Bilkent University, Ankara, Turkey
e-mail: melihc@[Link]
G. Kiziltas Sendur
Mechatronics Engineering Program, Faculty of Engineering and Natural Sciences,
Sabanci University, Istanbul, Turkey
e-mail: gkiziltas@[Link]
U. Durak
Institute of Flight Systems, German Aerospace Center (DLR), Braunschweig, Germany
e-mail: [Link]@[Link]

© Springer International Publishing AG 2017 39

S. Mittal et al. (eds.), Guide to Simulation-Based Disciplines,
Simulation Foundations, Methods and Applications,
DOI 10.1007/978-3-319-61264-5_3
40 M. Cakmakci et al.

Keywords Engineering design cycle V-process Waterfall model Hardware-




in-the-loop simulations Feature-in-the-loop simulations Component-in-the-loop



simulations Continuum mechanics Computational electromagnetics Partial




differential equations (PDE) Finite element method (FEM) Finite-difference



method (FDM) Multi-scale methods Lumped parameter models Model-based




control system design Vehicle dynamics models Networked control systems





Discretized systems Quantization Observer models Iterative learning

3.1 Introduction

3.1.1 Overview of the Engineering Design Process

The ultimate objective of all engineering projects is to reach a functional design

without violating any of the performance, cost, time, and safety constraints often
optimizing the design for one of them. Generally, in the beginning of each project
high-level requirements for the system is developed. These high-level requirements
can be as literal as “The fuel consumption of the vehicle shall be 40 mpg or more.”
or comparative such as “The new CNC machine will be as precise as our com-
petitors.” Then, these high-level requirements are cascaded down to the lower
levels of the system design steps to obtain the well-defined engineering design
problems.
Engineering design problems are concrete problem constructs that contain
quantifiable performance and constraint metrics. The inputs to the engineering
problems are the performance constraints, design parameters, external conditions.
The output of the engineering design process is the design communicated in tech-
nical terms such as materials, dimensions, and algorithms. Usually, a lesser focused
output of an engineering project is the operation recommendations, lifecycle
maintenance, and storage instructions. In general, main steps of the engineering
design process can be given as, requirements analysis, design, implementation,
verification, and maintenance.
In Fig. 3.1, inputs and outputs of the engineering development process is given
as discussed previously. It is also important to note that this process can be applied
at the component, the sub-system (i.e., group of interrelated components), and at the
system level.
Over time, two primary approaches emerged to approach the solution of com-
plex engineering design projects.
The early approach also known as the “Waterfall Design Process”, the sub-
problems can be tackled and solved sequentially. Even though, it provides a
structured method to perform design and testing tasks, its sequential nature fails to
catch design-related errors early in the development process.
Inspired from the approaches in development of software intensive systems, as
an extension of “Waterfall Design Process”, a new engineering design approach has
emerged which is called the “V-process” where relations of design and validations
3 Simulation-Based Engineering 41

Fig. 3.1 Engineering development process

System System

Sub-System Sub-System

Component Component

DESIGN VALIDATION

Fig. 3.2 Engineering V-process

steps are stressed. Simulations of varying resolution and fidelity become important
tool in the V-process, in order to conduct validations as well as the evaluations of
design decisions before the actual prototype of the system can be build.
In Fig. 3.2, the typical steps of the V-process is given based on (Ulsoy et al.
2012). The essence of the V-process is to cascade from the system level to the
smaller scale such as the component level and the level-based validation of the
work to catch problems at early stages. The different levels of validation and design
work in the V-process increase the importance of effective simulations throughout
the whole process.
Today almost all of the engineering community is using the iteration based
V-diagram process. One of the early hesitation points regarding the engineering
V-process was also its strongest feature, namely the existence of stepwise iterations
and the cost they bring to the overall development. However, the evolution of
42 M. Cakmakci et al.

advanced simulation techniques made iterations more manageable minimizing the

cost of rework during development.
One of the most important topics in the simulation development process is the
decision of the feature content and their fidelity. Too much content or dynamics and
the simulation will be consuming too much computational resources generally
resulting in time and cost problems. Too little detail in simulations will result in
misguided simulations, missing important modes of the target systems and taking
away the benefit of simulation-based iteration in the development process.
When the engineering V-process is considered, the level of validation increases
as the project progresses in time as shown in Fig. 3.2. The decomposition of the
requirements and development process requires a proof-of-concept simulation first,
which are detailed in physics but abstract at the interface level. When the test phase
starts these component-based simulations are combined to produce sub-system and
system-level simulations that are also advanced in terms of the mechanical
and electronics interactions (the interface) of the system. The sub-system and sys-
tem-level simulations are more detailed than the component-level simulations.
In most engineering development projects, the understanding of the target sys-
tem improves with the progress of the project. Therefore, the resolution and fidelity
of the simulations can also be improved using the new data and understanding of
system of interest.
Generally, in the early stages of the engineering development process a proto-
type of the target system does not exist. However, a concept emulating simulation
of the system can be developed using existing models from the company’s
resources or from the existing technical literature. When these simulations are
functional, the new feature of the system can be included in the model and the
simulations can be used to make early predictions about the performance of the
target system with fairly good confidence. These simulations can be used to verify
the feature-based requirements in the V-process. These simulations are usually
known as feature-in-the-loop simulations.
Figure 3.3 shows a simulation case where a new feature (Feature A) in the
system is simulated with the already validated features (Features B-E). Even though
the system (Features A-E) may have more than one component, in feature-in-the
loop simulations the physical boundaries and interfaces are not considered.
After enough confidence is gained about the feature-in-the-loop simulations for
the new features of the components, component-in-the-loop simulations can be
developed. Component-in-the-loop simulations usually contain all the newly
developed features of the system. Feature-in-the-loop simulations are generally
done separately for easy troubleshooting and for de-coupling of individual
contributions.
Component-in-the-loop simulations contain all the new and carry-over features
as well as the actual electronic and mechanical interface of the system. These
simulations are used in the component-level testing for the engineering develop-
ment process. When prepared properly with the actual system-level interface they
can be directly used in system simulations that include all components of the system
both electric and mechanical.
3 Simulation-Based Engineering 43

Fig. 3.3 Feature-in-the-loop

simulations
Simula

Fig. 3.4
Component-in-the-loop Simula
simulations

In Fig. 3.4, a component-in-the-loop simulation scenario is given. The work is

done all in the simulation environment, however, a notion of the component
physical boundaries exists that forces the interaction of Features A and B through a
common interface as compared to the feature-in-the-loop simulation given in
Fig. 3.3. Generally, this interface is developed as the proposed physical and elec-
trical interface of the component with the rest of the system.
One of the most important challenges of an engineering development process is
to work in a task-based team environment, where different teams are in charge of
different features/components/sub-systems of the project. The development cycle of
different targets can be at different stages at different times, which makes it difficult
to validate functionality with the complete configuration of the system. Testing with
hardware-in-the-loop simulations is an approach developed by engineers to over-
come this problem. In hardware-in-the loop simulations part of the system is
44 M. Cakmakci et al.

emulated using computers using simulations and part of the system is the actual
hardware, which already designed or carried over from the previous version of the
system. In many cases the benefit of the HIL simulations are bidirectional in the
sense that they can both be used for improving the quality of the simulations using
it against the actual hardware or testing a specific prototype hardware for func-
tionality while emulating the rest of the system.
Figure 3.5 shows a hardware-in-the-loop scenario for the system and features
given in Figs. 3.3 and 3.4. This time the actual hardware of the component that
includes Features A and B are run against the rest of the system (Features C–E) all
simulated in the computer environment. It is also important to note that the
preparation of the simulations in the earlier stages help to build successive versions
of the feature, component and hardware-in-the-loop simulations. For example the
physical and electrical based build of the interface in the component level increase
the reuse of the component representation in the hardware-in-the-loop simulations.
A good example of the simulation-based V-process development is the so-called
mode-based controller development process (MBCD) in the automotive industry. In
(Ulsoy et al. 2012), a technical requirements development method is shown for a
specific battery control module example. This example shows how the vehicle
100,000-mile requirement affects specific features (control problems) for a partic-
ular vehicle application. The effect of this requirement and others define the feature
control problem to be solved. The solutions obtained from all of the features
represent the control algorithm for a vehicle.
In the design step, first the control design problem is formulated based on the
given performance requirements and developed mathematical formulation. There
will be more than one control design approach, which will provide a solution for the
control problem. By using analytical methods and/or computer simulations, the best
alternative among these candidate algorithms is selected. If the control problem is
similar to an earlier application, development teams often prefer to start with an

Fig. 3.5
Hardware-in-the-loop
simulations

Simula
3 Simulation-Based Engineering 45

existing control algorithm and try improve the solution by building upon the
existing (and proven) solution.
Then the design is implemented on the actual hardware. During the imple-
mentation phase the objective is to develop a real-time application, which will be
executed in the control module using the desired control algorithm. While devel-
oping the executable code the real-time constraints of the target hardware (i.e., the
controller module) should also be considered. Software implementation of the
algorithm should be matched to the computing resources available and if there are
overruns during the real-time execution simplifications in the algorithm should be
made, or new target hardware should be selected. In today’s modern vehicles,
controller modules also communicate with other controllers via communication
networks. The effects of the loss of this communication with one or more contacts
or the cases of limited communications should be investigated and necessary
modifications should be made.
Testing in the MBCD process starts as early as in the algorithm development
step. By testing the algorithms open-loop (Fig. 3.6a) developers can feed in simple
test vectors and analyze the test output for expected functionality. These simple
algorithms can also be tested against the simpler conceptual vehicle models, which
are available in the earlier stages of the program (Fig. 3.6b). These models are later

Fig. 3.6 Different types of testing in automotive industry

46 M. Cakmakci et al.

Engine
HIL Simulator

Dyno
or

Plant Model

Desktop Computer
Engine Controller Module

Fig. 3.7 Different types of testing in automotive industry (cont’d)

fortified with improvements based on component and vehicle testing data, which
makes them suitable for more complex testing procedures such as module, com-
ponent and vehicle in the loop types of testing.
In the later stages of the vehicle development process, a hardware-in-the-loop
simulation can be run to see the proper operation of the vehicle controller using part
real hardware and part simulations run in the computer environment as shown in
Fig. 3.7.

3.1.2 Source of Models

Models are purposeful abstractions of the real world. With abstraction while certain
aspects of the system are explicitly represented, other aspects are omitted that are
not of concern (Topçu et al. 2016). They can be physical, mathematical, and/or
logical (Sokolowski and Banks 2010). The scaled aircrafts that are used in wind
tunnels are very good examples of physical models. When they are not physical,
models are composed of a series of mathematical equations and/or logical
expressions. These models can be physics-based, data-based, or hybrid (combined).
Physics-based models can be defined as the ones which are essentially mathe-
matical and the governing equations are based on physical principles such as
thermodynamics laws or Newton’s law of motion.
The application of physics-based models in engineering domain is so common.
Since early the days of engineering, Newton’s law of motion has been used for
modeling rigid bodies. Dynamics of machinery is an engineering field that deals
with forces and moments and their effects on the motion. The theory of machines
studies the relative motion of machine elements under the effects of external forces
(Khurmi and Gupta 1976).
Modeling the mechanical behavior as a continuous mass is the topic of con-
tinuum mechanics. It is concerned with the stress in the continuous medium (solids,
liquids or gases) and their deformation or flow (Malvern 1969). Continuous as an
adjective is used to express the approximation that assumes the mass without gaps
3 Simulation-Based Engineering 47

and empty spaces thereby representing the mathematical functions as well as their
derivatives are continuous. This hypothetical medium is called continuum. The
governing physical laws in this case are conservation of mass, momentum, and
energy. These equations will be summarized in Sect. 3.2.6. The motion of viscous
fluids is mostly computed by applying Navier–Stokes equations which encom-
passes time-dependent equations for conservation of mass, momentum, and energy.
Euler equations are well-employed simplification of Navier–Stokes equations
which neglects the effects of viscosity (Schetz and Fuhs 2013). Computational
Fluid Dynamics (CFD) is the area of study which applies numerical methods like
finite difference or finite volume to solve the approximations of these equations.
Physical model of heat has also been built considering it as a fluid inside the matter.
Heat equation is a partial differential equation that concerns the distribution of heat
in material over time (Widder 1976). Solid mechanics deals with the behavior of
solid materials under load. While elasticity is the study of body that retains its
original state after releasing the load, plasticity governs the nonreversible defor-
mation of solid. Euler–Bernoulli beam equation and plate theory are well-applied
simplifications in modeling and simulation of elastic behavior. They both define the
relations between the applied forces and the resulting deflections (Fung 1965).
Finite Element Method (FEM) as will be discussed later in Sect. 3.2.3., is com-
monly employed for approximating partial differential equations within Navier–
Stokes equations, heat equation and Euler–Bernoulli beam equation (Dhatt et al.
2012). It promotes using simple approximation of unknown variables to transform
partial differential equating to algebraic equations.
Data-based models utilize the data that describes the particular aspects of the
system that is subject to modeling. It is also named as empirical modeling since the
model depends on empirical observations rather than mathematical equations
(Sokolowski and Banks 2010). While the computing power as well as the optimized
implementations of finite element analysis and computational fluid dynamics soft-
ware getting better, engineering design optimization of complex systems like air-
crafts or cars requires long lasting simulations which are sometimes unacceptable in
practice (Wang and Shan 2007). Additionally, sometimes it is required to incor-
porate data from the real world into the simulation. Data-based models, or simply
called metamodels approximate computation-intensive functions or real-world data
to analytical models. The modeling process starts with data collection using sam-
pling methods such as fractional refactoring, or Latin hypercube. Then the model is
constructed is a particular method of choice. Polynomial equations, splines,
Multivariate Adaptive Regression Splines (MARS), artificial neural networks are
some of these methods. Model fitting is done with an appropriate approach like
least squares or backpropagation.
Hybrid modeling combines previously mentioned two modeling paradigms.
While a part of a physical process is approximated using data models, the rest of the
physical process is modeled using equations that represent the law of physics. In
modeling and simulation of air vehicles, it is a common practice to develop data
models for the aerodynamics modeling, where the flight dynamics is modeled using
Newton’s laws of motion (Jategaonkar et al. 2004). The aerodynamics data may be
48 M. Cakmakci et al.

collected from the flight experiments, wind tunnel tests or with CFD runs.
Nowadays, the design of complex multidisciplinary systems such as aircrafts,
automobiles and similar is carried out using hybrid models within a Multi-
Disciplinary Design Optimization (MDO) framework (Martins and Lambe 2013).
Such procedures allow designers to incorporate all relevant disciplines simultane-
ously. The optimum of the coupled problem is superior to the design found by
optimally designing each module sequentially, since it can exploit the synergistic
coupling between them. However, this concurrent consideration results in a much
more complex problem. Therefore, systematic structuring, modeling, and approx-
imation tools have to be employed within MDO, which has been applied with
successfully to the design of many commercial products.

3.2 Simulation of Continuum

The term modeling refers to the development of a mathematical representation of a

physical situation whereas simulation refers to the procedure of solving the equa-
tions that resulted from model development (Ashby 1996). With the development
of mathematical models it was possible for scientists to integrate research into
natural phenomena within their investigations. Analysis of these models was only
possible via existing analytical or numerical methods which by then were only
tackled for specific problems. Each of these methods in literature is known by the
great scientist who developed them such as Euler, Newton, and Gauss.
Despite major contributions by various outstanding scientists, the main issues
concerning the theoretical and physical understanding of the equations in contin-
uum mechanics are still being worked on. Continuum mechanics has changed
dramatically since the late nineteenth century, so that theoretical studies are now
coined with numerical experimentation and simulation. Furthermore, progress in
the computational speed and power allowed researchers to develop mathematical
models for much more complex physical problems some of which will be discussed
in the multi-scale and multi-physics sections. After the invention of calculus, many
advanced PDE’s were introduced to describe the physics of systems from different
disciplines such as solid mechanics, fluid mechanics, and elastodynamics.
Important contributions were initially made by Euler, Lagrange, and Cauchy and
these were followed by the application of PDE’s to describe the physics of elec-
tromagnetic (EM) theory by Maxwell, Heaviside, and Hertz, and finally to quantum
mechanics with major theoretical work by Schrodinger. These equations are
descriptive of the time evolution and relationship of various fields in a three
dimensional space. Major efforts of continuum simulation in the areas of Solid
Mechanics, Fluid Mechanics, and Electromagnetics will be described in the next
sections to follow.
The introduction of efficient and powerful platforms enabled researchers to solve
the constitutive laws of continuum in mechanics in combination with the laws of
conservation of mass, energy, and momentum. The same is valid for other fields
3 Simulation-Based Engineering 49

including EM. Some of the most popular methods used for this purpose are the
Finite Element Method (FEM), Finite Volume Methods (FVM), Finite Difference
Methods (FDM), and Boundary Element Methods (BEM). These methods are
applied to the simulation of matter in all forms, i.e., solids, liquid, and gas, based on
a major assumption of continuum media, thus Computational Mechanics of
Continua. Namely, the term continuum describes the nonseparability of the con-
sidered domain and validity of continuity between any points in the domain so that
differentiation is possible. Therefore, continuity between elements in any
continuum-based numerical technique is maintained as well. Unlike analytical exact
solutions of differential equations, which allow the solution at every point, the
numerical solution is only calculated at chosen finite number of nodes, yielding in
turn a reduction in complexity of the system. Well-known methods to conduct
continuum-based simulation are described in the next section.

3.2.1 Finite-Difference Method

One of the earliest and widely used numerical method for solving PDE’s within
continuum mechanics is the Finite-Difference Method (FDM). The main idea of
FDM is based on replacing the differential terms with respective to the spatial
coordinates with the so-called finite differences over small enough distances based
on the Taylor’s series approximation. For that purpose, the domain of interest needs
first to be discretizedinto vertical and horizontally located nodes, on which finite
differences are defined. Several finite difference integration schemes exist known as
forward, backward, and central difference schemes. It is worth noting that the FDM
is equally applicable to time differentiations. As a result of discretizing the domain
into nodes, a system of algebraic equations in terms unknowns at the chosen nodes
are constructed. Each algebraic equation belonging to its corresponding node is
expressed as a combination of function values at its own node and its neighboring
nodes. Next step is to impose boundary conditions, which leads to the solution step
of the equation system using either direct or iterative solution methods. Finally,
unknowns at each node are solved. This solution is only an approximate solution
since the finite differences are first-order approximations of the partial derivatives.
The FDM when compared with the FEM or BEM allows for a direct discretization
of the equations and does not rely on the use of interpolation functions. Therefore, it
is one of the most direct and intuitive techniques that exist for the solution of PDEs.
Moreover, for material nonlinearities, the FDM proves to be favorable as it allows
their simulation without the need of iterative techniques. However, it suffers from
relying on regular noded discretization scheme which makes modeling of irregular
geometries a challenging task. This also results in difficulties when heterogeneous
material compositions and unusual boundary conditions are present. However, the
FDM has been generalized to overcome related shortcomings through methods
based on irregular node/grid structures with methods such as irregular quadrilateral,
triangular, and Voronoi grids.
50 M. Cakmakci et al.

3.2.2 Finite Volume Method

The Finite Volume Method is similar to the FDM method and evolved as its
successor to solve PDE’s with one major difference: these differential equations are
expressed in integral form. Its formulation leads to the concept of finite volumes,
which essentially correspond to volumes around and encompassing each node in a
mesh. Similar to the FDM, algebraic equations of unknowns at nodes are built by
replacing the integrals and by considering boundary and initial conditions. Thereby,
the system of equations to be solved is constructed. The FVM, similar to the FDM
has certain advantages such as allowing the usage of irregular unstructured mesh
and modeling capabilities of nonhomogeneous material compositions.

3.2.3 Finite Element Method

The Finite Element Method (FEM) was introduced in the 1960s as an alternative
method to FDM for the numerical solution of stress concentration. More impor-
tantly, it is the first numerical solution method which was capable of dealing with
complexities such as nonlinearities, nonhomogeneous materials, complex geome-
tries, and sophisticated boundary conditions. As a result, FEM was soon recognized
as the most popular numerical method in continuum mechanics, mainly so because
unlike FDM, it allowed for nonuniform discretization. The method was found more
extensive and used a decade later with the theoretical developments made by Bathe
(2006) and Zienkiewicz and Taylor (2005). Many researchers have contributed to
the development of the method which is by far the most favorite method for the
approximate solution of many sophisticated continuum mechanics problems of
dynamic, anisotropic, and inelastic behavior. It is a generic numerical solution
technique for boundary value problems coming from various disciplines. The main
principle rests on the idea of dividing the problem domain into smaller subregions
(areas or volumes) called finite elements. This is followed by typical steps of
defining local element approximations, performing assembly of finite elements and
ultimately solving the resulting global matrix equation. More specifically, the
unknown function (e.g., displacement field, temperature field, electric field, velocity
and pressure fields) is approximated via trial/interpolation functions of the nodal
values (or edge unknowns in EM problems) using polynomial functions. Numerical
integration is performed in each element using Gauss quadrature points. After
assembly, the algebraic global system of equations is obtained. Because of con-
tinuum assumptions, standard FEM methods cannot be directly and efficiently
applied to discontinuum problems involving cracks, damage-induced discontinu-
ities or singularities and failure analysis.
In addition to the well-known superiority of the FEM which is well suited for
complex analysis of systems composed of heterogeneous materials and irregular
geometries owing to the possibility of using an irregular mesh, it also proved to be
3 Simulation-Based Engineering 51

an appropriate tool for modeling various nonlinear geometries and inelastic material
behavior and nowadays material hardening and softening. Moreover, it has the
capability of representing geometric nonlinearities, contact mechanisms,
fluid-structure interaction, multi-scales, etc., as will be discussed in separate sec-
tions below. Therefore, the FEM will stand out as the mostly used and diverse
numerical method in continuum mechanics.

3.2.4 Meshless Methods

The bottleneck in applying FEM to complex engineering problems with intricate

geometries, unusual material properties, and complex boundary conditions is the
mesh generation process, which usually in 3D problems is an extremely challenging
task mostly comparable to the problem solution itself. Another disadvantage of
FEM relates to numerical instability due to a distorted mesh. Both of these problems
can be avoided by another class of methods known as ‘meshless methods’, which as
the name implies does not rely on elements but interpolation functions are gener-
ated from neighboring nodes within a domain of influence. More specifically, nodes
are created across the domain without the need of a fixed element topology defi-
nition. As a result, the interpolation functions obtained are no longer polynomial
functions leading to more difficult numerical integration when compared with the
FEM where Gauss integration points are used. Moreover, meshless methods suffer
from increased computational requirements but do not rely on standard mesh
generators and are able to easily represent more complicated geometries. In liter-
ature, methods such as smoothed particle hydrodynamics, diffuse element method
(DEM), element-free Galerkin method, reproducing kernel particle methods,
moving least squares reproducing kernel method, hp-cloud method, the method of
finite spheres, and finite point method stand out.

3.2.5 Multi-scale Methods

All products whether man-made or natural are composed of multiple scales. Taking
an example from the aeronautical industry, the Airbus A380 consists of many
thousands of structural components and many more sub-structural details.
Similarly, its fuselage consists of 750,000 holes and cutouts with different structural
and material scales. When viewed at the roughest material scale, fuselage com-
posites’ part consists of woven/textile composite and laminate scales; at the inter-
mediate scale, it is composed of a tow or yarn, which consists of a bundle of fibers.
When looked at a more discrete scale, including atomistic and ab initio scales, the
aircraft’s metal part consists of a polycrystalline scale, a single crystal scale, a
discrete dislocation scale, and also time atomistic and ab initio scales.
52 M. Cakmakci et al.

Computations and simulations in the aforementioned multi-scales have been

identified as areas of utmost importance to advance the future in nanotechnology.
One of the obvious fundamental challenges associated with such a multi-scale
approach relates to the increased uncertainty and complexity introduced by these
finer scales. However, the application of any multi-scale approach has to be care-
fully evaluated. For instance, considering metal matrix composites with fibers
arranged in a periodic fashion, finer scales could prove useful because the bulk
material typically does not obey normality rules, and the development of a phe-
nomenological coarse-scale constitutive model would be extremely difficult. This
would also allow a better understanding of each phase and the overall material
response could be extracted from its fine-scale constituents via homogenization
techniques. However, for brittle ceramic matrix composites, with microcracks that
exist in a random distribution and complex interface properties, are difficult to
characterize, a multi-scale approach would not be an appropriate alternative.
There are two main categories of multi-scale approaches in literature, namely,
hierarchical or concurrent. In the former approach, the fine-scale response is
idealized/approximated and its overall/average response is integrated into the coarse
scale. In the latter approach, fine and coarse-scale resolutions are simultaneously
employed in different portions of the same problem domain, and the exchange of
information occurs through the interface. The sub-domains which present them-
selves at different scale resolutions can be either overlapping or disjoint.
Various hierarchical multi-scale methods have been labeled by different names,
including upscaling, coarse-graining, homogenization, or simply multi-scale
methods. There are also subcategories of the above definitions, such as system-
atic upscaling, operator upscaling, variational multi-scale, computational homoge-
nization, multigrid homogenization, numerical homogenization, numerical
upscaling, and computational coarse-graining, just to mention a few. Moreover,
different definitions are used to indicate various scales. If the structure exists in two
scales, however, the fine scale is often referred to as a micro-scale, unresolvable
scale, atomistic scale, or discrete scale; the coarse scale is often defined as the
macroscale, resolvable scale, component scale, or continuum scale. For more than
two scales, the additional scales may be termed as mesoscales.
One alternative approach to the homogenization of artificial structures avoiding
the limitations associated with earlier analytical homogenization models (Milton
2002) is the theory of a mathematical homogenization approach. It is based on the
asymptotic expansion also known as two-scale homogenization, which is a well
established concept in the theory of PDEs with rapidly oscillating periodic coeffi-
cients (Bensoussan et al. 1978). Its main advantage is that the method is generalized
enough and unlike analytical techniques, can handle unit cells with inclusions of
arbitrary geometry and any number of phases with no additional computational
cost. Also, instead of formulating the problem as an eigenvalue problem, two-scale
homogenization works directly on the original form of the governing equations and
is therefore able to result in expressions valid for effective constitutive tensors.
In their studies, (El-Kahlout and Kiziltas 2011) further developed this approach
by applying two-scale homogenization method to Maxwell’s equations and
3 Simulation-Based Engineering 53

extracting the effective parameters of periodic dielectric and magnetic materials,

that can be in their most generalized form lossy and are made of inclusions with
arbitrary shapes and multi-phase material constituents. The numerical solution of
the resulting PDE is carried out using a commercial FEA based solver, namely
COMSOL Multiphysics, where the effective tensors are evaluated at a single fre-
quency, for both isotropic and anisotropic effective material tensors with isotropic
constituents. This is the first study where numerical material model based on
two-scale homogenization is used to synthesize the microstructure of EM material
with desired material matrices using formal design techniques such as topology
optimization. Results of this design study are shown in Fig. 3.8 which was also
fabricated (El-Kahlout and Kiziltas 2011) using novel Dry Powder Deposition
techniques as demonstrated in Fig. 3.9.
Similar to EM materials, most heterogeneous materials, such as composites,
polycrystals, and soils consist of constituents/phases with clear-cut boundaries that
display different mechanical and transport properties. The use of homogenization of
continuum allows a better understanding of the physical governing equations of
individual phases, including their geometry and constitutive equations at the
fine-scale phases, or at least a better grasp than at the coarse-scale phases. Put in
other way, the process of homogenization provides a mathematical means by which
coarse-scale equations can be deduced from well-defined fine-scale equations.
Moreover it allows the determination of heterogeneous material behavior, at least
theoretically without the need of testing, which is usually a very expensive
endeavor. Also, through homogenization, one can estimate the full multiaxial
properties and responses of heterogeneous materials, which present themselves as
anisotropic materials and are most of the time extremely difficult to measure
experimentally. In addition to describing the overall behavior of heterogeneous
materials, the act of homogenization leads to local fields via the process known as
downscaling given coarse-scale fields, phase properties, and phase geometries. This
information is of critical importance in understanding and describing material
damage and failure.

140 140

2 120 5 120

100 10 100
4
80 15 80
6
60 20 60
8 40 25 40
20 20
2 4 6 8 5 10 15 20 25

Fig. 3.8 Optimal material distribution (dielectric ranges from 20 to 140) of designed unit cell
(left) and array (right) for a desired permittivity tensor of e = [45.0;’0.70] using mathematical
homogenization and topology optimization.[reproduced courtesy of The Electromagnetics
Academy]
54 M. Cakmakci et al.

Fig. 3.9 Automated fabrication of design in Fig. 3.8 using dispensing machine within DPD in
action (left) and resulting desired deposited substrate (right). [Reproduced courtesy of The
Electromagnetics Academy]

3.2.6 Solid Mechanics

Two main branches exist within the application of the principles of mechanics to
bulk matter: the mechanics of solids and fluids. When viewed from a global per-
spective, the common subject is that of continuum mechanics. More specifically,
continuum mechanics conceives the useful model of matter as continuously
divisible, and does not make any reference to its discrete structure at microscale,
which is well below those scales of the phenomenon of interest. Solid mechanics is
concerned with stresses, deformation, and failure of structures and solid matter.
A material is called a solid and not a fluid if it is able to support significant amount
of shear force over a certain time period of a natural process or technological
application of interest.
The main equations of continuum physics can be presented by separating them
into global and local laws. The former serve as the foundations of continuous media
theory and are summarized here (Muntean 2015). In all these formulations, X′
(t) denotes arbitrary configuration of partial volume B′ of B. More specifically,
global balance laws for the five major conservation principles are presented here for
mass, linear and angular momentum, energy, and entropy.

Mass:
The conservation of mass is expressed in its most general form as

d
mðX0 ðtÞ; tÞ ¼ 0 ð3:1Þ
dt

for all X′(t)  X(t), where m(t) stands for the total mass in X′(t), i.e.,

mðX0 ðtÞ; tÞ ¼ dlm ¼ qdx ð3:2Þ

Z Z

X0 ðtÞ X0 ðt Þ
3 Simulation-Based Engineering 55

with q(t, x) denoting the density. Assuming that there is no internal mass pro-
duction, (3.2) states that the total mass of any material partial volume is conserved.
Linear Momentum:
The conservation of linear momentum or balance of forces is expressed in its
most general form as: For every part X′(t)  X(t) we have

d
‘ðX0 ðtÞ; tÞ ¼ F ð3:3Þ
dt

where the linear momentum can be defined via

‘ðX0 ðtÞ; tÞ ¼ v qdx ð3:4Þ

Z

X0 ðtÞ

More specifically, the time rate of change of total linear momentum of l in B′ is

equal to the force F exerted on B′. The force F consists of the contribution of the
internal body forces per unit of volume qfb and contact or surface forces per unit of
area t acting on the boundary ∂B′ of B′. Here, t is the stress vector or traction.

Angular Momentum and Moment of Momentum:

The conservation of angular momentum or balance of moments is expressed in its
most general form as
For every part X′(t)  X(t), we have

d
aðX0 ðtÞ; tÞ ¼ M ð3:5Þ
dt

where the angular momentum can be defined via

aðX0 ðtÞ; tÞ ¼ x  v q dx ð3:6Þ

Z

X0 ðtÞ

More specifically, the time rate of change of total angular momentum a of B′ is

equal to the moment M of the force F exerted on B′.
Energy:

The conservation of energy balance is expressed in its most general form as: The
time rate of change of the total energy within B′ which is composed of the kinetic
energy K and internal energy E and is equal to the rate of work, say P, done by both
the body force and the contact force, plus the heat supply Q from internal heat
production and heat fluxes across the boundary of B′. So for every part X′(t) 
(t), this can be written as
56 M. Cakmakci et al.

d
ðK ðtÞ þ EðtÞÞ ¼ PðtÞ þ QðtÞ ð3:7Þ
dt

where
jvj2 jvj2
K ðt Þ ¼ dlm ¼ q dx; ð3:8Þ
Z Z

X0 ðtÞ
2 X0 ðtÞ
2

E ðt Þ ¼ e dlm ¼ e q dx; ð3:9Þ

Z Z

X0 ðt Þ X0 ðtÞ

PðtÞ ¼
Z
v
ðTnÞdr ¼
Z
~
f
~
v qdx; ð3:10Þ
X0 ðtÞ X0 ðtÞ

Q ðt Þ ¼ fHeat q dx þ q
ndr: ð3:11Þ

Z Z

X0 ðt Þ X0 ðt Þ

In Eq. (3.9), e represents the inner energy density. The first term in Q(t) accounts
for the heat source. The measure lm in the equations of K(t) and E(t) corresponds to
the mass measure associated with the material body B.
Entropy:
The entropy increase within B is greater than or equal to the internal entropy
supply, i.e., internal heat source over h, which is the absolute temperature, plus the
entropy flux across the boundary of B′, which can be expressed as follows:
0 1
d@ Z fHeat q
n
s dlm A  ð3:12Þ
Z Z
dlm dr:
dt X0 ðtÞ X0 ðt Þ
h 0 h
@X ðtÞ

Here s represents the entropy density. It is explicitly noted that all conservation
laws are in term of extensive quantities. More specifically, global balance laws can
only be written in terms of extensive quantities. However, the intensive quantities
are related to local balance laws expressed in terms of PDEs and inequalities as well
as boundary conditions and can be derived based on global laws of the preceding
section (Muntean 2015).

3.2.7 Fluid Mechanics

Various theories govern the physics of fluid mechanics and different methods are
proposed and used in literature to provide numerical solutions/simulations primarily
depending on the spatial and temporal scale of the phenomenon. Instead of going
3 Simulation-Based Engineering 57

Fig. 3.10 Typical numerical methods used in fluid mechanics based on temporal and spatial scale

into detail with all methods, these theories and typical numerical methods employed
according to the temporal and spatial scales are summarized in Fig. 3.10. As
depicted in the graph, continuum mechanics prevails for above microscale and
below tens of meters with a time scale between 1 s and hours. When the continuum
assumption breaks down, the fluid has to be described by an atomistic point of
view, such as the molecular dynamics as a microscale method or statistical rules
govern the molecular group behavior, i.e., kinetic theories as mesoscopic methods
for larger scales. On the spatial and time scale limit, if the characteristic length is
smaller than 1 nm or the characteristic time is shorter than 1 fs, the quantum effect
may not be negligible for the system of interest and quantum mechanics has to be
brought into describe the transport. In fact, modeling at a smaller scale may present
a more accurate description of the problem, but is likely to cause a much higher
computational cost. Therefore, as always in numerical simulations, in engineering
an appropriate tradeoff is considered when trying to determine in an accurate and
fast way, the fluid behavior of interest.
Despite the emergence of high-speed platforms and the advances in efficient and
accurate numerical methods, some computational fluid dynamics (CFD) problems
still present themselves as challenging problems for the practical solution via
numerical simulation techniques. For example, NASA has recently modified its
aerospace design codes for earth science applications, thereby speeding up super-
computer simulations of hurricane formation (Kazachkov and Kalion 2002). An
example of such a CFD simulation using a 512-processor supercomputer is referred
to in (Kazachkov and Kalion 2002). More specifically, actual data from a variety of
different sources and climate models were integrated to generate high fidelity
simulations so as to reproduce a hurricane forming in the Gulf of Mexico. As a
result, engineers were able to simulate the formation and movement of a hurricane.
However, the weather forecast of global earth based on CFD atmospheric and ocean
58 M. Cakmakci et al.

simulation is still a challenging problem, and calls for even larger amount of
computing power and more accurate data. Overall, this is a multi-phase CFD
problem with very complex geometry and dynamic boundary conditions.

3.2.8 Electromagnetics

In electromagnetics (EM), matrix systems with a few millions of unknowns known

as dense matrix systems have been solved numerically for ten years now. Today,
the number of unknowns that can be solved via simulations is on the order of a
billion of unknowns (Gurel and Ergul 2007). This impressive improvement is
attributed to the synergistic progress between hardware and algorithm design. It is
also noted that for sparse matrix systems resulting from specific simple EM
problems within electrostatic and magnetostatics, even larger scale problems can be
addressed such as the World-record algorithm from Jülich calculating over three
trillion particles (World-Record Algorithm from Jülich Calculates Over Three
Trillion Particles—Research in Germany 2011).
As stated earlier, the physical response of many fields including EM, can be
analyzed via differential equations. Therefore, PDE’s were used for more than four
centuries and continue to set the standard for modeling the physics of different
media today. There are three main groups of differential equations, namely
hyperbolic, parabolic, and elliptic, which describe fields with various physics. The
Laplace equation or the Poisson equation given below is a well-known example of a
generalized simple elliptic PDE. These are encountered in the numerical modeling
of EM problems in the static regime and various transport problems. They are
known for characterizing fields or potentials associated with no singularities distant
from the source location, or equivalently these equations are differentiable functions
and therefore do not allow for any singularity propagation.

qðrÞ
r2 uðrÞ ¼ : ð3:13Þ
e

Typical examples of parabolic equations, the second group of PDE’s, are the
Schrodinger and diffusion equations. These equations are characterized by their first
time derivative and second space derivative. These are fundamental equations in
quantum mechanics and heat transfer as well as low-frequency EM propagation in
conductive media, respectively. A diffusion equation in standard form is

1 @
r2 uðrÞ uðrÞ ¼ 0: ð3:14Þ
c s @t

The third class PDE’s refers to hyperbolic equations, and an example belonging
to this group is the wave equation. It has second-order space and time derivatives.
3 Simulation-Based Engineering 59

1 @2
r2 uðrÞ uðrÞ ¼ 0: ð3:15Þ
c2 @t2

Solution of differential equations as earlier denoted are carried out by three major
methods: a subspace projection method (e.g., FEM), the FDM, and the pseudo-
spectral method. Various basis/interpolation functions are introduced to fit the
unknown field (Chew 1995) in the subspace projection method. It covers a sub-
space of the larger space that the field is defined over due to the finite characteristics
of basis functions. Thereby, the PDE is easily converted to a time dependent
ordinary differential equation. For the solution of the equation via time stepping or
marching, the derivatives can further be approximated using finite difference or the
subspace projection method. Or as an alternative, time domain Fourier transform
can be used to remove the time derivatives resulting in a matrix equation to be
solved via iterative or inversion techniques.
A major alternative exists to the numerical solution of the governing Maxwell’s
equations in EM expressed in PDE. Specifically, initially a point source response
called Green’s function can be introduced. Based on linear superposition and as a
result of an arbitrarily distributed source, the unknown field is obtained via spatial
convolution of the distributed source expressed via Green’s function. This corre-
sponds to the equivalence principle (Harrington 2001) which allows the field in a
given region to be expressed as Green’s operator acting on the sources. Hence, the
resulting equations are of integral equation type (IE). When compared with PDE, IE
have an important advantage where the EM unknowns correspond to only surface
unknowns, or to volume unknowns that occupy only a spatial finite region.
Therefore, the number of unknowns in the IE formulation may be much less than
those in the PDE formulation. More importantly, this IE formulation leads to the
automatic satisfaction of the radiation condition if a suitable Green’s function is
chosen. However, in the PDE formulation absorbing boundary conditions or the
so-known boundary integral equations replace the radiation condition. Additionally,
using the subspace projection method (Harrington 2001), these IE can be converted
into matrix equations. Equivalently, operators of the integral are replaced with
matrix operators. However, the matrix representation of Green’s operator corre-
sponds to a matrix system which is dense because of its non-local nature. Hence,
the computational storage and operations such as matrix vector products with that
type of a matrix system can be computationally expensive. In literature, some
methods have been developed to overcome these expensive matrix solutions. These
include fast Fourier transform based methods, fast-multipole-based methods,
rank-reduction methods, the nested equivalence principle algorithm, recursive
algorithms, etc. (Weng et al. 2001).
As a final class of numerical techniques for EM radiation and scattering prob-
lems, hybridized versions of the two main classes combining their advantages have
been developed. The FE-Boundary Integral (BI) method is one of the most pow-
erful techniques belonging to this class. More specifically, it offers the flexibility of
the FEM to analyze structures with highly complex geometrical and material details
60 M. Cakmakci et al.

Fig. 3.11 Design results of novel material distributions of a patch antenna via integration of
FE-BI method and topology optimization (Kiziltas et al. 2003)

but at the same time imposes a rigorous boundary condition via the use of the BI
formulation. This tool’s efficient and accurate analysis capability has allowed
researchers to conduct numerous designs (Volakis et al. 2006).
It is especially noted that these efficient and accurate codes allowed for the first
metamaterial-based antenna design using topology optimization based techniques
as shown in Fig. 3.11 (Kiziltas et al. 2003). The design developed from scratch as
shown in Fig. 3.11 was based on 5 individually textured layers which were also
fabricated and measured. The agreement between measurements and calculations is
truly impressive for the complex dielectric design. Above all, the threefold
improvement in bandwidth is a clear demonstration of the remarkable potential of
efficient and accurate numerical techniques in delivering novel designs not only in
EM but also in other engineering disciplines.
Multi-scale problems as discussed in Sect. 3.2.5, present themselves in circuits,
packages, and chips at various levels of complexity. Similarly, they exist also in
antennas on complex platforms, in nano-optics and nanolithography applications.
Therefore, multi-scale solutions of problems are critical for many applications.
Similar to other applications, the size evaluation of the EM multi-scale problem is
of great importance. More specifically, one needs to evaluate the multi-scale
structures relative to the wavelength to determine which physics of the three to
apply for their solution: circuit physics, wave physics, or optics physics. Avoidance
or identification of ill-conditioned numerical systems plays a great role in the
effective solution of multi-scale EM problems.
It is finally noted, that one of the biggest challenges today in the numerical
solutions of EM problems is the model size of realistic problems which deems
high-performance computing a vital necessity. Significant speedups have been
achieved by hardware scaling and additional efforts have resulted in three main
types of HPC platforms: (1) supercomputers, (2) computer clusters, and (3) cloud
3 Simulation-Based Engineering 61

computing. It is quite evident that, computational EM and large-scale computing

will continue to evolve given that these are indispensable tools for EM analysis and
design. Not only will it allow for efficient and practical performance evaluation and
novel designs but it is expected to continue the enhancement of our thorough
understanding of the physics within highly complex systems.

3.2.9 Multi-physics Methods

Many realistic problems present themselves as very complex problems due to their
multi-physics nature. Scientists and engineers from various fields have been
working on the combination of different numerical techniques with the goal of
addressing these elaborate physical processes, such as the transition from contin-
uum to discontinuum (e.g., fracture processes) or the interaction of multi-phases of
matter (e.g., hydrofracture processes). As a result, a new class of numerical methods
called hybrid/multi-physics methods evolved. It is due the developments in
high-performance computing and computational science and computer hardware
that this group of methods evolved. Major examples are: Combined Finite-Discrete
Element Method (F-DEM), Hybrid Lattice Boltzmann-FEM, Lattice
Boltzmann-DEM, etc. Areas of interest include algorithms and novel solutions for:
– Coupling of FEM and DEM simulations
– Coupling of FEM and/or DEM with CFD solvers
– Coupling of different solvers of continuum mechanics, e.g., FEM-FVM.
– Coupling of continuum and discontinuum mechanics solvers, e.g., FEM-DEM,
FEM-MPM, FEM-LBM, etc.
– Coupling of solid and fluid mechanics solvers, e.g., FEM-LBM, FEM-FVM,
etc.
– Coupling of discontinuum mechanics solvers, e.g., DEM-SPH, DEM-LBM, etc.
– Coupling of solvers for different scales, e.g., coupling of FEM-DEM.

3.3 Simulation of Machinery

In many simulation studies developers represent the components of the target

system based on their dominant energy-based properties. Although various linear
and nonlinear extensions exist, basic energy-based properties can fundamentally be
given as inertia, storage (spring), and dissipation (damping). For example, in many
engineering simulations a bearing is represented as a damping element and an axis
slider in a manufacturing system is represented as inertia. It is also noted that a
bearing component brings a negligible rotational inertia to the system as well as a
slider which is somewhat flexible and its dimensions may change very slightly
under heavy operating conditions. But these are not considered as dominant
properties for these components.
62 M. Cakmakci et al.

The approach of representing complex and spatially distributed physical systems

based on their dominant energy properties is known as the “lumped parameter
modeling”, implying the dominant energy-based characteristic(s) of a component
are represented by using specific and predetermined elements (Karnopp et al. 2000).
The use of lumped parameter systems approach results in a more structured
approach of developing simulations for complex engineering systems. Using
energy rather than other physical features (force, current, etc.) also makes it possible
to use this approach in multi-domain systems.

3.3.1 Single Degree of Freedom Systems

Generally, in engineering, the classification of lumped parameter systems is given

based on the number of the inertia elements. In many cases especially for
mechanical systems, the freedom of motion of the component represented as inertia
is important. For example, if a component can move in both x and y axis and/or can
also rotate about z axis, these motion properties are all represented as separate
inertial elements and flow variables.
Single degree of freedom systems are systems represented with one inertial
element and have one variable governed by fundamental physical equations.
A good example for a single degree of freedom system is the longitudinal motion
simulation of vehicles as shown in Fig. 3.12.
The system represented in this figure is given in (3.17) as a mathematical
relationship (i.e., model) based on Newton’s second law, where wheel traction
forces, F, are the input and the vehicle acceleration, ax, is the output.

max ¼ ðW=gÞax ¼ Fxr þ Fxf W sinH Rxr Rxf DA þ Rhx ð3:17Þ

This mathematical model is straightforward to apply in simulation environments

such as MATLAB/Simulink. The single degree of freedom longitudinal simulation
can be used in basic fuel economy and traction (acceleration/braking) studies as

Fig. 3.12 Longitudinal

motion of a vehicle
3 Simulation-Based Engineering 63

reported in Rajamani et al. (2000), Ulsoy et al. (2012). However for many vehicle
engineering studies such as axle-based traction control (Cakmakci et al. 2011;
Dokuyucu and Cakmakci 2016) more complicated representations (i.e., higher
fidelity simulations) that are also suitable for V-process development model dis-
cussed in Sect. 3.1 is needed.

3.3.2 Multi-Degree of Freedom Systems

One way to improve the fidelity of the simulations is to increase the degree of
freedom of its underlying mathematical model. This can be done by increasing the
number of flow variables representing the inertia element, or adding more inertia
elements to the system simulation. As an example of increasing the fidelity of the
model by adding new flow variables to the inertia representing is the half-car model
for vertical motion given in Fig. 3.13.
In Fig. 3.13, the vertical motion of a vehicle is represented with two degrees of
freedom (translation and rotation about the center of mass) rather than only the
vertical motion of the mass of the vehicle. A detailed mathematical model describing
this system can be found in (“Automotive Suspension—MATLAB Simulink
Example,” n.d.) 2017, using road elevations, q, as input and vertical movement of
the center of mass, z, and body rotation, h as outputs. With this representation in
simulations, the vertical motion of the occupant area can be studied as well as the
wheel based vertical road force, which is critical for traction control studies such as
wheel-based braking, and acceleration with so-called load transfer.
Another way of increasing the content is to increase the number of inertial
elements considered in simulations. In this case, rather than using a single system
boundary, where all of the components were lumped together before, can be broken
into components and their relative interaction can be studied.
A good example for this kind of situation is the quarter car model shown in
Fig. 3.14. In this model, a quarter of vehicle vertical dynamics is studied using
quarter of the mass of the vehicle with the suspension system represented by ks and
cs, f, tire parameters, kus and cus and vertical motion variables, z.
The mathematical equations representing this simulation is given in (3.18) based
on Newton’s second law:

ms€zs þ cs ð_zs z_ us Þ þ ks ðzs zus Þ ¼ f

ð3:18Þ
mus€zus þ cs ð_zus z_ s Þ þ ks ðzus zs Þ þ cus ð_zus z_ 0 Þ þ kus ðzus z0 Þ ¼ f

Fig. 3.13 Vehicle vertical

dynamics θ Vehicle dynamics
z

K1 C1 K2 C2
Road
excitation q1 q2
64 M. Cakmakci et al.

Fig. 3.14 Quarter car model

for vertical motion
simulations

It is important to note that by adding a new inertial element representing the

mass of the wheel hub and suspension frame, the stroke motion of the suspension
can be studied including the effects from the tire and the vehicle inertia which is not
possible with the system given previously in Fig. 3.13.
The quarter car model given in (3.18) is an example of a multicomponent
simulation in single axis of motion. More complicated models can be used to study
wheel based multi degree of motion as given in Figs. 3.15 and 3.16 respectively for
vehicle vertical motion (Rajamani and Hedrick 1995) and full vehicle motion.
In engineering simulation studies, as the number of components and the ele-
ments representing these components increase, the number of mathematical equa-
tions representing used in these simulations also increase. Therefore, the
appropriate simulation content should be chosen to do the correct analysis with
optimal computation time. For example, the full car model given in Fig. 3.16 is
executed by solving 18 nonlinear equations per simulation step size as compared to
the longitudinal model given in Fig. 3.12 contains only one ordinary differential
equation. Both of the models can be used for fuel economy studies and both of them
will contain uncertainties in inertia, spring and storage parameters.

Fig. 3.15 Standard half-car

model
3 Simulation-Based Engineering 65

Fig. 3.16 Full car (18DOF) Axle

model Wid
th th
Leng

W. H
ub

W. H
ub

on
a
l ev
a dE
Ro

3.4 Simulation of Multi-domain Systems

In modern engineering systems, components that operate primarily in different

domains (such as mechanical, electrical, digital) work together to complete com-
plicated tasks. Therefore, a realistic simulation of the system should include the
elements from different domains such as mechanical parts, power electronics to
energize these mechanical parts, and digital components to monitor and/or control
the operation.
Two representative cases of multi-domain simulations improving the perfor-
mance can be given as the multi-domain simulations of electromechanical systems
with controllers and online simulations to improve smart mechatronic/robotic
systems.

3.4.1 Control Systems

Generally, algorithms in control systems are designed to have a certain dynamic

behavior which can be represented in terms of a transfer function (in Laplace
Domain) or a state space model (Ogata 1990). When these algorithms are actually
implemented in actual systems their performance shows variations (usually degra-
dations) due to the effect of execution in digital medium. These variations are due to
the digitization of the algorithm, lack of realistic representation of the control system
hardware and the effect of a communication in common medium such as networks.
Many control algorithms are developed by using a frequency or a time
domain-based structured method based on their dynamic properties as discussed in
many sources in control literature (Chen 1995; Ogata 1995). Once the development
is finished, the resulting output is a fractional function that represents the
input/output relationship of the control algorithm called controller transfer function,
66 M. Cakmakci et al.

C(s), where s is the Laplace variable. The dynamic controller relationship can also
be represented by a matrix equation pair generally given in the form x_ ¼
Ax þ Bu; q ¼ Cx þ Du where u is the controller input, q is the controller input, and
x is the controller states. These representations both imply a continuous system
where calculations or events take place instantaneously using a specific order.
However, when implemented in a real-time control system, algorithm computations
take certain amount of time to finish before an updated command can be issued. For
many systems with fast dynamics, the effect of implementation generates a defi-
ciency in performance since the optimal performance was designed for a medium
where events take place instantaneously.
More realistic and predictable results can be obtained by using discrete control
systems that take into consideration of the digital timing in their formulation
(Franklin et al. 2009; Ogata 1995). Algorithms can be designed and simulated as
digital controllers using adaptations of the continuous methods. Alternatively,
continuous controller functions can be digitized afterwards using simple methods.
For example by using a direct conversion approach, a controller transfer function C
(s) can be converted to its discrete version by replacing the s operator with
ðz 1Þ=Tz using a backward difference transformation. In this recipe z is the dis-
crete variable and T is the sampling period.
Another important aspect of implementation of an algorithm in the digital world
is the effect of quantization. Real numbers can have infinitely large digits during
calculations, however for computers, it is more practical and maintainable to do
operations in chunks of bits causing the calculations to take place in limited digits,
which cause round off errors [Franklin et al.]. The effects of digitization and
quantization can both be included in simulations to predict possible controller
performance degradation in engineering systems.
Another important aspect in implementation of the control systemsis to include
the hardware related properties such as un-modeled sensor/actuator dynamics (S
(s) and A(s) respectively) and the effect of sampling as shown in Fig. 3.17.
In many controller development activities, the controller is designed based on
the plant dynamics P(s) only without including the sensor (S(s)) and actuator (A(s))
dynamics. The dynamic response effects provided by actuators and sensors can be
included in simulations by using time delays, noise, and offsets. The effect of digital
to analog conversion in the actuator is modeled as a zero-order hold (ZOH) element
that keeps the value of the actuator output constant for one time step. This element
also represents the fact that the actuator has internal dynamics and cannot change its
output instantaneously. A sampler element is used at the sensor to represent the

Fig. 3.17 Feedback system Actuator A(s) Sensor S(s)

Plant
with device boundaries and
sampling u T ys
r C(z) ZOH S(s)P(s)A(s)
3 Simulation-Based Engineering 67

Fig. 3.18 Dedicated digital (a) r u Controlled

C A System (P) S
communications (a) versus Dedicated
networks (b, c) y Digital Line

(b) C A
Controlled
S
System (P)
u r,s u y
Network

(c) Controlled
System

A1 S1 A2 S2 A3 S3

C1 C2 C3

Network

analog to digital sampling with rate T. This models the behavior of the sensor that it
can only report plant outputs in every T seconds. Adding these effects to the overall
simulation of the system provides more realistic performance studies.
Finally, in today’s engineering applications, a common approach is to use
communication networks instead of dedicated digital communication lines as
shown in Fig. 3.18a, b. In fact, the benefit of this networked structure is being able
to integrate as many components together with the capability of increased resources
and easy maintenance as shown in Fig. 3.18c (Cakmakci and Ulsoy 2009).
However, with the introduction of networks, the communication among system
components can experience delays (or even loss of contact) as reported and studied
by many researchers (Lian et al. 2002; Walsh et al. 2002). To remedy this effect, the
overall system can be simulated using worst communication delays possible to
measure the performance using step size-based delay elements and the controllers
are calibrated accordingly.

3.4.2 Robotics and Cyber-Physical Systems

One of the important use of simulations that predict system performance after the
product design phase is to employ them as observers and/or monitoring threads in
actual systems running in parallel and making predictions/modifications to improve
system performance.
A good example of this type of utilization is the friction observers in robotic
locomotion devices such as the one developed in Ristevski and Cakmakci (2015) as
shown in Fig. 3.19. Many non-wheeled robotic systems observe the friction force
during translation. Inside the controller, a simulation of the whole system based on
the dynamic force balance is run to predict the effective friction force called the
friction observer. The friction predictions from this observer is used to update level
68 M. Cakmakci et al.

(a) Vibration based translational system

(b) Translation Controller using a Friction Observer

Fig. 3.19 Translating mechatronic device and friction observer

of the actuator force given to the system as an offset in parallel with its feedback
controller so that the response performance can be improved almost 25%.
Another application of after-design simulation work in engineering systems is
the pre-analysis and optimization of inputs embedded in computers of the manu-
facturing systems. Manufacturing of small parts can be costly and cumbersome
since it often requires trial and error of adjustment of the machine settings.
However, a remedy to this can be found by use of virtual iterative learning as
reported in (Türeyen et al. 2016). A simulation of the additive-manufacturing
system can be developed and used in parallel with a learning algorithm on the
dimensional error of the final part before the real production is actually ran as
shown in Fig. 3.20a. Researchers report using this method can improve the
dimensional accuracy of a representative part up to 75% (Fig. 3.20b).
It is finally noted that similar to MDO based efforts for designing multidisci-
plinary systems such as automotive and aerospace products, there has been a con-
tinuous effort to design controlled mechanical systems using co-design strategies
(Patil et al. 2010). The ultimate goal in these studies is to develop design frameworks
that allow to reach system optimal designs from both the control and the mechanical
design perspectives. Toward that goal, one such recent study is performed in
(Kamadan 2016) where co-design strategies are proposed for robotic systems.
3 Simulation-Based Engineering 69

First Run

5th Run Cure Simulations

(a) Algorithm and Simulations

(b) Production Example

Fig. 3.20 Improving manufacturing quality with preproduction iterations

Robotic systems designed using domain-specific conventional approaches result in

underperforming systems, i.e., are not system-optimal. This work introduces for the
first time a unified framework of system-optimal designs of nonlinear controlled
robotic systems driven by compliant actuators spanning a range of designs.

3.5 Conclusions and Outlook of the Topic

The ultimate objective within nearly all engineering projects is to reach a functional
design without violating any of the performance, cost, time, and safety constraints
while optimizing the design with respect to one of these metrics. Generally, in the
beginning of each project, wish list like high-level requirements for the msyste
70 M. Cakmakci et al.

performance are specified. Then, high-level requirements are cascaded down to the
lower levels of the system allowing systematic design steps to be applied to these
well-defined engineering design problems. The resulting problems are concrete
problem constructs that contain quantifiable performance and constraint metrics. In
time, two primary approaches emerged for the solution of complex engineering
design projects. With the early approach also known as the “Waterfall Design
Process”, the subproblems can be tackled and solved sequentially. In recent years,
as an extension of “Waterfall Design Process”, a new approach has emerged called
the “V-model” where scalable and varying fidelity simulations plays an important
role before the actual prototype of the system can be build.
A good mathematical model is at the heart of each powerful engineering sim-
ulation being a key component in the design process. These models can be obtained
by using physics-based methods, empirical collections and analysis or a combi-
nation of these two for balanced fidelity and complexity. Another important aspect
of developing simulations is its resolution, or in other words its building blocks. In
the simulation of the continuum, systems can be built from their smallest elements
using the most fundamental forms of the governing equations. Sometimes, a
lumped parameter-based simulation of machinery approach can be taken to simplify
the simulations and the forthcoming engineering work such as in the case of the
model based control system design.
The introduction of efficient and powerful platforms enabled researchers to
solve/simulate the constitutive laws of continuum in mechanics in combination with
the laws of conservation of mass, energy, and momentum. The same is valid for
other fields including fluid mechanics and EM. Some of the most popular methods
used for this purpose are the Finite Element Method (FEM), Finite Volume
Methods (FVM), Finite-Difference Methods (FDM), and Boundary Element
Methods (BEM). These methods are applied to the simulation of matter in all forms,
i.e., solids, liquid, and gas, based on a major assumption of continuum media, thus
Computational Mechanics of Continua. Namely, continuum describes the non-
separability of the considered domain and validity of continuity between any points
in the domain so that differentiation is possible. Therefore, continuity between
elements in any continuum based numerical technique is maintained as well.
Many realistic problems present themselves as very complex problems due to
their multi-physics and multi-scale nature. More specifically, scientists and engi-
neers from various fields have been working on the combination of different
numerical techniques with the goal of addressing these complex and elaborate
physical processes, such as the transition from continuum to discontinuum (e.g.,
fracture processes) or the interaction of multi-phases of matter (e.g., hydrofracture
processes) at micro–macro scales. As a result, new classes of numerical methods
called hybrid or multi-physics and multi-scale methods evolved.
In today’s world, developing multidisciplinary systems such as for instance
cyber-physical systems that consist of both mechanical and electrical components
constitute a key part of the engineering projects. These types of systems contain
software algorithms, digital sampling and power electronics, mechanical compo-
nents as well as communication networks that have to be developed concurrently to
3 Simulation-Based Engineering 71

work coherently. This critical need for coherence has rapidly increased the
importance of developing multi-domain simulations and engineers capable of
supporting multidisciplinary analysis and design methodologies.

Review Questions
1. What are the primary phases of the engineering design cycle and how is sim-
ulation work used in each of them?
2. Name some of the numerical methods that exist to solve continuum problems.
3. Which type of problems can be classified as multi-physics and multi-scale
continuum problems?
4. What challenges exist today when solving continuum problems?
5. What do all numerical methods suited to solve continuum problems have in
common?
6. What are the effects of not-considering supporting hardware (networks, digital
computers, etc.) to simulation performance in engineering systems?
7. What are possible uses for the simulations developed for the design cycle after
product release to improve performance?

References

Ashby, M. F. (1996). Modelling of materials problems. Journal of Computer-Aided Materials

Design, 3(1–3), 95–99.
Automotive Suspension - MATLAB & Simulink Example. (n.d.). Retrieved February 15,
2017, from [Link]
html?requestedDomain=[Link].
Bathe, K.-J. (2006). Finite element procedures. Klaus-Jurgen Bathe.
Bensoussan, A., Lions, J.-L., & Papanicolaou, G. (1978). Asymptotic analysis for periodic
structures (Vol. 5). North-Holland Publishing Company Amsterdam.
Cakmakci, M., & Ulsoy, A. G. (2009). Improving Component-Swapping Modularity Using
Bidirectional Communication in Networked Control Systems. IEEE/ASME Transactions on
Mechatronics, 14(3), 307–316. [Link]
Cakmakci, M., Li, Y., & Liu, S. (2011). Model-in-the-loop Development for a Fuel Cell Vehicle.
In Proceedings of the American Control Conference (pp. 2462–2467).
Chen, C.-T. (1995). Linear system theory and design. Oxford University Press, Inc.
Chew, W. C. (1995). Waves and Fields in Inhomogeneous Media (Vol. 522). New York: IEEE
Press.
Chew, W. C., Michielssen, E., Song, J. M., & Jin, J.-M. (2001). Fast and efficient algorithms in
computational electromagnetics. Artech House, Inc.
Dhatt, G., Lefrançois, E., & Touzot, G. (2012). Finite element method. John Wiley & Sons.
Dokuyucu, H. I. H. I., & Cakmakci, M. (2016). Concurrent Design of Energy Management and
Vehicle Traction Supervisory Control Algorithms for Parallel Hybrid Electric Vehicles. IEEE
Transactions on Vehicular Technology, 65(2), 555–565. [Link]
2405347.
El-Kahlout, Y., & Kiziltas, G. (2011). Inverse synthesis of electromagnetic materials using
homogenization based topology optimization. Progress In Electromagnetics Research, 115,
343–380. doi:10.2528/PIER10081603 [Link]
72 M. Cakmakci et al.

Franklin, G. F., Powell, J. D., & Emami-Naeini, A. (2009). Feedback Control of Dynamic Systems
(6th ed.). Prentice Hall. Retrieved from [Link]
Fung, Y.-C. (1965). Foundations of solid mechanics. Prentice Hall.
Gurel, L., & Ergul, O. (2007). Fast and accurate solutions of extremely large integral-equation
problems discretised with tens of millions of unknowns. Electronics Letters, 43(9), 499–500.
Harrington, R. F. (2001). Time-Harmonic Electromagnetic Fields. New York. IEEE Press.
Jategaonkar, R. V, Fischenberg, D., & Gruenhagen, W. (2004). Aerodynamic modeling and
system identification from flight data-recent applications at dlr. Journal of Aircraft, 41(4),
681–691.
Kamadan, A. (2016). Development of Co-design frameworks for optimal variable compliant
actuation. Sabanci University.
Karnopp, D. C., Margolis, D. L., & Rosenberg, R. C. (2000). System Dynamics- Modeling and
Simulation of Dynamic Systems (Third Edit). Wiley-Interscience.
Kazachkov, I. V, & Kalion, V. A. (2002). Numerical Continuum Mechanics. Lecture notes. KTH.
Khurmi, R. S., & Gupta, J. K. (1976). Theory of machines. Eurasia.
Kiziltas, G., Psychoudakis, D., Volakis, J. L., & Kikuchi, N. (2003). Topology design optimization
of dielectric substrates for bandwidth improvement of a patch antenna. IEEE Transactions on
Antennas and Propagation, 51(10), 2732–2743.
Lian, F., Moyne, J., & Tilbury, D. (2002). Network Design Consideration for Distributed Control
Systems. IEEE Transactions on Control Systems Technology, 10(2), 297–307.
Malvern, L. E. (1969). Introduction to the Mechanics of a Continuous Medium.
Martins, J. R. R. A., & Lambe, A. B. (2013). Multidisciplinary Design Optimization: A Survey of
Architectures. AIAA Journal, 51(9), 2049–2075. [Link]
Milton, G. W. (2002). The theory of composites (Cambridge monographs on applied and
computational mathematics).
Muntean, A. (2015). Continuum Modeling: An Approach Through Pratical Examples. Springer.
Ogata, K. (1990). Modern Control Engineering. Prentice Hall.
Ogata, K. (1995). Discrete-Time Control Systems (2nd ed.). Prentice Hall. Retrieved from
[Link]
Patil, R., Filipi, Z., & Fathy, H. (2010). Computationally Efficient Combined Design and Control
Optimization using a Coupling Measure. IFAC Proceedings Volumes, 43(18), 144–151.
[Link]
Rajamani, R., Choi, S. B., Law, B. K., Hedrick, J. K., Prohaska, R., & Kretz, P. (2000). Design
and Experimental Implementation of Longitudinal Control for a Platoon of Automated
Vehicles. Journal of Dynamic Systems, Measurement, and Control, 122(3), 470–476.
[Link]
Rajamani, R., & Hedrick, J. K. (1995). Adaptive observers for active automotive suspensions:
theory and experiment. IEEE Transactions on Control Systems Technology, 3(1), 86–93.
[Link]
Ristevski, S., & Cakmakci, M. (2015). Mechanical design and position control of a modular
mechatronic device (MechaCell). In 2015 IEEE International Conference on Advanced
Intelligent Mechatronics (AIM) (Vol. 2015–August, pp. 725–730). IEEE. [Link]
1109/AIM.2015.7222623.
Schetz, J. A., & Fuhs, A. E. (2013). Fundamentals of fluid mechanics. John Wiley & Sons.
Sokolowski, J. A., & Banks, C. M. (2010). Modeling and simulation fundamentals: theoretical
underpinnings and practical domains. John Wiley & Sons.
Topçu, O., Durak, U., Oğuztüzün, H., & Yilmaz, L. (2016). Distributed Simulation. Cham:
Springer International Publishing. [Link]
Tureyen, E. B., Karpat, Y., & Cakmakci, M. (2016). Development of an iterative learning
controller for polymer based micro-stereolithography prototyping systems. In 2016 American
Control Conference (ACC) (pp. 852–857). IEEE. [Link]
Ulsoy, A. G., Peng, H., & Çakmakci, M. (2012). Automotive Control Systems [Hardcover].
Cambridge University Press. Retrieved from [Link]
Systems-Galip-Ulsoy/dp/110701011X.
3 Simulation-Based Engineering 73

Volakis, J. L., Mumcu, G., Sertel, K., Chen, C.-C., Lee, M., Kramer, B., … Kiziltas, G. (2006).
Antenna miniaturization using magnetic-photonic and degenerate band-edge crystals. IEEE
Antennas and Propagation Magazine, 48(5).
Walsh, G. C., Ye, H., & Bushnell, L. (2002). Stability Analysis of Networked Control Systems.
IEEE Transactions on Control Systems Technology, 10(3), 438–446.
Wang, G. G., & Shan, S. (2007). Review of metamodeling techniques in support of engineering
design optimization. Journal of Mechanical Design, 129(4), 370–380.
Widder, D. V. (1976). The heat equation (Vol. 67). Academic Press.
World-Record Algorithm from Jülich Calculates Over Three Trillion Particles - Research in
Germany. (2011). Retrieved February 15, 2017, from [Link]
research-landscape/news/news-archive/2011/08/2011-08-02-world-record-algorithm-from-
[Link].
Zienkiewicz, O. C., & Taylor, R. L. (2005). The finite element method for solid and structural
mechanics. Butterworth-heinemann.

Author Biographies

Melih Cakmakci is Assistant Professor of Mechanical Engineering at Bilkent University in

Ankara, Turkey. He received his B.S. in Mechanical Engineering from M.E.T.U Ankara in 1997.
He received his M.S. and Ph.D. in Mechanical Engineering from University of Michigan in 1999
and 2009 respectively. His research areas include modeling, analysis and control of dynamic
systems. Prior to joining Bilkent University, he was a Senior Engineer at the Ford Scientific
Research Center at Dearborn, Michigan, USA. He is a member of ASME, IEEE, and SAE.

Gullu Kiziltas Sendur is Associate Professor with the Mechatronics Engineering Program,
Sabanci University. She received the [Link]. and [Link]. in Mechanical Engineering from the Middle
East Technical University, Ankara, in 1995 and 1998, respectively, and the Ph.D. in mechanical
engineering from the University of Michigan, Ann Arbor, in 2003. She was a Postdoctoral
Researcher with both the Electro-Science Laboratory, Ohio State University and the University of
Michigan. She also coordinated with the Ceramic Research Group, Material Science and
Engineering Department, University of Michigan, on the advanced fabrication of dielectric
composites. Her research efforts center on the design, analysis, and fabrication of complex
engineering systems, such as miniaturized electromagnetic, electromechanical, and biomedical
devices and multidisciplinary design optimization techniques. She has published a book chapter,
numerous journal and conference articles and was honored with Career Award by TUBITAK
(2007) and Young Scientist Award by the Science Academy of Turkey (2008). She is a member of
the IEEE, ASME, and ACerS.

Umut Durak is a research scientist in the Institute of Flight Systems at the German Aerospace
Center. He is also Adjunct Lecturer in the Department of Informatics at Clausthal University of
Technology. His research interests include model-based simulation engineering, simulation-based
systems engineering and ontologies in simulation. He received his BS, MS and Ph.D. degrees in
Mechanical Engineering from Middle East Technical University (METU) in Turkey. He is a
member of the Society for Computer Simulation International (SCS) and Arbeitsgemeinschaft
Simulation (ASIM) and a senior member of the American Institute of Aeronautics and
Astronautics (AIAA). And, he is a part of AIAA Modeling and Simulation Technical Committee
and editorial teams of the Simulation: Transactions of the Society for Modeling and Simulation
International; and the International Journal of Modeling, Simulation, and Scientific Computing.