FE Finite Elements

FINITE ELEMENTS MODEL OF THE HUMAN BODY: GEOMETRY AND NON-LINEAR MATERIAL PROPERTIES
Niels C.C.M.Moes
Delft University of T echnology
Dept.OCP/IO/DE/ICA
C.C.M.Moes@IO.TUDelft.nl
Imre Horv´ath
Delft University of T echnology
Dept.OCP/IO/DE/ICA
I.Horvath@IO.TUDelft.nl
ABSTRACT
To optimize the shape of the interface between the skin and a sitting support a Finite Elements model is needed.The geometry of the model was based on(i)in vivo measured skin shape data of the up-per leg and the buttock area,and the position of a set of bony landmarks,/nr in vitro measured bone shape data of the Visible Human Dataset,(ii)aux-illary surfaces to bridge the skin and the bone data. First a closed surface FE model was constructed,fol-lowed by hexmeshing to obtain the solid model.The boundary and the contact conditionswere derived from the spatial,the continuum and the force appli-cation boundaries.To quantify the material proper-ties,a one parameter version of the Mooney model was used(neo-Hookean).The actual coefficient was estimated by reference with predicted values of the maximum pressure.
KEYWORDS
FEA,ergonomics,CAD,sitting,biomechanics
NOMENCLATURE
{C10,C01}Mooney coefficients
ect ectomorphic index
E Ergonomics Goodness Index
F sitting force
FE Finite Elements
F n maximum nodal force
G neo-Hookean elasticity constant
{I1,I2,I3}deformation invariants
m body mass
MCS measuring coordinate system
W strain energy function
WCS working coordinate system
λprinciple stretch ratio πPressure value
πϕPhysiologically acceptable pressure
ΠPressure distribution
ΠE Externally applied pressure distribution σstress
strain
If in the text the word mesh is used it always refers to a FE mesh and never to a geometric mesh.The following anatomical terms are used:
caudal,cranial towards the head
distal towards the feet
epicondyles landmarks indicating the width
of the knee
femur thigh bone
ischial tuberosity seating bone
SIAS Anterior Superior Iliac Spine SIPS Posterior Superior Iliac Spine
1.INTRODUCTION
Many consumer products are physically handled by the user.For such products the shape of the con-tact area is among the most important aspects of product design.To generate the optimal shape spe-cific ergonomics requirements are needed.Such re-quirements usually relate to the pressure distribition during the expected product usage.Often such re-quirements are available in the form of specific mea-sures,guidelines or rules that are applicable for a specific product group,type of usage,or user group (Bullinger&Solf,1979;Sanders&McCormick, 1993).The development of an engine that gener-ates suggestions for the shape of the contact area for a large range of products needs a fundamental approach,based on fundamental,including medical and physiological,criteria.
Such approach calls for a deep knowledge as well as
a more general formulation of how different shapes reflect the fulfillment of basic criteria.Such knowl-edge is not yet available to a suffiently deep level.
The objective of this research is to apply shape and pressure optimization,controlled for improving the ergonomics quality of an initial shape of the contact area.The elaboration will be done for a sitting sup-port and standardized circumstances.
Any optimization process modifies a set of working variables,and needs an objective functional to con-trol the direction of the modification.A typical func-tional for this context must be build on a set of algo-rithms that quantify the deviation from a set of pre-defined ergonomics criteria.The actual functional be used in this research is called the Ergonomics Goodness Index,E,which depends on the differ-ence between the loaded model characteristics and the corresponding physiological linits for the same characteristics,for instance the internal stress.The E is not discussed in the current paper,but explained fully in(Moes,2001b).Its validity for this optimiza-tion process was inherently assumed.
The physical body of the initial configuration to be optimized is represented as a valid1FE model.If a pressure distribition is applied to the FE model, the model will deform.Now two situations can be compared.First,the physical body is supported by a given shape,which results in a pressure distribition,Π.Second,the FE model is loaded by the sameΠ, which results in a deformed model,outwardly rep-resented as a set of geometric points coinciding with the surface nodes of the contact area.The assump-tion was made that the shape of the support is not significantly different from the shape that is derived from the mentioned surface points.
The goal of this part of the research is the construc-tion of the FE model.including(i)the geometric definition of anatomical components the model,(ii) the scaling and the assembly of the components,(iii) the generation of the surface elements,(iv)the gen-eration of the solid mesh,and(v)the validation of the model for material characteristics.
Earlier attempts to build a FE model of the human buttock area were reported by,among others,(Chow &Odell,1978).They constructed a axisymmetric linearly elastic model,Young’s modulus15kPa,to calculate the internal stress distribution for several types of load.But to elaborate the last,high deform-1Representative for the user group.ing part of the analysis,the value of the Young’s modulus was increased to300kPa.The material was considered almost incompressible(Poisson’s ra-tio0.49).They found for a300Nflat frictionless sitting load that the maximum pressure at the inter-face was ca.200kPa,and von Mises stress contours up to400.The hydrostatic pressure caused minimal deformation,while the pressure gradient caused the von Mises stress.
Another linear analysis was carried out by(Todd& Tacker,1994),who tried to relate the interface pres-sure with the internal pressures at the ischium.They constructed the geometry from MRI images.The support was a model of a cushion,which was inti-mately integrated with the buttock model.The lin-ear Young’s modulus for the male and the female subjects were64.8kPa and47.5kPa The pressure values at the interface and the element just beneath the ischial tuberosity were17kPa and74kPa for the male subject,and15kPa and205kPa for the female subject.
The procedures were elaborated and implemented for a sitting support.The initial,not optimized sup-port was a horizontal,flat and hard surface,without support for the arms and the back.This configu-ration is also used to check the validity of the FE model,see section3.5.So far only the areas of the upper leg and the buttock area are included in the model of the body,and the arm rests and the back-rest omitted.
2.CONCEPTUAL SOLUTION
2.1.Reasoning model
Since no theoretical relationships are known that can be used to generate an ergonomically,and even more specific for pressure distribution criteria,op-timal shape the application of optimization proce-dures is inevitable.Therefore the basic reason-ing model is based on structural optimization ually structural optimization procedures consider modifications of mechanical structures(the independent variables such as the number of sup-porting beams,their dimensions,cross section and positioning)that support a specific type of load,with the eyemark of optimizing a specific criterion such as the total mass or weight.In the current context the selected variable for modification is the externally applied pressure distribition,while the optimization criterion is a compilation of the differences of the calculated internal(continuum)pressure values with
the pressure values that are acceptable from a phys-iological and medical point of view.Because of the continuum-mechanical relationships between stress and deformation the optimization of the pressure distribition leads to an external product shape,that coincides with the shape of the contact area.Now the essential assumption is done that this shape rep-resents the ergonomics optimal(as far as ergonomics considers pressure distribition)shape for the given circumstances.
2.2.General overview of the
optimization procedures
The configuration to be optimized was defined as a person sitting upright on a support.Only the body parts that are most relevant for the pressure distribi-tion were considered:the upper leg and the buttock area.Assuming symmetry only one body half was considered.
The shape optimization process contains the follow-ing steps,seefigure1.
1.In thefirst step the shape of the unloaded,
thus undeformed,body,Ξ0,is generated from
a vague geometric model for a set of body char-
acteristics.
2.The initial shape of the support,and thus the ini-
tial,deformed shape of the contacting surface of the body,isflat and horizontal,ΞF.The accord-ing pressure distribition is called theflat pres-sure distribution,ΠF.
3.The variables to be modified are the externally
applied pressure values(in the contact area).
The variables that define the E are the phys-iologically acceptable internal pressure values,πϕ,and the actual internal pressure values.
4.The direction of the pressure modification is cal-
culated for each of the FE nodes that were se-lected for external pressure application(contact nodes).
5.The configuration is modified to improve the er-
gonomics quality.
2.3.Vague geometric model of the
body
Since human properties such as anthropometry and force exertion are subject to uncertainty they are usually described in terms of statistics.A user group can therefore not be defined by crisp measures,and the shape of the unloaded human body must be for-mulated by vague expressions.
Although properties of a user group are statisti-cally distributed,and any analysis that is based on such data should consider vagueness,the FEA of a vaguely defined FE model is,however,not yet pos-sible with currently available techniques.
The usual escape to such problems is to analyse a set of statistically distributed instances,which has, as a matter of fact,nothing to do with analysis that is based on vagueness.The selection of such in-stances must then fulfill the requirements of strati-fication2.The stratification is carried out for the in-dependent variables,the body characteristics.From such a stratified sample of the predicting variables the corresponding properties of the FE model are de-rived,such as the geometry,the material properties, the contact conditions etc,see section2.4.
In section3this type of modelling and the gener-ation of instances will be explained briefly.For a detailed discussion see(Moes,2001a;Moes et al., 2001).
2.4.FE model
The actual calculation of the deformation of the con-tinuum,the deformation of the shape of the con-tact area,and the internal pressure distribition for the progressively varying external pressure distribi-tion,ΠE within the contact area,is done by FEA.
A complete,perfect FE model includes all tissues that possibly participate in the continuum mechan-ical behaviour of the body,such as skin,muscles, tendons,ligaments and fat,and it describes the mor-phological properties.In addition each tissue is de-fined by geometry,material properties(anisotropic, non-linear,time-dependent),and contact conditions. In the current research a much simpler model is con-structed for several reasons.
•No data are available that describe the shape of the involved tissues.
•A perfect FE model would require an unaccept-able amount of CPU-time with hardware cur-rently available.
2A stratified sample has a frequency distribution that re-flects the statistical distribution functions of the independent variables,in this context the body characteristics.
Figure1Overview of the main procedures that were developed to optimize the shape of the product.
•The construction of a simpler model is already a very complex job,that requires much modelling to be done by hand.
The FE model that developed for this research con-sists of three tissues:skin,bone and soft tissue.The soft tissue replaces all deformable tissues between the skin and the bony part of the model.The mo-bility of the joints between the different bones is re-duced to zero.
2.5.Validation of the FE model
In order to keep the FE model within manageble lim-its of complexity and required processing time of the FEA,the model includes a number of approxima-tions,assumptions and simplifications.Therefore a validation of the FE model is needed.The devel-oped validation procedure considers the result of a standardized load on the pressure distribition in the contact area and the deformation.In this research the deformation by aflat support is selected,result-ing in a corresponding pressure distribition,ΠF.A regression model predictingΠF for a set of body characteristics,is discussed in(Moes,2000c;Moes, 2000a).
2.6.Optimization procedure
In general any optimization procedure optimizes a set of independent variables of a starting configu-ration to obtain better values for a criterion vari-able,which is usually an objective functional.In the current research the starting configuration is the FE model loaded by aflat,rigid,horizontal support.The independent variables are the external pressure val-
ues caused by the support.This pressure is exerted
in the contact area.The objective functional is the E. The E is a functional of the deviations of the actual
continuum pressure values fromΠϕ←{πϕ}i.The definition of the E is discussed in detail in(Moes,
2001b).Although its definition allows for many criteria,currently only the pressure within the con-tinuum is considered.The optimization procedure varies the pressure in the contact area so that the E improves.The selected search method is according to the gradient search.The direction of the search is thus derived from the gradient∇E={∂E/∂πE}i. The search is done in sufficiently small steps so that the region of the optimum can be found.The search is terminated if the E does not improve for continued search in the currently defined direction.Then it is decided if a new search will be performed in a dif-ferent direction,or if the optimization process will be terminated.
2.7.Postprocessing for surface
generation
After termination the spatial positions of the nodes of the part of the FE model,that represent the con-tact area,will be processed further to create the sur-face of the sitting support.The method to describe the surface makes use of singularity hints(Horv´a th &Vergeest,1998).These singularities are used to generate the natural and physical surfaces,that are needed for the actual physical prototyping.
3.DETAILED SOLUTION
This section discusses the details of the construc-tion of the FE model.Only marginal attention is paid to the procedures to create a vague geometric model and a geometric instance,or to the details of the postprocessing stage for prototyping,which have been reported in(Moes,2001b;Moes,2001a;Moes et al.,2001).
3.1.Generation of geometry
The geometry of the FE model consists of four ge-ometric bodies,the skin,the femur,the sacrum and the pelvis.
Although it is in the line of the current research to generate the shape of the skin from a vague model for a specific set of body characteristics,a different approach was followed.Because the vague mod-elling program still needs further optimization,espe-cially for the construction of the distribution vectors for concave curvatures,the geometry of the skin was therefore obtained from the measured data of one of the subjects(male).
Since measuring the actual shape of the bones of the subjects requires invasive techniques such as CT-scanning(X-radiation,which shows medical risc), or MRI-scanning(expensive and slow),a different approach was followed.The data of the Visible Hu-man Project(VHP,1997)were used to obtain the shape of the pelvis,the sacrum and the femur of the male VHP-subject,for which permission was granted by the National Library of Medicine.
The shape of the bones is thus given as a set of coun-tour lines lying in a transversal3plane.The thick-ness of the slices of the original data set is1mm. but depending on the local curvatures the density of the contours may be decreased.
The surface geometry was used to create the FE surface meshes of the four components(skin,three bony parts).These surface meshes were combined to one surface model,and processed further to a solid mesh.
3.2.Surface meshing
In order to apply solid meshing,the surface mesh of the overall model must enclose a volume,while showing no defects such as holes,free edges or crossing elements.The actual continuum,where the 3Perpendicular to the length axis of the body internal pressure values will be evaluated for the E, is the soft tissue enclosed by the femur and the skin. At the distal end(knee)the continuum is closed by a
cross sectioning surface,see A infigure17(thisfig-ure shows thefinished surface FE model).At the caudal boundary the skin surface ends in the end curve(b),which is a closed curve along the iliac
crest,the inguinal ligament,the pubic symphysis, the buttock fold,and via the SIPS back to the iliac crest.The boundary between the soft tissue of the pelvis and the abdominal region is a surface coin-
ciding with the pelvic inlet surface(C).The medial boundary is the intersection with the medio-sagittal surface x=0(D).
Since the bodies originate from different subjects, the surfaces will probably notfit exactly.There-fore auxillary surfaces must be generatedfinalize the closing operation.
3.3.Solid meshing
The solid meshing of the continuum will be done using a hexahedral element mesh.The HEXMESH option of the FE modeller4applies a hexmesh within
the volume,leaving after thefirst pass a gap with the inner aspect of the surface mesh.This gap is then filled with elements in afinal meshing operation. The parameters that are used to create the hexmesh
will be discussed in section4.7.
When the solid mesh is generated,due to the proper-ties of the FE modeller all information about the sur-
face elements has been lost,including the sets defi-nitions,the morphological structure and the original tissues the elements belong to.Therefore,to apply
the boundary conditions for the FEA,a number of sets of elements have to be generated.
The set of the skin elements,the skin set is needed
to apply the mechanical properties of the skin.A subset of the skin set is needed to apply for con-tact conditions.The set of the elements enclosing the volume,that was formerly occupied by the bony
parts will be used to apply thefixation in space(ge-ometric boundary condition).Actaully the surface nodes of the bone set will befixed in space.This set is called the bone set.The medial closing plane,see
D infigure17,is used to inhibit transversal transla-tion(symmetry condition for the left and right body parts).The distal closing plane is used to inhibit
translation along the longitudinal axis of the femur 4Mentat2001,MSC
(local symmetry).
3.4.Material properties
Material properties were considered for soft tissue and skin.In the current model the bony parts are fixed in space and only used as part of the bound-ary of the soft tissue;thus its material properties are irrelevant.Since the soft tissue is subject to large de-formations,the elastic behaviour is necessarily non-linear.The simplest nonlinear approach is the neo-Hookean elastic behaviour
σ=G ((1+ )−(1+ )−2)
(1)
which is a special case of the Mooney-Rivlin defini-tion for incompressible (I 3=1)rubber materials
W =C 10(I 1−3)+C 01(I 2−3)
(2)
where I 1,I 2and I 3are the invariants I 1= i =1,2,3λ2
i ,I 2=λ1λ2+λ2λ3+λ3λ1,and I 3=λ21λ22λ23.W is the strain energy function and λ1,2,3are the principle stretch ratios.
If this is elaborated for the elements below the ischial tuberosity,assuming uni-axial compression
and λ3the principle stretch ratio in the z -direction,then λ3=λ,and λ1=λ2=1/√
λ.The strain energy function can then be written as W =
C 10(λ21+λ22+λ23)=C 10(λ2)+2/λ−3.The
coefficient C 10=G/2,where G is the shear mod-ulus for λ≈1.The engineering stress function σ=dW/dλ=G (λ−1/λ),and the true stress (Cauchy stress)τ=σ/(1/λ)=G (λ2−1/λ).In figure 2the function W ,σand τare shown for G =15kPa and 0.2≤λ≤2.
(Rollh¨a user,1950)reported the elongation of the skin without dermis constituents for several age groups.In figure 3the stress-strain curves and the derived Young’s modulus (σ/ )are shown,as they were derived from the graphs of the original paper.(Todd &Tacker,1994)describe a method to ob-tain an estimation of the effective normal Young’s Modulus for the sitting load of the tissues below the ischial tuberosities while the buttock were unsup-ported.The lowest value was found for the supine position:11.9kPa for a female and 15.2kPa for a male subject.The largest value was found for the normal seated position:47.5kPa and 64.8kPa re-spectively.
The mechanical properties of the human skin vary with the relative humidity of the environment,the wetness of the skin,the direction of the applied
load,
Figure 2The strain energy,the engineering stress and
the Cauchy stress functions.
age,and even the season of the year.Moreover,the layers of the skin vary strongly in this respect.The elastic behaviour is non-linear:the skin stiff-ness increases with the applied tension by a factor of 50(!);a mathematical formulation for the coeffi-cient of elasticity was derived by (Manschot &Brak-kee,1987b;Manschot &Brakkee,1987a).They found that the elastic characteristics of the skin de-pend on the direction of the strain (anisotropic),and that many other variables have a significant influ-ence,such as the season of the year.The charac-teristics of the underlying soft tissues were not in-vestigated.
3.5.Model validation
Firstly the support is brought into contact with the body until the total supporting force is balanced with the body weight of the upper body including the upper legs.The deformation of the body follows the stages of increasing contact force.The result-ing pressure distribition in the contsct area usually shows a few areas of high pressure.The pressure distribition in these areas will be used to optimize the FE model for its material properties.
The Mooney-Rivlin constants C 10and C 01of the soft tissue and the skin will be approached using an optimization procedure.This procedure com-
Figure3The stress-strain relationships of the skin for different age groups(Rollh¨a user,1950).
pares the calculated pressure distribition in the con-tact area with the predicted values of the maximum pressure and the pressure gradient below the region of the ischial tuberosities as they were described in (Moes,2000c).These values are predicted by a set of body characteristics for aflat,hard and horizontal support,which allows using the FE model,after ap-plying theflat deformation,for the materials model validation.
The criterion for the optimization of the elastic prop-erties of the model evaluates the difference between the calculated and the predicted values of the maxi-mum pressure and the pressure gradient.
3.6.Shape optimization
Firstly a small modification of the external pressure distribition is calculated that improves E.Then this pressure correction is iteratively applied until one of the termination criteria is fulfilled.After each cor-rection step the corresponding shape is calculated by FEA.Finally the optimized shape of the support is derived from the position of the external nodes of the elements,that are enabled for contact conditions.4.IMPLEMENTATION
4.1.Geometry of the skin
The shape of the skin of a male subject was scanned and defined by a set of longitudinal lines.The main relevant(to predict the sitting force,the maximum pressure and the pressure gradient)body character-istic of the subject are the body mass,m,and the ectomorphic index,ect.For the subject selected for this research:m=77kg and ect=6.These mea-surements included the position of a set of palpa-ble bony landmarks:the lower aspect of the ischial tuberosity,the SIASes,the major trochanter,and the lateral and medial epicondyles of the femur,seefig-ure4.For a detailed discussion see(Moes,
2000b). Figure4Example of the measured lines and landmarks.
4.2.Geometry of the bones
The femur was scanned using slices 1856to 2341of the data set of the male subject (VHP,1997),the pelvis bones slices 1733to 1949,and the sacrum 1760to 1907.The program used for the scan-ning (SurfDriver,3.5.5)is described in (Moody &Lozanoff,1999).Figure 5shows an example of a slice,selected at level 1810.The green and the blue curves show the scanned contours of the right and the left pelvis bones,the red curve the contours of the sacrum.At this level the femur is not
visible.
Figure 5An example of a photograph of a slice of the
Visible Human Project.The slice is number 1810from the male set.
The distance between the selected contours varied from 1to 5mm,depending on the observed shape singularities.The vertices of the contours were mea-sured by hand.Care was taken to include sufficient points to determine significant shape
singularities.
Figure 6The scanned contour lines of the femur.
Figure 6shows the scanned lines of the femur.The head,the major trochanter and the knee areas are
shown at the right.The knee area has been scanned,but in the final model it is omitted,since its bound-ary coincides with the distal end of the skin model.The scan density was increased at the level of the head to improve the accuracy.Along the shaft the scan density was 1–2cm −1The auxillary,logitudi-nal lines were added by the ACIS4.0formatter for geomtric solid models,that are used for import in the FE
modeller.
Figure 7The scanned lines of the pelvis bones.
Figure 7shows a lateral view of the scanned lines of the right pelvis bones,containing the ala (iliac bone),the pubic bone and the ischial tuberosity.The right side shows the pubis and the SIAS.At the bot-tom the ischial tuberosity is indicated,and at the left side the greater sciatic notch.The acetabulum can been recognized by the increased curvature in the lower middle area,see the encircled area.Since the obturator foramen is actually closed by a strong liga-mentous structure,it has been omitted in the model.Figure 8shows the scanned lines of the sacrum.The promontorium,which is the upper front part,and the coccyx are indicated.Also in the sacrum model not relevant cavities has been omitted.
4.3.Geometric preprocessing of the
data
Geometric preprocessing of the generated shape data of the skin was needed to prepare the data for the FE modeller.
Figure8The scanned lines of the sacrum.•Measured points that typically resulted from uncontrolled body movements of the observer must be removed.
•Incidentally a small correction for crossing scan lines was necessary.The correction simply means exchanging the crossing parts.
•A check is needed whether the measured data set is sufficient and complete.This check considers (i)the distribution density of the scan lines,and (ii)the extend of each contour(running from the knee to the depth of the buttock fold).An unsuf-ficient density of the scan lines is usually caused by tiredness of the investigator or the subject during the measurements.In either case doing additional measurements does not make sense. The extend of a scan line is unsufficient if it was not continued until the end line B(figure17).•A coordinate transformation is applied so that the measurement coordinate system(MCS)co-incides with the working coordinate systen (WCS).Since the thickness of the skin was as-sumed3mm,the origin of the WCS was defined 3mm above5the midpoint of the measured lo-cations of the ischial tuberosities.The x-axis runs from left to right,parallel to the line that connects the SIASes.The y-axis runs in sagit-tal direction,and the z-axis in upward direction. In(Moes et al.,2001)a detailed explanation is given.
5In the+z direction.
•The bony parts must be correctly positioned(i) to each other,and(ii)to the skin.This position-ing was performed by translation,rotation and scaling.The pelvis rotation angle6is set to ca 30
◦,which conforms the average position dur-ing sitting upright(Moes,1998),seefigure
9.
Figure9The pelvis rotation angle was set to30
◦.
The SIAS and the lower aspect of the ischial tuberosity are positioned in register with the cor-responding measured landmarks.The centres of the spheres of the head of the femur and the acetabulum must coincide.The epicondyles (knee)of the femur are positioned in register with the corresponding measured landmarks.To optimize thefit with the skin shape the bones were scaled in different directions.A detailed discussion is given in(Moes,2001a;Moes et al., 2001)
•For the exchange of the graphics(geometric) data from the CAD program to the FE modeller the ACIS4.0format for solid body definitions was used.The conversion of the solid,ACIS-formatted data to geometric curves will be per-formed within the FE modeller.
4.4.Generation of FE surface meshes
of the parts
In the FE modeller the imported edges of the geo-metric solid models were converted to closed curves, corresponding to the original measured curves. Then a surface was generated by lofting the scanned lines(Mentat:SKIN-mode).From this geometric surface a surface mesh was generated.For compli-cated areas such as skin folds,that show increased curvature,the surface mesh was partly corrected by hand or even completely(re)generated by hand di-rectly from the vertices(points)of the curves.
6The angle of rotation rotation of the pelvis is defined as the angle between the plane through the two SIASes and the front edge of the pubis and the z=0-plane.。