Document [original]

This version is available at https://doi.org/10.14279/depositonce-10661

right to use is granted. This document is intended solely for

personal, non-commercial use.

Marinkovic, D., Zehn, M. & Rama, G. (2018). Towards real-time simulation of deformable structures by

means of co-rotational finite element formulation. Meccanica, 53(11–12), 3123–3136.

https://doi.org/10.1007/s11012-018-0868-5

This is a post-peer-review, pre-copyedit version of an article published in Meccanica. The final

authenticated version is available online at: http://dx.doi.org/10.1007/s11012-018-0868-5.

Marinkovic, D.; Zehn, M.; Rama, G.

Towards real-time simulation of

deformable structures by means of co-

rotational finite element formulation

Accepted manuscript (Postprint)Journal article |

Meccanica manuscript No.

(will be inserted by the editor)

Towards real-time simulation of deformable structures by

means of co-rotational finite element formulation

Dragan Marinkovic ·Manfred Zehn ·Gil Rama

Received: date / Accepted: date

Abstract Whereas typical Finite Element (FE) com-

putations are performed off-line, many virtual-reality

(VR) applications put a demand for interactive simula-

tions involving deformable objects. Interactive simula-

tion implies real-time or nearly real-time computation

and graphical representation of modeled deformable ob-

jects. The growing computational power of modern con-

ventional hardware calls for FEM developments in this

direction. Depending on specific VR applications, the

developments need to account for different aspects of

physical behavior, with geometrically nonlinear defor-

mations emerging as one of most important and fre-

quent. This paper proposes a simplified co-rotational

FE formulation that considers the overall finite ele-

ment motion as a superposition of rigid-body rotation

and deformation described by a linear model with re-

spect to the co-rotated reference frame. By neglecting

the stress stiffening effect and the dependence of the

element stiffness matrix on the deformational displace-

ments, the formulation aims at meeting the objectives

of highly efficient simulation, under certain conditions

even real-time simulation, and with acceptable devia-

tion in accuracy compared to the rigorous nonlinear FE

formulation. Computations with both solid and shell

elements are addressed. A set of examples is provided

to illustrate and discuss the aspects of accuracy and

achievable simulation speed.

D. Marinkovic

Department of Structural Analysis, Berlin Institute of Tech-

nology, Germany

Strasse des 17. Juni 135, 10623 Berlin

Tel.: +49-30-31421483

E-mail: dragan.marinkovic@tu-berlin

M. Zehn

E-mail: manfred.zehn@tu-berlin ·G. Rama

E-mail: gil.rama@tu-berlin

Keywords Real-time simulation ·Co-rotational

FEM ·Geometric nonlinearity; Solid ·Shell

1 Introduction

In the last several decades, the tools for computer aided

engineering (CAE) have offered invaluable assistance to

engineers. Their continuous development aims at pro-

viding the maximum in terms of performance and reli-

ability in a variety of engineering tasks. An important

ingredient of those tools is the simulation of physical

processes and phenomena. The conventional approach

implies off-line simulations followed by the assessment

and analysis of the obtained results in postprocessors

in order to gain an insight into the physical quanti-

ties of interest that characterize the processes. However,

the strong development of advanced visualization sys-

tems, followed by further development of already exist-

ing software components and necessary hardware, have

enabled the integration of a novel concept into many ar-

eas of engineering and other fields of application. The

concept is based on the Virtual Reality (VR) or Aug-

mented Reality (AR) technology and real–time simula-

tions, thus offering the possibility of manipulating and

analyzing a 3D virtual world.

The finite element method (FEM), as an established nu-

merical method in the field of structural analysis, is typ-

ically addressed to compute the behavior of deformable

structures. Although numerically demanding, with the

increasing computational power of modern hardware it

is also the method of choice for interactive and real-time

simulations, either used directly for the computation or

as a part of particularly developed techniques. The pos-

sibility of FEM-based real-time simulations has been a

subject of interest in a number of fields.

2 Dragan Marinkovic et al.

With the purpose of enhancing structural analysis with

AR technologies, Huang et al. [17] proposed a system

that integrates sensor measurement and real-time FEA

simulation into an AR-based environment and superim-

poses the FEA results directly onto real-world objects.

In order to speed up the computation, the authors intro-

duce the assumption of linear and quasi-static behav-

ior, while the real-time solver is based on the concept of

pre-computing the inverse of the stiffness matrix. Sim-

ilarly, Fiorentino et al. [15] presented an AR-based ap-

plication that visualizes the linear FEM results overlaid

over the real model for interactive teaching of dynamic

stress/strain distribution in engineering education. Cer-

racchio et al. [6] implemented the linear inverse FEM

for real-time reconstruction of the deformed structural

shape using in situ strain measuremenets. Adopting lin-

ear behavior in the aforementioned applications allowed

significant numerical savings and therewith acceleration

of computational processes. If, however, the considera-

tion of nonlinear behavior is a necessity for reasonable

simulation accuracy, the numerical effort and the com-

plexity of simulation increase notably.

A number of authors have addressed the difficulty of

performing real-time FEM computations of nonlinear

structural deformations by suggesting approaches based

on neural-networks to obtain highly efficient, real-time

solutions. Such approach consists in using the nonlinear

FEM computations to provide a set of results for the

training phase and can thus be classified as model re-

duction technique. It was applied to a number of prob-

lems such as real-time simulation of impact [16], geo-

metrically nonlinear deformation of a cantilever beam

[38] and a shell structure [8], etc. Dulong et al. [12] pro-

posed a similar approach for real-time interaction be-

tween a designer and a virtual prototype as a support

to design optimization. They also used a set of nonlin-

ear FEM computations in what they refer to as training

phase, while the actual deformation for a specific load

case is determined by an interpolation method proposed

by the authors. Furthermore, a real-time control system

can also benefit from reduced but sufficiently accurate

models of the system being controlled to improve the

predictor part of the controller. Kalkkuhl et al. [20]

applied the nonlinear FEM-based neural-network ap-

proach to a longitudinal vehicle dynamics control. Gen-

erally speaking, the approach based on a pre-computed

set of deformed configurations is promising but the va-

lidity of the resulting model is strongly dependent on

the deformation set density and range of deformation

covered in the training phase.

In the field of multi-body system (MBS) dynamics with

flexible bodies, real-time simulation is not a strict re-

quirement, but high simulation efficiency belongs to the

priorities. Though apparent savings are made by con-

sidering parts mainly as rigid-bodies, accounting for de-

formational behavior of some parts becomes crucial in

certain simulations. In such a case, model reduction is

one of the basic ideas on how to cope with the com-

putational burden. On the contrary to the above men-

tioned technique of model reduction (training of neural-

networks), a kind of direct FE-model reduction is intro-

duced here based on modal coordinates. It implies that

orthogonal mode shapes, calculated in a step prior to

simulation, represent the degrees of freedom, in terms of

which the elastic behavior of the body is described. The

solution used by commercial MBS software packages is

the Component Mode Synthesis (CMS) technique, par-

ticularly the Craig-Bampton method [9]. The flexible

behavior remains linear with respect to the floating ref-

erence frame attached to the body. Some modal-based

solutions are proposed to account for moderate geomet-

ric nonlinearities by means of the geometric stiffness

matrix [41],[26],[27] or modal warping [27],[7],[29]. It is

also possible to divide a model into a linear part that is

further reduced and a part that considers the geometri-

cally nonlinear effects in order to efficiently account for

local nonlinearities [4]. The success of these solutions is

strongly dependant on the character and range of de-

formation and modes provided, so they need to be used

with special care.

Many solutions for interactive VR-simulations have been

developed and most advances have been motivated by

the demand for providing real-time deformation of hu-

man tissue in surgical simulations. Pioneer works in

the 90s were based on linear FE models. Bro-Nielsen

and Cotin [35] used simple tetrahedral elements in the

framework of linear elasticity. In order to further reduce

the numerical effort during simulation, the stiffness ma-

trix was condensed so that only the surface nodes were

involved and the so-obtained stiffness matrix was ad-

ditionally inverted prior to simulation. All these steps

done to produce high numerical efficiency reflect the

limitations of available hardware at that time. Model

reduction techniques were also exploited in this field,

combined with nonlinear visco-elastic material mod-

els [11] or with the extended finite element method

(XFEM) [36]. Certain developments rely on use of mod-

ern hardware tools to provide the necessary computa-

tional power for the demanding real-time simulations

[18]. Cueto and Chinesta [10] gave a survey of devel-

oped solutions for real-time simulation in surgery in-

cluding techniques that involve supercomputing facili-

ties, parallel implementations on GPUs (graphics pro-

cessor units), model order reduction, etc.

Towards real-time simulation of deformable structures by means of co-rotational finite element formulation 3

2 Setting objectives and the choice of method

Strictly speaking, real-time simulation could be under-

stood as a simulation, in which the computer model

runs at the same rate, or even faster, than the actual

physical system. The standard DIN ISO/IEC 2382 [40]

defines real-time simulation in a somewhat looser man-

ner, by essentially stating that it refers to the operation

of a computing system, whereby the programs process

data in such a manner that the processing results are

available within a predetermined period of time. Ac-

cording to the standard, the processing time is not nec-

essarily the same or faster than the real-time. Indeed,

in certain fields of application, a simulation would still

do what it is meant for, even if it performed somewhat

slower than the real time, but the resulting simulation

could be played at an interactive frame rates.

One can easily identify the means that can be addressed

to achieve the objective of real time simulation:

–Powerful hardware components,

–Optimization and parallelization of computer pro-

grams,

–Appropriate formalisms and algorithms for the com-

putation of physical processes.

In the recent years, hardware components have seen

a rapid pace of development offering ever increasing

computational power. This refers to both the central

processing units (CPU) and graphics processing units

(GPU). Particularly modern GPUs can be used as a

modified form of stream processors providing a massive

computational power [18], [37].

Program optimization with respect to the requirements

of real time simulation implies modifications in the pro-

gram so that it executes more rapidly. Modern compil-

ers offer significant assistance in this aspect. As modern

processors contain multiple cores, the possibility of par-

allelization is another important aspect. This may also

affect the choices made in the third group, because dif-

ferent solver types possess different parallelization po-

tential. By the choice of adapted formalisms to describe

and compute deformational behavior of structures, one

can significantly affect the numerical efficiency and pos-

sibly make a trade-off between the efficiency and accu-

racy. Certain solutions even allow to perform this trade-

off dynamically, i.e. during the simulation, based on im-

mediate requirements.

This paper is focused on the third aspect. It aims at

a FEM-based solution that keeps full fidelity FE mod-

els (no model reduction is performed), covers geomet-

ric nonlinearities to a large extent, performs fast even

on conventional hardware tools and offers acceptable

accuracy, the definition of which depends on specific

application. Keeping the focus on developments with

similar objectives, available literature offers several ap-

proaches. The simplest one would be a model based

on mass-spring system, which substitutes a continuum

with a collection of point masses connected by a net-

work of massless springs, which are reminiscent of early

days of discretization methods. The main advantages

are the simplicity and ability to cope with geometric

nonlinearities with a relatively small numerical effort

and hence the interest in the approach for the purpose

of real-time simulation in various simulators [30],[22].

But it also brings serious drawbacks related to achiev-

able accuracy and ambiguity of model building [32].

Due to rapid hardware development the attention was

turned to more demanding FE models. The rigorous ge-

ometrically nonlinear FEM is computationally very de-

manding and, hence, carefully chosen simplifications of

the rigorous geometrically nonlinear FEM are required

in order to meet the aforementioned objectives. The

potential of the co-rotational approach was recognized

in a number of works, particularly in the field of com-

puter graphics. Capell et al. [5] proposed division of

objects into sub-domains, while their local rigid-body

rotations were considered by means of local coordinate

frames. A similar but extended idea was proposed by

Mueller et al. [31] who implemented local coordinate

frames at nodes. The idea was modified by Etzmuss et

al. [13] by using local coordinate frames attached to tri-

angular elements in order to model cloth behavior. The

authors of the present work consider the possibility of

using a simplified co-rotational FEM formulation with

element-based rotation in combination with shell and

solid elements to meet the above set objectives.

3 A simplified co-rotational FEM formulation

The total and updated Lagrangian formulations [2] are

known as classical FE formulations for geometrically

nonlinear analysis used in most commercially available

FE software packages. The only difference between the

two formulations lies in the choice of different reference

configurations and, provided the appropriate constitu-

tive tensors are employed, they yield identical results

[2]. They provide the accuracy necessary for engineer-

ing tasks, but are also numerically quite demanding and

may also exhibit convergence issues.

A high-quality overview of different co-rotational for-

mulations including a detailed analysis of their prop-

erties is given by Felippa and Haugen [14]. Having in

mind the set objectives, the authors of the present work

propose a rather simplified co-rotational FE formula-

tion that combines the advantages of the linear and

geometrically nonlinear FE analysis. The idea of the

4 Dragan Marinkovic et al.

approach is to cover geometric nonlinearities to a signif-

icant extent through the consideration of the element-

based rigid-body rotation and to neglect further effects

such as the geometric stiffness and the dependence of

the element stiffness matrix on the deformational dis-

placements. The idea can be seen as a refinement of

the approach used in MBS programs. Within MBS the

overall flexible body motion is decomposed into a rigid-

body motion and (typically small, linear) deformable

motion. These simplifications mean that the formula-

tion is not suitable for certain types of analyses, such as

buckling analysis, which is primarily due to the lack of

the geometric stiffness matrix. The same idea is used in

this work, but it is applied on the element level. Hence,

in the present approach the elastic behavior of each ele-

ment is considered to be linear with respect to the local

element frame that is attached to the element and fol-

lows its rigid-body motion. In what follows, the compu-

tation of internal forces and tangential stiffness matrix

together with some details of handling translational and

rotational degrees of freedom within the framework of

the applied co-rotational formulation will be addressed.

Given the rotation matrix at time t,tRe, which de-

scribes the rigid-body rotation of an element between

its initial configuration, 0xe, and its current configura-

tion, txe, the rotation-free translations between the two

configurations, read (see Fig. 1):

0uR=tRe

Ttxe−0xe(1)

where the left superscript denotes the time at which

the quantity of interest is given, while the left sub-

script (where applicable) refers to the element orien-

tation with respect to which the quantity is given. The

left subscript is omitted with the rotation matrix as it

is always given with respect to the initial element con-

figuration (at time t= 0). As already mentioned, the

Fig. 1 Co-rotational concept rotation-free translations in

the original element orientation

element behavior is considered to be linear with respect

to the local reference frame and, hence, the element in-

ternal forces with respect to the initial configuration

can be simply computed by means of a pre-computed

element stiffness matrix, 0Ke:

0Fe=0KetRe

Ttxe−0xe(2)

To proceed with the solution, the internal forces are

rotated by means of tReto the current element orien-

tation:

tFe=tRe

0Fe=tRe0Ke

tRe

Ttxe−tRe0Ke0xe(3)

It should be noticed that the term 0Ke0xecan be com-

puted in a pre-step, i.e. prior to simulation and further

used during the simulation. Also, Eq. (3) reveals the

updated element stiffness matrix, which is simply given

by rotating the initial, linear element stiffness matrix,

0Ke, through tRe:

tKe=tRe0Ke

tRe

T(4)

As solids have only translational degrees of freedom and

nodal forces as loads, Eqs. (1-4) give the essence of the

co-rotational FEM formulation for this type of finite el-

ements. Typically, shell elements also employ rotations

as degrees of freedom and the procedure needs to be

adequately extended. The rotations and translations do

not share the same properties and the update of rota-

tional degrees of freedom is more demanding. The incre-

mental nodal rotations computed in each time-step are

used to update the shell normals at each node starting

from normals in the previously determined configura-

tion. This is done by computing the incremental rota-

tion matrix of a shell normal as [1]:

t−∆tQi=

I+sin t−∆tγi

t−∆tγi

t−∆tSi+1

2 sin t−∆tγi/2

(t−∆tγi/2) !2

t−∆tSi

2(5)

where

t−∆tγi=qt−∆t∆θ2

i1+t−∆t∆θ2

i2+t−∆t∆θ2

i3(6)

and

t−∆tSi=



0−t−∆t∆θi3t−∆t∆θi2

t−∆t∆θi30−t−∆t∆θi1

−t−∆t∆θ2

t−∆t∆θi10

(7)

with t−∆t∆θi1,t−∆t∆θi2and t−∆t∆θi1denoting the 3

incremental global nodal rotations between the config-

urations at times t−∆t and t, while the index irefers

to node i. The rotation matrix t−∆tQiis further used

to update the orientation of the shell normal at node i,

ni:

tni=t−∆tQt−∆tni(8)

Towards real-time simulation of deformable structures by means of co-rotational finite element formulation 5

Rotations are further handled in a similar manner as

translational degrees of freedom. The updated node nor-

mal is rotated backwards to the initial element config-

uration by means of tRe

tnRi =tRe

Ttni(9)

and further compared to the original node normal in

order to determine the deformational nodal rotations,

free of rigid-body rotation, and with respect to the ini-

tial configuration, Fig. 2. Now, the same approach al-

ready elaborated in Eqs. (2) and (3) is applied. The

deformational translations and rotations are used with

the linear stiffness matrix to obtain the internal forces

and moments with respect to the initial configuration,

which are finally rotated to the current element config-

uration.

Fig. 2 Co-rotational concept handling of shell normals

Obviously, extraction of the element rotation matrix

from the overall motion is an important aspect in the

co-rotational formulation. Generally speaking, each in-

cremental volume of the structure may exhibit different

rigid-body rotation. But in the framework of the present

co-rotational FEM, a single rotation matrix is used per

element which implies averaging of the rigid-body ro-

tation over the element domain. The rotation matrix

at a material point can generally be obtained by po-

lar decomposition of the deformation gradient matrix

at that point. For relatively simple elements, such as

the linear tetrahedron or linear triangular element, the

deformation gradient matrix is constant over the whole

element domain which fits nicely into the formulation.

With more complex elements (e.g. quadratic elements),

the polar decomposition of the deformation gradient

matrix at the element centroid is used here.

Hence, the formulation uses the element linear stiffness

matrix, which is computed only once, in a pre-step. In

further computation, the element rotation matrix be-

tween the initial and current configurations is deter-

mined in order to update the element stiffness matrix

and compute the internal forces and moments (where

applicable). The global stiffness matrix is re-assembled

using the rotated element stiffness matrices and the vec-

tor of internal forces is updated, so that either static or

dynamic nonlinear analysis may proceed. It should be

noticed that the formulation neglects the influence of

the change in the element shape and initial stress state

onto the element stiffness matrix.

4 Aspect of accuracy static shell examples

As emphasized above, the presented formulation imple-

ments certain simplifications, i.e. neglects certain ge-

ometrically nonlinear effects. An important question

rises: how is the accuracy of obtained results affected

by those simplifications?

The achievable accuracy in terms of numbers by means

of the presented co-rotational FEM in combination with

the linear tetrahedral element has been considered by

Marinkovic et al. [28]. Additionally, in the same refer-

ence, a technique was proposed to improve the accu-

racy in cases involving moderate strains, but it implies

a larger numerical effort because the co-rotational for-

mulation is extended by a secant approach with an up-

date of the element stiffness matrix. Another approach

to improve the accuracy with some additional numer-

ical effort is based on the so-called projector matrix

[14],[34].

Thin-walled structures are known for their susceptibil-

ity to large rotations, whereby the strains remain small.

This is why they represent a good candidate to consider

this aspect. A couple of illustrative examples are chosen

for this consideration. Those examples have also been a

subject of interest of other authors (e.g. [23],[21],[19]).

It is not intended here to study this aspect exhaustively

but rather to get a general impression on what can be

expected. In typical applications involving interactive

simulations, displacements are the result of primary in-

terest.

Two shell elements were implemented with the pre-

sented co-rotational formulation. The full biquadratic

nine-node shell element, denoted here by S9 (Fig. 3,

left), belongs to the family of degenerated shell ele-

ments. It was originally developed as a piezoelectric

shell element by Marinkovic et al. [24] and tested in the

commercially available FEM software package ABAQUS

[33]. The linear triangular shell element, denoted by S3

(Fig. 3, right), is the mechanical part of the electro-

mechanical shell element presented and tested by Marinkovic

and Rama [25] and Rama [39]. A few remarks on the

shell elements would be worthwhile at this point. Both

elements implement the Mindlin-Reissner kinematics

and, hence, include transverse shear. Due to the flat

shape of the S3 element, a FE discretization of a curved

6 Dragan Marinkovic et al.

Fig. 3 Full biquadratic quadrilateral and linear triangular

shell elements

geometry is represented as a set of facets thus demand-

ing a finer mesh. However, the element itself is compu-

tationally very efficient and, additionally, the finer dis-

cretization also talks in favor of the applied co-rotational

formulation (rigid-body rotation accounted for element-

wise). The full biquadratic shape functions of the S9

element facilitate covering of a relatively large range

of curvatures. Although this implies that reasonably

rough meshes can be used to discretize a complex geom-

etry with the S9 element (compared to linear elements),

this advantage may easily be lost with substantial geo-

metrically nonlinear deformations. The results obtained

with the implemented shell elements are compared with

those obtained using ABAQUS 3-node linear shell el-

ement (A-S3) and ABAQUS 8-node biquadratic shell

element (A-S8R) and the rigorous geometrically non-

linear (updated Lagrangian) FEM formulation.

4.1 Curled cantilever beam

A clamped beam exposed to a bending moment at the

free tip is considered in the first example. If the bending

moment is slowly increased, the beam will bend until

it gets curled into a complete circle. The required over-

all bending moment for this deformation is analytically

determined as M=πEbh3/6l[3], where adenotes the

length, bthe width, hthe thickness, and Eis the Youngs

modulus.

The geometry is chosen here so that l= 10 m , b= 0.5 m

and h= 0.01 m, see Fig. 4. As for material properties,

the Youngs modulus is taken to be E= 2 ·1011 N/m2

and in order to use the analytical solution for the beam,

the Poissons coefficient is taken as equal to 0. Hence,

in this specific case the edge distributed bending mo-

ment is obtained as Ml= (π/3)104N. The mesh of

10×1 elements was used with the quadratic quadrilat-

erals, while the discretization with the linear triangu-

lar elements was finer in order to capture the bend-

ing adequately, so that 20 element-stripes were used

along the length with each ’stripe containing 4 elements

across the width, hence 80 elements altogether. These

meshes give converged results with the implemented

S9 and S3 element. Figure 5 shows the initial and de-

formed geometries computed with the S9 element and

co-rotational formulation, with the intermediate con-

Fig. 4 Geometry of the cantilever beam and loading

figurations depicted upon each 20 The computations

Fig. 5 Intermediate configurations at constant load incre-

ments of 20%

with the S9 and S3 elements were performed using in-

crements of 10%. The analysis in ABAQUS was done

using the same meshes (for respective elements) and it

was set to use an automatic increment size with the

initial increment size of 10%. With the A-S8R element

ABAQUS performed 252 increments, but aborted the

analysis at the load level of 87.28% due to convergence

issues. The computation with A-S3 element was com-

pleted (100% of the load) and was done in 314 incre-

ments. As representative results, the displacements of

the beam tip in the x- and y-directions are chosen and

can be seen in Figs. 6 and 7, respectively. All the results

are in very good agreement, with the given remark that

A-S8R element failed to complete the analysis with the

used mesh.

It should also be emphasized that, for better accu-

racy of the analysis, the size of deformation in this

case would also require to account for the change in

the thickness (the thinning in tension and the thick-

ening in compression) [42], which also gives rise to the

neutral surface shift during the deformation. However,

in this analysis the Poisson coefficient was set to zero.

But even if this was not the case, this effect is neglected

with all the used elements as well as in the analytical

Towards real-time simulation of deformable structures by means of co-rotational finite element formulation 7

Fig. 6 Displacement of the beam tip in the x-direction

Fig. 7 Displacement of the beam tip in the y-direction

solution that yielded the bending moment required to

produce the considered deformation.

4.2 Pinched hemisphere

In the second case, a structure with a doubly curved

geometry is considered. It is a hemispherical shell with

an 18◦-hole at the top, radius r= 10 m and thickness

h= 0.04 m, Fig. 8. The Young’s modulus of the mate-

rial is taken to be E= 6.825·1010 N/m2and the Poisson

ratio ν= 0.3. The structure is exposed to two pairs of

concentrated forces acting on the bottom edge of the

shell. The magnitude of each of the forces is 150 kN.

The pair of forces that is symmetric with respect to the

yz-plane acts so as to stretch the hemisphere along the

x-axis, while the pair symmetric with respect to the xz-

plane compresses it along the y-axis.

The double symmetry allows modeling of only a quar-

ter of the hemisphere with the adequate kinematical

boundary conditions applied (Fig. 9, left). In order to

capture the local rotations appropriately, finer meshes

were needed. After convergence analysis, the results

Fig. 8 Pinched hemisphere geometry and loads

obtained with the FE mesh containing 400 (20 ×20)

quadrilateral elements (Fig. 10, left) and the FE mesh

containing 800 triangular elements (Fig. 9, right) were

taken as representative.

For a better insight into the size of deformation the

Fig. 9 Pinched hemisphere - FEM mesh with quadrilateral

(left) and triangular (right) elements

structure is exposed to, Fig. 10 depicts the original and

deformed configurations at the full load from different

perspectives and without scaling, i.e. the scale factor

equals 1.

The diagrams in Figs. 11 and 12 show the development

of the displacements of points A and B (Fig. 9) along

the x- and y-axes, respectively, with the increasing load.

A good agreement of the results by all four elements can

be noticed. It was important in this case to use a suffi-

ciently fine mesh so that the local rigid-body rotations

8 Dragan Marinkovic et al.

Fig. 10 Initial and deformed hemisphere (scale factor equals

are adequately captured in the co-rotational formula-

tion. Whereas in the first considered case the element

rigid-body rotation was practically around one axis, the

overall displacement field in this case is more complex.

Fig. 11 Displacement of point A in the x-direction

Fig. 12 Displacement of point B in the y-direction

5 Aspect of efficiency

In various applications that involve interactive simu-

lation, the simulation speed is the primary aspect. A

distinctive example of such an application is a surgi-

cal simulator, e.g. for the laparoscopic surgery. Due to

the high level of complexity of this type of surgery, sur-

geons frequently emphasize their wish for high-quality

VR-based simulators for the purpose of training. Obvi-

ously, the real-time requirement has the highest priority

in such an application, while the accuracy requirement

is not defined in terms of numbers and may be vio-

lated as long as the on-screen behavior appears realis-

tic to human perception. Though in many commercially

available devices this is achieved by means of the sim-

plest approach, namely the mass-spring systems, the

presented co-rotational FEM formulation represents a

more sophisticated alternative.

Beside the formalism used for the computation of de-

formable body behavior, there are some further com-

putational parameters that play a very important role

with respect to the defined objective. They are not in

the very focus of this paper, but need to be addressed

briefly for adequate comprehension of the results given

below.

The FE mesh plays an important role in every FEM

analysis. In typical engineering FE analyses it needs

to be set with very special care. However, if the re-

quired minimum accuracy is defined as plausible be-

havior, then the idea of coupled meshes can be applied.

Namely, a relatively fine mesh of surface vertices can

be used to represent the actual object geometry in suf-

ficient detail, while a relatively rough FE mesh can be

used for the computation. The surface vertices are con-

nected into triangles to give a triangulated surface rep-

resentation, see Fig. 13, right. The two meshes are cou-

pled to each other based on the local element coordi-

nates of the surface vertices with respect to the finite el-

ements. A vertex that is inside an element is assigned to

that element. Upon deformation the global vertex coor-

dinates are recovered based on its element local coordi-

nates and global position of the element nodes. Hence,

the FE mesh does not have to fit the actual geome-

try, but can instead only resemble the actual geometry

with more (Fig. 14b) or less precision (Fig. 14c). In this

manner the computational burden can be adjusted to

the available hardware and other simulation parame-

ters and requirements.

Time integration, solver type and time-step size are im-

portant aspects in dynamic FE computations. Whereas

the explicit time-integration schemes are very efficient

in computing a single time-step, provided a lumped

mass and damping matrix is used, the time-step size

Towards real-time simulation of deformable structures by means of co-rotational finite element formulation 9

Fig. 13 Liver model geometry: 2598 vertices and 5192 faces

Fig. 14 a) Geometric liver model with two FE meshes: b)

FE mesh A (846 nodes, 2945 elements), c) FE mesh B (660

nodes, 2640 elements)

is on the other hand severely limited due to the condi-

tional stability. Therefore, the explicit integration schemes

are used for highly nonlinear (no iterations) and fast dy-

namics problems. In an interactive simulation, the dy-

namics that is visible to the human eye is of interest, i.e.

relatively slow dynamics. This fact permits to simply

filter high frequency dynamics out, which is effectively

done by using sufficiently large time-steps. This talks in

favor of an implicit time integration with a time-step

that is sufficiently large not only to compensate for the

larger numerical effort per time-step, but also to provide

better overall numerical efficiency compared to explicit

time-integration. An implicit time-integration scheme

keeps the system of equations coupled [17]:

Mt+∆t ¨u +tCt+∆t ˙u(k)+tKt+∆t∆u(k)

=t+∆t Fext −t+∆t F(k−1)

int (10)

where Kand Care the structural stiffness and damp-

ing matrices, respectively, uare the nodal displace-

ments, dots above udenote the time derivatives, i.e.

velocities (one dot) and accelerations (two dots), ∆de-

notes the increment, while Fext and Fint are the nodal

external and internal forces, respectively. The implicit

time-integration of a nonlinear problem demands itera-

tions and (k) in the right superscript denotes the iter-

ation number. In each iteration, the linearized system

of equations can be resolved by a direct or iterative

solver. An iterative solver, such as the method of pre-

conditioned conjugate gradients (PCG), is an interest-

ing choice due to several reasons. In dynamics, partic-

ularly relatively slow dynamics, the velocities do not

change dramatically within a time-step and, hence, a

good choice of the starting vector is already available,

thus reducing the number of iterations notably. A fur-

ther reason would be the fact that a trade-off between

the numerical efficiency and accuracy can be performed

easily by limiting the number of iterations, and this can

even be done dynamically, during a simulation. Finally,

the PCG solver has a great parallelization potential.

Interactive simulations of three solid test objects/models

are considered below to provide an insight into the as-

pect of efficiency by means of the presented formulation

in combination with the linear tetrahedral element and

the PCG solver. The time step of 0.01s is used. The

simulation is set so that 4 configurations are computed

before the current configuration is shown on the screen,

which implies that 25 frames per second would corre-

spond to the real-time computation (25×4×0.01 = 1).

The test objects, i.e. models used in the interactive sim-

ulations are:

–a torus-like model, FE mesh: 2487 elements, 1956

degrees of freedom (DOFs), see Fig. 15a;

–a spleen model, FE mesh: 2226 elements, 1656 DOFs,

the geometric model: 962 vertices, 1920 faces, see

Fig. 15b;

–the liver model depicted in Figs. 13 and 14, with

the finer FE mesh (Fig. 14b): 2945 elements, 2538

DOFs; the geometric model: 2598 vertices, 5192 faces,

see Fig. 15c.

Obviously, the latter two models are motivated by a

possible application in the field of surgical simulation

and they make use of the coupled mesh technique. The

spleen model uses a uniform FE mesh that only roughly

resembles the geometric model, while for the liver model

the adaptive FE mesh A depicted in Fig. 14b is used.

This analysis represents an extension of a similar anal-

ysis given in [28], in which only uniform meshes to-

gether with older hardware configurations were consid-

ered. Figure 16 gives several screen-shots from the in-

teractive simulations with the aforementioned models

demonstrating the ability of the co-rotational FEM to

cover large displacements and rotations producing plau-

sible behavior. Based on the accuracy demonstrated in

the previous section, plausible behavior is expectable,

but one should also have in mind that the coupled mesh

technique is applied in the later two examples. The

same settings are used for an interactive simulation of

two shell models, both meshed using the S3 element:

10 Dragan Marinkovic et al.

Fig. 15 Models for interactive simulations: a) torus-like

structure, b) spleen and c) liver

Fig. 16 screen-shots from an interactive simulation: a) torus

model b) spleen model; c) liver model

–an annular slit plate, with one side of the slit clamped

and the other one free, meshed by 234 elements, 960

DOFs, see Fig. 17a;

–the hemisphere from Section 4.2 with the edge of

the 18-hole clamped, meshed by 1040 elements, 3432

DOFs, see Fig. 17b.

Fig. 17 Shell models for interactive simulations: a) slit an-

nular plate, b) hemisphere with 18-hole

Proceeding in the similar way as with the solid struc-

tures, Fig. 18 gives several screen-shots from the in-

teractive simulations with the considered shell struc-

tures involving rather large deformations computed by

the presented co-rotational formulation. The results re-

garding the simulation speed are summarized in Tables

1 and 2 for different hardware configurations. Three

hardware configurations have been used, all of which

can be described as conventional personal computers.

Information on the CPU and GPU in each configura-

tion is given as those are the components with the main

impact onto the simulation efficiency. The considered

hardware configurations involve multi-core CPUs, but

only one core is used for the computation, i.e. paral-

Fig. 18 Screen-shots from an interactive simulation: a) slit

annular plate; b) hemisphere with the 18-hole

Table 1 Simulation pace of different hardware configura-

tions with solid models

Hardware configuration:

processor graphic card

Sets

Parameter Torus Spleen Liver

Elements 2487 2226 2945

DOFs 1956 1656 2538

Vertices - 962 2598

Faces - 1920 5192

Intel E8500 (3.16 GHz)

NVidia 740 GT

Ratio 1.29 1.51 1.18

F/s 32 38 30

Intel i7-870 (2.93 GHz)

NVidia 650 GTX

Ratio 1.56 1.82 1.43

F/s 39 45 36

Intel i3-2120 (3.3 GHz)

NVidia 750 GTX

Ratio 1.75 2.05 1.62

F/s 44 51 41

Table 2 Simulation pace of different hardware configura-

tions with shell models

Hardware configuration:

processor graphic card

Sets

Parameter Slit Plate Sphere

Elements 234 1040

DOFs 960 3432

Intel E8500 (3.16 GHz)

NVidia 740 GT

Ratio 3.05 0.92

F/s 76 23

Intel i7-870 (2.93 GHz)

NVidia 650 GTX

Ratio 3.69 1.11

F/s 92 28

Intel i3-2120 (3.3 GHz)

NVidia 750 GTX

Ratio 4.14 1.25

F/s 103 31

lelization is not used. The main results in the tables

are ratio and the number of frames per seconds (F/s).

Ratio denotes the ratio between the clock pace of the

simulation and the real time needed for the simulation.

If its value is greater than 1 it means the simulation

time runs at a pace faster than real time, and vice versa.

Obviously, with the chosen time-step and computing 4

time-steps before the screen is refreshed, the simulation

with all considered models can run on each of the used

hardware configurations at a pace faster than real time

or quite close to real time

Towards real-time simulation of deformable structures by means of co-rotational finite element formulation 11

6 Conclusions

VR- and AR-based environments call for adequate so-

lutions for the background physics of deformable ob-

jects in the framework of interactive simulations that

are supposed to run at the same speed as a real clock or

somewhat slower but without a delay that would affect

the aimed functionality. The essence of the presented

solution is the simplified co-rotational FEM formula-

tion. The implemented simplifications obviously make

the formulation unsuitable for certain analysis types,

such as buckling analysis. It was shown that, despite

the simplification implemented in the formulation, it

can cover the geometric nonlinearities to a significant

extent with accuracy that can be acceptable in certain

fields of applications. The examples with the shell el-

ements were focused on this aspect. But it should be

emphasized that the achievable accuracy strongly de-

pends on the character of geometrically nonlinear be-

havior. If it is dominated by large local rigid-body rota-

tions whereby strains and stress stiffening effects remain

small, then the presented formulation, which keeps the

element elastic behavior linear with respect to the co-

rotational reference frame, is expected to yield good ac-

curacy. The prerequisite for this is adequate FE mesh-

ing as the rigid-body rotation is accounted for element-

wise. This implies that the areas characterized by sub-

stantial gradients of rigid-body rotation upon deforma-

tion require finer meshing (of course, standard criteria

for FE meshing apply as well).

It was shown that FE models with several thousands

elements may run in real-time with conventional hard-

ware configurations. Some of the hardware configura-

tions used in the tests contain components that could

even be described as outdated. Shell elements are ob-

viously more demanding as they employ both nodal

translations and rotations, whereby handling the lat-

ter is numerically more expensive. A suitable solution

may be sought in solid-shell elements [43] that use only

translations as nodal degrees of freedom. Certainly, the

mentioned FE model sizes can be described as modest

compared to the FE models used in typical engineering

tasks, but the objectives of engineering FE models and

performed simulations are also substantially different.

In addition, the formulation can be enriched by the cou-

pled mesh technique, which together with the adequate

choice of the integration scheme (time-step size) and

solver provides vast options for the trade-off between

the numerical effort and accuracy.

As examples suggest, the approach can be applied for

the purpose of surgical simulation. If the real-time re-

quirement is not strict, i.e. if it is relaxed to, say, nu-

merically very efficient then the approach may also find

its application in MBS dynamics for flexible bodies ex-

periencing geometrically nonlinear deformation with re-

spect to the floating reference frame. This would elimi-

nate the need to superpose the large rigid-body motion

with the deformational motion, as the former is readily

incorporated into the co-rotational FEM. Finally, con-

sidering the computational power of commercially avail-

able hardware components (personal computers meant

here), the approach may represent a well-balanced al-

ternative to the available solutions for the current needs

of interactive simulations. Beside the considerations re-

lated to shells, further work should also tackle the chal-

lenges of implementing efficient solutions for material

nonlinearities, contact and tear and adequate estima-

tion of the average element rigid-body rotation.

Compliance with ethical standards

Conflict of Interest The authors declare that they

have no conflict of interest.

References

1. Argyris, J.: An excursion into large rotations. Computer

methods in applied mechanics and engineering 32(1-3),

85–155 (1982)

2. Bathe, K.J.: Finite element procedures. Klaus-Jurgen

Bathe (2006)

3. Bathe, K.J., Bolourchi, S.: Large displacement analysis of

three-dimensional beam structures. International Journal

for Numerical Methods in Engineering 14(7), 961–986

(1979)

4. Bro-Nielsen, M., Cotin, S.: Real-time volumetric de-

formable models for surgery simulation using finite el-

ements and condensation. In: Computer graphics forum,

vol. 15, pp. 57–66. Wiley Online Library (1996)

5. Capell, S., Green, S., Curless, B., Duchamp, T., Popovi´c,

Z.: Interactive skeleton-driven dynamic deformations. In:

ACM Transactions on Graphics (TOG), vol. 21, pp. 586–

593. ACM (2002)

6. Cerracchio, P., Gherlone, M., Tessler, A.: Real-time dis-

placement monitoring of a composite stiffened panel sub-

jected to mechanical and thermal loads. Meccanica

50(10), 2487–2496 (2015)

7. Choi, M.G., Ko, H.S.: Modal warping: Real-time sim-

ulation of large rotational deformation and manipula-

tion. IEEE Transactions on Visualization and Computer

Graphics 11(1), 91–101 (2005)

8. ´

Cojbaˇsi´c, ˇ

Z., Nikoli´c, V., Petrovi´c, E., Pavlovi´c, V.,

Tomi´c, M., Pavlovi´c, I., ´

Ciri´c, I.: A real time neural

network based finite element analysis of shell structure.

Facta Universitatis-Series Mechanical Engineering 12(2),

149–155 (2014)

9. Craig, R., Bampton, M.: Coupling of substructures for

dynamic analyses. AIAA journal 6(7), 1313–1319 (1968)

10. Cueto, E., Chinesta, F.: Real time simulation for com-

putational surgery: a review. Advanced Modeling and

Simulation in Engineering Sciences 1(1), 11 (2014)

12 Dragan Marinkovic et al.

11. Dogan, F., Serdar Celebi, M.: Real-time deformation sim-

ulation of non-linear viscoelastic soft tissues. Simulation

87(3), 179–187 (2011)

12. Dulong, J.L., Druesne, F., Villon, P.: A model reduc-

tion approach for real-time part deformation with non-

linear mechanical behavior. International Journal on In-

teractive Design and Manufacturing (IJIDeM) 1(4), 229

(2007)

13. Etzmuß, O., Keckeisen, M., Straßer, W.: A fast finite ele-

ment solution for cloth modelling. In: Computer Graphics

and Applications, 2003. Proceedings. 11th Pacific Confer-

ence on, pp. 244–251. IEEE (2003)

14. Felippa, C.A., Haugen, B.: A unified formulation of small-

strain corotational finite elements: I. theory. Computer

Methods in Applied Mechanics and Engineering 194(21-

24), 2285–2335 (2005)

15. Fiorentino, M., Monno, G., Uva, A.: Interactive touch

and see fem simulation using augmented reality. Interna-

tional Journal of Engineering Education 25(6) (2009)

16. Hambli, R., Chamekh, A., Salah, H.B.H.: Real-time de-

formation of structure using finite element and neural

networks in virtual reality applications. Finite elements

in analysis and design 42(11), 985–991 (2006)

17. Huang, J., Ong, S.K., Nee, A.Y.: Real-time finite element

structural analysis in augmented reality. Advances in

Engineering Software 87, 43–56 (2015)

18. Joldes, G.R., Wittek, A., Miller, K.: Real-time nonlin-

ear finite element computations on gpu–application to

neurosurgical simulation. Computer methods in applied

mechanics and engineering 199(49-52), 3305–3314 (2010)

19. Jung, W.Y., Han, S.C.: A refined element-based la-

grangian shell element for geometrically nonlinear analy-

sis of shell structures. Advances in Mechanical Engineer-

ing 7(4), 1687814015581272 (2015)

20. Kalkkuhl, J., Hunt, K.J., Fritz, H.: Fem-based neural-

network approach to nonlinear modeling with application

to longitudinal vehicle dynamics control. IEEE Transac-

tions on Neural Networks 10(4), 885–897 (1999)

21. Kang, L., Zhang, Q., Wang, Z.: Linear and geometrically

nonlinear analysis of novel flat shell elements with rota-

tional degrees of freedom. Finite Elements in Analysis

and Design 45(5), 386–392 (2009)

22. Kawamura, K., Kobayashi, Y., Fujie, M.G.: Basic study

on real-time simulation using mass spring system for

robotic surgery. In: International Workshop on Medi-

cal Imaging and Virtual Reality, pp. 311–319. Springer

(2008)

23. Levy, R., Gal, E.: The geometric stiffness of thick shell

triangular finite elements for large rotations. Inter-

national journal for numerical methods in engineering

65(9), 1378–1402 (2006)

24. Marinkovi´c, D., K¨oppe, H., Gabbert, U.: Numerically effi-

cient finite element formulation for modeling active com-

posite laminates. Mechanics of Advanced Materials and

Structures 13(5), 379–392 (2006)

25. Marinkovi´c, D., Rama, G.: Co-rotational shell element for

numerical analysis of laminated piezoelectric composite

structures. Composites Part B: Engineering 125, 144–

156 (2017)

26. Marinkovi´c, D., Zehn, M.: Geometric stiffness matrix in

modal space for multibody analysis of flexible bodies with

moderate deformation. Proc. of ISMA, Leuven, Belgium

pp. 3013–3020 (2010)

27. Marinkovic, D., Zehn, M.: Consideration of stress stiff-

ening and material reorientation in modal space based

finite element solutions. 20(5) (2017)

28. Marinkovic, D., Zehn, M., Marinkovic, Z.: Finite element

formulations for effective computations of geometrically

nonlinear deformations. Advances in Engineering Soft-

ware 50, 3–11 (2012)

29. Marinkovi´c, D., Zehn, M., Marinkovi´c, Z.: Modal space

based solutions including geometric nonlinearities for

flexible multi-body systems. Civil-Comp Proceedings

100, 4–7 (2012)

30. Mosegaard J Herborg P, S.T.: A gpu accelerated spring

mass system for surgical simulation. Studies in health

technology and informatics 111, 342–348 (2005)

31. M¨uller, M., Dorsey, J., McMillan, L., Jagnow, R., Cut-

ler, B.: Stable real-time deformations. In: Proceedings of

the 2002 ACM SIGGRAPH/Eurographics symposium on

Computer animation, pp. 49–54. ACM (2002)

32. Nealen, A., M¨uller, M., Keiser, R., Boxerman, E., Carl-

son, M.: Physically based deformable models in computer

graphics. In: Computer graphics forum, vol. 25, pp. 809–

836. Wiley Online Library (2006)

33. Nestorovi´c, T., Marinkovi´c, D., Shabadi, S., Trajkov, M.:

User defined finite element for modeling and analysis of

active piezoelectric shell structures. Meccanica 49(8),

1763–1774 (2014)

34. Nguyen, V.A., Zehn, M., Marinkovi´c, D.: An efficient

co-rotational fem formulation using a projector matrix.

Facta Universitatis-Series Mechanical Engineering 14(2),

227–240 (2016)

35. Niroomandi, S., Alfaro, I., Cueto, E., Chinesta, F.: Real-

time deformable models of non-linear tissues by model

reduction techniques. Computer methods and programs

in biomedicine 91(3), 223–231 (2008)

36. Niroomandi, S., Alfaro, I., Gonzalez, D., Cueto, E.,

Chinesta, F.: Real-time simulation of surgery by reduced-

order modeling and x-fem techniques. International jour-

nal for numerical methods in biomedical engineering

28(5), 574–588 (2012)

37. Nutti, B., Marinkovi´c, D.: An approach to efficient fem

simulations on graphics processing units using cuda.

Facta Universitatis-Series Mechanical Engineering 12(1),

15–25 (2014)

38. Ordaz-Hernandez, K., Fischer, X., Bennis, F.: Model re-

duction technique for mechanical behaviour modelling:

efficiency criteria and validity domain assessment. Pro-

ceedings of the Institution of Mechanical Engineers, Part

C: Journal of Mechanical Engineering Science 222(3),

493–505 (2008)

39. Rama, G.: A 3-node piezoelectric shell element for linear

and geometrically nonlinear dynamic analysis of smart

structures. Facta Universitatis-Series Mechanical Engi-

neering 15(1), 31–44 (2017)

40. Scholz, P.: Softwareentwicklung eingebetteter systeme:

grundlagen, modellierung, qualit¨atssicherung. Springer-

Verlag (2006)

41. Schwertassek, R., Wallrapp, O.: Dynamik flexibler

Mehrk¨orpersysteme: Methoden der Mechanik zum rech-

nergest¨utzten Entwurf und zur Analyse mechatronischer

Systeme. Springer-Verlag (2017)

42. Sussman, T., Bathe, K.J.: 3d-shell elements for struc-

tures in large strains. Computers & Structures 122, 2–12

(2013)

43. Sze, K., Chan, W., Pian, T.: An eight-node hybrid-stress

solid-shell element for geometric non-linear analysis of

elastic shells. International Journal for Numerical Meth-

ods in Engineering 55(7), 853–878 (2002)