CN106202738B - Method for building up based on hyperelastic solids phase behaviour articular cartilage tow phase model - Google Patents

Method for building up based on hyperelastic solids phase behaviour articular cartilage tow phase model Download PDF

Info

Publication number
CN106202738B
CN106202738B CN201610555178.1A CN201610555178A CN106202738B CN 106202738 B CN106202738 B CN 106202738B CN 201610555178 A CN201610555178 A CN 201610555178A CN 106202738 B CN106202738 B CN 106202738B
Authority
CN
China
Prior art keywords
phase
articular cartilage
formula
model
equation
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Expired - Fee Related
Application number
CN201610555178.1A
Other languages
Chinese (zh)
Other versions
CN106202738A (en
Inventor
王沫楠
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Harbin University of Science and Technology
Original Assignee
Harbin University of Science and Technology
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Harbin University of Science and Technology filed Critical Harbin University of Science and Technology
Priority to CN201610555178.1A priority Critical patent/CN106202738B/en
Publication of CN106202738A publication Critical patent/CN106202738A/en
Application granted granted Critical
Publication of CN106202738B publication Critical patent/CN106202738B/en
Expired - Fee Related legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Prostheses (AREA)

Abstract

Based on the method for building up of hyperelastic solids phase behaviour articular cartilage tow phase model, the present invention relates to the method for building up of articular cartilage tow phase model.In order to solve the problems, such as that the solid phase by cartilage is regarded as the model of linear elastic materials for having very big difference with actual conditions when articular cartilage generation moderate finite deformation, even failing.The present invention obtains the geometrical model for the tissue that needs are studied according to the CT data reconstructions in joint first, then the initial data that the element number of model and node coordinate input as deformation calculation is obtained, the articular cartilage tow phase model based on theory of mixtures and the governing equation based on v p variables are established, articular cartilage mechanical balance equation is then established and carries out FEM calculation.The present invention is applied to the foundation and emulation of articular cartilage tow phase model.

Description

Method for building up based on hyperelastic solids phase behaviour articular cartilage tow phase model
Technical field
The present invention relates to the method for building up of articular cartilage tow phase model.
Background technology
Articular cartilage is mainly made up of water and compound organic matter.Compound organic matter mainly includes collagenous fibres and albumen is more Sugar, mainly contain water in interstitial fluid, the flowings of these fluids not only to the mechanical property important of cartilage, also with nothing The nutriment transmission close relation of vascular tissue.Meanwhile articular cartilage makes joint have extremely low coefficient of friction, plays Good lubrication.The ability that the lubrication of articular cartilage is moved to movable joint and bears load has extremely important meaning Justice.Therefore, the mechanical characteristic of articular cartilage is studied, recognizes its stress in motion process, strain conditions and its during exercise The changing rule of each mechanical quantity, the exploitation to medical research, clinical diagnosis and novel biomaterial, there is very important theory Meaning and practical value.
At research initial stage of articular cartilage mechanical modeling, people regard articular cartilage as isotropy linear elasticity single-phase material Material.In following period of time, researcher receives viscoelastic single-phase articular cartilage model, however, passing through this constitutive relation Obtained simulation result always has very big difference with the actual conditions and experimental data of cartilage.At present, articular cartilage is considered as Tow phase model.But mostly regard the solid phase of cartilage as linear elastic materials in these researchs, linear elastic model is only to small deformation Effectively, the solid phase of cartilage is regarded as line bullet can produce moderate finite deformation during the application of power the problem of for articular cartilage Property material model will have very big difference with actual conditions, or even failure.
The content of the invention
The model that the present invention is regarded as linear elastic materials for the solid phase solved cartilage produces larger change for articular cartilage There is the problem of very big difference, even fail with actual conditions when shape.
Based on the method for building up of hyperelastic solids phase behaviour articular cartilage tow phase model, comprise the following steps：
Step 1, geometrical model and mesh generation：
Using joint as research object, by obtaining the CT data in joint, obtained data are exported and stored up in dicom format Deposit in a computer；The region of articular cartilage in CT data is separated from its hetero-organization by function of image segmentation in MIMICS Out, the geometrical model for the tissue that needs are studied is obtained by three-dimensional reconstruction；
From DICOM data to generation grid model and obtain model element number and node coordinate it is defeated as deformation calculation The initial data entered, acquisition process are as shown in Figure 1；
Step 2, establish the articular cartilage tow phase model based on theory of mixtures and the governing equation based on v-p variables, bag Include following steps：
Step 2.1, the tow phase model for establishing articular cartilage：
Articular cartilage is considered as to the two-phase medium mixture being made up of hyperelastic solids and perfect fluid, and two-phase has solely The vertical characteristics of motion, solid phase is represented with s, f represents liquid phase, then φsRepresent solid phase volume fraction, φfRepresent liquid phase Volume fraction；In initial configuration, the initial volume fraction of two-phase is φs0、φf0, the initial density of two-phase isAnd It is uniformly distributed；
Obtained by the volume integral numerical expression under saturation：
φs+φf=1 (1)
Mass balance equation is：
Wherein,Represent vector differential operator symbol；vsIt is solid phase velocity, vfIt is liquid phase velocity；
In quasi-static problem, momentum balance equation is converted into the equation of static equilibrium, and acceleration a is equal to zero；In this research Ignore the effect of external volume power, in two-phase mixture, the frictional resistance shape as caused by interstitial fluid flows through porosu solid Into momentum-exchange be Ps、Pf, and have
Ps=-Pf=K ' (vf-vs) (3)
In formula, K ' is diffusional resistance coefficient；
In formula, κ is cartilage permeability；Solid phase momentum is Ps, liquid phase momentum Pf；
Make σsAnd σfThe respectively cauchy stress tensor of solid phase and liquid phase, articular cartilage is in the momentum balance side of current configuration Cheng Wei：
The power of articular cartilage elastic solid (Hookean body) phase and inviscid flow body phase is obtained by the Incoercibility and mechanical characteristic of cartilage Learning constitutive equation is：
σs=-φspI+σe (7)
σf=-φfpI (8)
Wherein, I is unit matrix, and p is pressure；σeDeflection phase one with solid phase is also referred to as the effective stress of solid phase The elastic stress tensor of cause；
It is derived from more than by the sheet of mass balance equation (2), momentum balance equation (5) (6) and two-phase mixture Structure equation (7) (8) constitutes the tow phase model of articular cartilage；
2nd Piola-Kirchhoff stress tensors S is understood according to the definition of the 2nd Piola-Kirchhoff stress tensors With elastic stress tensor σeRelation be：
σe=J-1FSFT=FSFT (9)
Wherein, F is deformation gradient matrix；J is elastic volume ratio, is the Jacobian of deformation gradient matrix F；S is 2nd Piola-Kirchhoff stress tensors, the parameter are solved in step 2.3, and articular cartilage is obtained by S connection Super-elasticity two-phase medium model；
In finite deformation, it is contemplated that the change of locus before and after object receiving force, make t=0 moment object in space institute The region V occupied0For initial configuration；Current configuration is expressed as in current time t, the area of space V occupied by object；
Using initial configuration as reference configuration, the change of volume elements, bin in finite deformation can be drawn using material formulation method, Volume in the volume and reference configuration of unit material in i.e. current configuration meets following relation：
DV=JdV0 (10)
Step 2.2, establish hyperelastic model：
The hyperelastic model of articular cartilage represents by Helmholtz strain energy functions, i.e.,：
α0、α1、α2, β be material mechanics parameter, from the fitting result of mathematical modeling and experimental data；I1、I2、I3 It is C first, second, third main invariant respectively；
For incompressible material, I3=det (C)=1, the hyperelastic model after simplifying are：
Step 2.3, solve the 2nd Piola-Kirchhoff stress tensors S：
The hyperelastic solids phase of articular cartilage, ess-strain no longer meet linear relationship, and elastic matrix has with strain tensor Close, the 2nd Piola-Kirchhoff stress tensor S and Green strain tensors E is the stress using initial configuration as reference configuration And strain tensor, during an instantaneous deformation, constitutive relation is：DS=DdE；
Green strain tensors E is：
Elasticity tensor D is expressed as：
In formula, the element in DSi″j″Represent the element of i-th " row jth " row in S, Ek″l″Represent in E Kth " row l " row element；
Order：
2α1+4α2(1+E22+E33)=G
2α1+4α2(1+E11+E22)=H
2α1+4α2(1+E11+E33)=B
Symmetrical matrix D is obtained, the symmetric part in the matrix is represented with sym：
Step 2.4, establish the governing equation based on v-p variables
Access speed v, pressure p are that unknown quantity establishes the finite element equation based on v-p variables in this research, and then express and close Save governing equation of the cartilage based on v-p variables；
In order to eliminate the speed of fluid phase, first formula (3) and formula (8) are taken in formula (6), obtained by formula (4)：
Under saturation conditions formula (1), formula (20) is brought into formula (2), obtains the relational expression for only including solid phase velocity：
Next can be obtained by momentum conservation equation (5) (6) and constitutive equation (7) (8) simultaneous：
Governing equation based on v-p variables is just constituted by formula (21) and formula (22)；
Step 3, establish articular cartilage mechanical balance equation and carry out FEM calculation：
Step 3.1, establish boundary condition：
Fluid flux is defined asThe border of governing equation then based on v-p variables Condition is：
In formula, spatial point when x refers at a time t in current configuration；For known solid phase displacement,For total load Load forces,For known pressure, n*It is border Γ outer normal vector；ΓpAnd ΓQLiquid phase volume domain Ω is represented respectivelyfSide Boundary,And ΓTSolid phase volume domain Ω is represented respectivelysBorder, and border meets condition：
Step 3.2, construction shape function：
Define finite element tetrahedron element shape function be：
N=[Ni′·I3×3 Nj′·I3×3 Nm′·I3×3 Np′·I3×3] (44)
Wherein,
In formula, VeIt is the volume of tetrahedron element, ai′、bi′、ci′、di′For with Ni′Corresponding coefficient；
Nj′、Nm′、Np′With Ni′It is similar, have
Step 3.3, establish articular cartilage mechanical balance equation：
Cell type based on finite element solving is tetrahedron element, constructs interpolating function；Pass through the control based on v-p variables Equation (21) processed, (22) and boundary condition (23)-(26), articular cartilage mechanical balance is established using the golden weighted residual method of gal the Liao Dynasty Equation；
Step 3.4, articular cartilage mechanical balance equation is solved, realize the emulation of articular cartilage tow phase model.
The present invention has advantages below：
1st, current articular cartilage is considered as two models.But mostly the solid phase of cartilage is regarded as in these researchs Linear elastic materials, linear elastic model are only effective to small deformation.Larger change can be produced during the application of power for articular cartilage The problem of shape, present invention proposition establish articular cartilage tow phase model based on hyperelastic solids phase behaviour, realize articular cartilage two The foundation and emulation of phase model.
2nd, the present invention is directed to the design feature of articular cartilage, and articular cartilage is described as by incompressible non-viscous liq The biphasic mixture mutually and with super-elasticity, the solid phase of transverse isotropy formed, and can also under conditions of deformation Cause the infiltrative change of cartilaginous tissue.
3rd, the super elastic characteristics of solid phase are defined according to Helmholtz strain energy functions, liquid phase is defined as perfect fluid, The model can describe the non-linear of articular cartilage, Incoercibility and permeability.Base is obtained by the golden weighted residual method of gal the Liao Dynasty In the finite element system equilibrium equation of v-p variables, numerical computations are carried out to equilibrium equation by finite difference calculus.
4th, the present invention produces the situation of moderate finite deformation, the pass that the present invention establishes for articular cartilage during the application of power Section cartilage tow phase model also can accurately react the situation of actual articular cartilage, compared with the situation of actual articular cartilage, Error is less than 13%.
Brief description of the drawings
Fig. 1 is to show from DICOM data to the flow of generation grid model and the element number and node coordinate that obtain model It is intended to.
Fig. 2 is the stress-speed experiments curve and FEM Numerical Simulation comparison diagram of uniaxial compression.
Embodiment
Embodiment one：
Based on the method for building up of hyperelastic solids phase behaviour articular cartilage tow phase model, comprise the following steps：
Step 1, geometrical model and mesh generation：
Using joint as research object, by obtaining the CT data in joint, obtained data are exported and stored up in dicom format Deposit in a computer；The region of articular cartilage in CT data is separated from its hetero-organization by function of image segmentation in MIMICS Out, the geometrical model for the tissue that needs are studied is obtained by three-dimensional reconstruction；
From DICOM data to generation grid model and obtain model element number and node coordinate it is defeated as deformation calculation The initial data entered, acquisition process are as shown in Figure 1；
Step 2, establish the articular cartilage tow phase model based on theory of mixtures and the governing equation based on v-p variables, bag Include following steps：
Step 2.1, the tow phase model for establishing articular cartilage：
Articular cartilage is considered as to the two-phase medium mixture being made up of hyperelastic solids and perfect fluid, and two-phase has solely The vertical characteristics of motion, solid phase is represented with s, f represents liquid phase, then φsRepresent solid phase volume fraction, φfRepresent liquid phase Volume fraction；In initial configuration, the initial volume fraction of two-phase is φs0、φf0, the initial density of two-phase isAnd It is uniformly distributed；
Obtained by the volume integral numerical expression under saturation：
φs+φf=1 (1)
Mass balance equation is：
Wherein,Represent vector differential operator symbol；vsIt is solid phase velocity, vfIt is liquid phase velocity；
In quasi-static problem, momentum balance equation is converted into the equation of static equilibrium, and acceleration a is equal to zero；In this research Ignore the effect of external volume power, in two-phase mixture, the frictional resistance shape as caused by interstitial fluid flows through porosu solid Into momentum-exchange be Ps、Pf, and have
Ps=-Pf=K ' (vf-vs) (3)
In formula, K ' is diffusional resistance coefficient；
In formula, κ is cartilage permeability；Solid phase momentum is Ps, liquid phase momentum Pf；
Make σsAnd σfThe respectively cauchy stress tensor of solid phase and liquid phase, articular cartilage is in the momentum balance side of current configuration Cheng Wei：
The power of articular cartilage elastic solid (Hookean body) phase and inviscid flow body phase is obtained by the Incoercibility and mechanical characteristic of cartilage Learning constitutive equation is：
σs=-φspI+σe (7)
σf=-φfpI (8)
Wherein, I is unit matrix, and p is pressure；σeDeflection phase one with solid phase is also referred to as the effective stress of solid phase The elastic stress tensor of cause；
It is derived from more than by the sheet of mass balance equation (2), momentum balance equation (5) (6) and two-phase mixture Structure equation (7) (8) constitutes the tow phase model of articular cartilage；
2nd Piola-Kirchhoff stress tensors S is understood according to the definition of the 2nd Piola-Kirchhoff stress tensors With elastic stress tensor σeRelation be：
σe=J-1FSFT=FSFT (9)
Wherein, F is deformation gradient matrix；J is elastic volume ratio, is the Jacobian of deformation gradient matrix F；S is 2nd Piola-Kirchhoff stress tensors, the parameter are solved in 2.3, and the superlastic of articular cartilage is obtained by S connection Property two-phase medium model；
In finite deformation, it is contemplated that the change of locus before and after object receiving force, make t=0 moment object in space institute The region V occupied0For initial configuration；Current configuration is expressed as in current time t, the area of space V occupied by object；
Using initial configuration as reference configuration, the change of volume elements, bin in finite deformation can be drawn using material formulation method, Volume in the volume and reference configuration of unit material in i.e. current configuration meets following relation：
DV=JdV0 (10)
Step 2.2, establish hyperelastic model：
The hyperelastic model of articular cartilage represents by Helmholtz strain energy functions, i.e.,：
α0、α1、α2, β be material mechanics parameter, from the fitting result of mathematical modeling and experimental data；I1、I2、I3 It is C first, second, third main invariant respectively；
For incompressible material, I3=det (C)=1, the hyperelastic model after simplifying are：
Step 2.3, solve the 2nd Piola-Kirchhoff stress tensors S：
The hyperelastic solids phase of articular cartilage, ess-strain no longer meet linear relationship, and elastic matrix has with strain tensor Close, the 2nd Piola-Kirchhoff stress tensor S and Green strain tensors E is the stress using initial configuration as reference configuration And strain tensor, during an instantaneous deformation, constitutive relation is：DS=DdE；
Green strain tensors E is：
Elasticity tensor D is expressed as：
In formula, the element in DSi″j″Represent the element of i-th " row jth " row in S, Ek″l″Represent in E Kth " row l " row element；
Order：
2α1+4α2(1+E22+E33)=G
2α1+4α2(1+E11+E22)=H
2α1+4α2(1+E11+E33)=B
Symmetrical matrix D is obtained, the symmetric part in the matrix is represented with sym：
Step 2.4, establish the governing equation based on v-p variables
Access speed v, pressure p are that unknown quantity establishes the finite element equation based on v-p variables in this research, and then express and close Save governing equation of the cartilage based on v-p variables；
In order to eliminate the speed of fluid phase, first formula (3) and formula (8) are taken in formula (6), obtained by formula (4)：
Under saturation conditions formula (1), formula (20) is brought into formula (2), obtains the relational expression for only including solid phase velocity：
Next can be obtained by momentum conservation equation (5) (6) and constitutive equation (7) (8) simultaneous：
Governing equation based on v-p variables is just constituted by formula (21) and formula (22)；
Step 3, establish articular cartilage mechanical balance equation and carry out FEM calculation：
Step 3.1, establish boundary condition：
Fluid flux is defined asThe border of governing equation then based on v-p variables Condition is：
In formula, spatial point when x refers at a time t in current configuration；For known solid phase displacement,For total load Load forces,For known pressure, n*It is border Γ outer normal vector；ΓpAnd ΓQLiquid phase volume domain Ω is represented respectivelyfSide Boundary,And ΓTSolid phase volume domain Ω is represented respectivelysBorder, and border meets condition：
Step 3.2, construction shape function：
Define finite element tetrahedron element shape function be：
N=[Ni′·I3×3 Nj′·I3×3 Nm′·I3×3 Np′·I3×3] (44)
Wherein,
In formula, VeIt is the volume of tetrahedron element, ai′、bi′、ci′、di′For with Ni′Corresponding coefficient；
Nj′、Nm′、Np′With Ni′It is similar, have
Step 3.3, establish articular cartilage mechanical balance equation：
Cell type based on finite element solving is tetrahedron element, constructs interpolating function；Pass through the control based on v-p variables Equation (21) processed, (22) and boundary condition (23)-(26), articular cartilage mechanical balance is established using the golden weighted residual method of gal the Liao Dynasty Equation；
Step 3.4, articular cartilage mechanical balance equation is solved, realize the emulation of articular cartilage tow phase model.
Embodiment two：
The detailed process that step 3.3 described in present embodiment establishes articular cartilage mechanical balance equation is as follows：
The cell type of finite element solving is tetrahedron element, selects power of the shape function as the differential equation and boundary condition Function, it is respectively N to make the weight function in Weighted Residual integral expressionJWith-NJ, corresponding Weighted Residual expression formula is obtained, is passed through Weighted Residual integral expression after the definition of divergence theorem of Gauss and gross tractive effort is simplified is：
Problem Areas is discretized into unit form, row interpolation is entered to speed and pressure：
V=Nve (30)
P=Npe (31)
Wherein, νeRepresent tetrahedron element speed, peRepresent tetrahedron element pressure；
Then the matrix form of the Weighted Residual expression formula on a finite element tetrahedron element grid is：
In formula, enFor a unit vector, due to enIt is non-zero constant vector, therefore can directly disappears；
Order：
Obtain final articular cartilage mechanical balance equation：
Yv+M=f (34)
In formula, Y is a capacity matrix, and M is the non-linear of related to elastic stress tensor nodal force in solid phase Flexible vector, f are plus load force vectors.
Other steps and parameter are identical with embodiment one.
Embodiment three：
The detailed process of solution articular cartilage mechanical balance equation described in step 3.4 described in present embodiment is as follows：
Numerical computations are carried out to articular cartilage mechanical balance equation (34) by finite difference calculus：
Method of addition is used first in time domain, turns to several time points by time t is discrete, i.e.,：
0=t0<t1<…<tn<tn+1<…<tN=t (35)
Incremental time is expressed as：
△tn+1=tn+1-tn (36)
Using complete Lagrangian method, initial time t is made0=0 configuration is as reference configuration, and the reference configuration is not Change with time and change；In tn+1Moment formula (34) is write as：
Yn+1vn+1+Mn+1=fn+1 (37)
By Newton-Raphson method, M linearised form is obtained by following recursive form：
In formula, Kn+1It is the tangent stiffness matrix of solid phase, Kn+1It isVariation；I represents iterations；
Iteration displacement is determined by trapezoidal methodRecurrence relation be：
In formula, u is modal displacement,For tn+1The modal displacement of moment ith iteration；ω is the time integral specified Parameter, 0≤ω≤1, and as ω >=1/2 is implicit integration, equation unconditional stability；
The incremental representation of speed is：
In the ith iteration of speed, equation (34) can be write as：
It can be obtained by (38) (39) (40) (41)：
In △ tn+1In time step, iteration is repeated to formula (42), until speed incrementMeet convergence criterion； In a time step, the convergence criterion of iteration is：
When current step inner iteration restrains generation, next time step is arrived using obtained solution as initial value renewal It is interior, and iteration is carried out until convergence occurs again；Iterate and carry out numerical computations, finally realize articular cartilage two The emulation of phase model.The process for iterating and carrying out numerical computations programs in fact under Visual Studio platform C language environment It is existing.
Other steps and parameter are identical with embodiment one or two.
Embodiment four：
Iterating described in present embodiment step 3.4 and carry out numerical computations iterative flow it is as follows：
(a) in △ t=tn+1-tnInitial velocity v known to definition on incrementnWith initial displacement un, and have
(b) iterations is set：I=1
(c) useThe matrix and remaining unknown quantity come in accounting equation (42), resolves this linear algebraic equation and obtains First speed increment
(d) terminate first time iteration, update displacement and velocity vector：
(e) iterations is updated to i
(f) useMatrix and vector are updated, solution linear equation obtains speed increment
(g) current iteration, renewal speed and motion vector are terminated：
Return to step (e) and continue new iteration, until the increment of speedMeet convergence criterion.
Other steps and parameter are identical with one of embodiment one to three.
Embodiment five：
V described in present embodiment step 3.2e、ai′、bi′、ci′、di′Concrete form it is as follows
In formula, (x 'i′, y 'i′, z 'i′)、(x′j′, y 'j′, z 'j′)、(x′m′, y 'm′, z 'm′)、(x′p′, y 'p′, z 'p′) respectively For four apex coordinates of tetrahedron element.
Other steps and parameter are identical with one of embodiment one to four.
Embodiment six：
Elastic volume ratio J=1 described in present embodiment step 2.1, because articular cartilage is considered as incompressible shaping Material, so J=1.
Other steps and parameter are identical with one of embodiment one to five.
Embodiment seven：
Cartilage permeability described in present embodiment step 2.1Wherein L is constant, subscript " 0 " Represent not deformed preceding initial configuration.Articular cartilage not only has super-elasticity, and can also cause tissue under conditions of deformation Infiltrative change.
Other steps and parameter are identical with one of embodiment one to six.
Embodiment eight：
Described in present embodimentIn φs0=0.2, κ0=2.519 × 10-15m4/ (Ns), L= 0.0848。
Other steps and parameter are identical with one of embodiment one to seven.
Embodiment nine：
α described in present embodiment step 2.20=0.1084Mpa, α1=0.592Mpa, α2=0.0846Mpa.
Other steps and parameter are identical with one of embodiment one to eight.
The present invention produces the situation of moderate finite deformation, the joint that the present invention establishes for articular cartilage during the application of power Cartilage tow phase model also can accurately react the situation of actual articular cartilage.As shown in Fig. 2 by with the reality under the same terms Comparative result is tested, Finite element analysis results all coincide on numerical value and trend with experimental data curve, and " experimental result " is in figure The experimental data of actual articular cartilage, " simulation result " are the foundation of the present invention based on hyperelastic solids phase behaviour articular cartilage The emulation data of tow phase model.Found by the contrast of simulation result and empirical curve：Based on hyperelastic solids phase behaviour joint Cartilage tow phase model can accurately describe the mechanical characteristic of articular cartilage, and compared with the situation of actual articular cartilage, error is low In 13%, the accuracy of model and the validity of finite element program are demonstrated.

Claims (9)

1. the method for building up based on hyperelastic solids phase behaviour articular cartilage tow phase model, it is characterised in that comprise the following steps：
Step 1, geometrical model and mesh generation：
Using joint as research object, by obtaining the CT data in joint, obtained data are exported and are stored in dicom format In computer；The region of articular cartilage in CT data is isolated from its hetero-organization by function of image segmentation in MIMICS Come, the geometrical model for the tissue that needs are studied is obtained by three-dimensional reconstruction；
To generation grid model and the element number of model is obtained from DICOM data and node coordinate inputs as deformation calculation Initial data；
Step 2, articular cartilage tow phase model of the foundation based on theory of mixtures and the governing equation based on v-p variables, including with Lower step：
Step 2.1, the tow phase model for establishing articular cartilage：
Articular cartilage is considered as to the two-phase medium mixture being made up of hyperelastic solids and perfect fluid, and two-phase have it is independent The characteristics of motion, solid phase is represented with s, f represents liquid phase, then φsRepresent solid phase volume fraction, φfRepresent liquid phase volume Fraction；In initial configuration, the initial volume fraction of two-phase is φs0、φf0, the initial density of two-phase isAnd uniformly Distribution；
Obtained by the volume integral numerical expression under saturation：
φs+φf=1 (1)
Mass balance equation is：
Wherein,Represent vector differential operator symbol；vsIt is solid phase velocity, vfIt is liquid phase velocity；
In quasi-static problem, momentum balance equation is converted into the equation of static equilibrium, and acceleration a is equal to zero；Ignore external volume The effect of power, in two-phase mixture, momentum-exchange that the frictional resistance as caused by interstitial fluid flows through porosu solid is formed For Ps、Pf, and have
Ps=-Pf=K ' (vf-vs) (3)
In formula, K ' is diffusional resistance coefficient；
In formula, κ is cartilage permeability；Solid phase momentum is Ps, liquid phase momentum Pf；
Make σsAnd σfRespectively the cauchy stress tensor of solid phase and liquid phase, articular cartilage are in the momentum balance equation of current configuration：
The mechanics sheet of articular cartilage elastic solid (Hookean body) phase and inviscid flow body phase is obtained by the Incoercibility and mechanical characteristic of cartilage Structure equation is：
σs=-φspI+σe (7)
σf=-φfpI (8)
Wherein, I is unit matrix, and p is pressure；σeFor the effective stress of solid phase；
It is derived from more than by this structure side of mass balance equation (2), momentum balance equation (5) (6) and two-phase mixture Journey (7) (8) constitutes the tow phase model of articular cartilage；
2nd Piola-Kirchhoff stress tensors S and bullet are understood according to the definition of the 2nd Piola-Kirchhoff stress tensors Property stress tensor σeRelation be：
σe=J-1FSFT=FSFT (9)
Wherein, F is deformation gradient matrix；J is elastic volume ratio, is the Jacobian of deformation gradient matrix F；S is second Piola-Kirchhoff stress tensors；
In finite deformation, it is contemplated that the change of locus before and after object receiving force, make t=0 moment objects occupied by space Region V0For initial configuration；Current configuration is expressed as in current time t, the area of space V occupied by object；
Using initial configuration as reference configuration, the change of volume elements, bin in finite deformation can be drawn using material formulation method, that is, is worked as Volume in the volume and reference configuration of unit material in preceding configuration meets following relation：
DV=JdV0 (10)
Step 2.2, establish hyperelastic model：
The hyperelastic model of articular cartilage represents by Helmholtz strain energy functions, i.e.,：
α0、α1、α2, β be material mechanics parameter；I1、I2、I3It is C first, second, third main invariant respectively；
For incompressible material, I3=det (C)=1, the hyperelastic model after simplifying are：
Step 2.3, solve the 2nd Piola-Kirchhoff stress tensors S：
The hyperelastic solids phase of articular cartilage, ess-strain no longer meet linear relationship, and elastic matrix is relevant with strain tensor, the Two Piola-Kirchhoff stress tensor S and Green strain tensors E are using initial configuration as the stress of reference configuration and answered Become tensor, during an instantaneous deformation, constitutive relation is：DS=DdE；
Green strain tensors E is：
Elasticity tensor D is expressed as：
In formula, the element in D Represent the element of i-th " row jth " row in S, Ek″l″Represent the kth in E " The element of row l " row；
Order：
2α1+4α2(1+E22+E33)=G
2α1+4α2(1+E11+E22)=H
2α1+4α2(1+E11+E33)=B
Symmetrical matrix D is obtained, the symmetric part in the matrix is represented with sym：
Step 2.4, establish the governing equation based on v-p variables
Access speed v, pressure p are that unknown quantity establishes the finite element equation based on v-p variables, and then express articular cartilage and be based on v- The governing equation of p variables；
First formula (3) and formula (8) are taken in formula (6), obtained by formula (4)：
Under saturation conditions formula (1), formula (20) is brought into formula (2), obtains the relational expression for only including solid phase velocity：
Next can be obtained by momentum conservation equation (5) (6) and constitutive equation (7) (8) simultaneous：
Governing equation based on v-p variables is just constituted by formula (21) and formula (22)；
Step 3, establish articular cartilage mechanical balance equation and carry out FEM calculation：
Step 3.1, establish boundary condition：
Fluid flux is defined asThe boundary condition of governing equation then based on v-p variables For：
In formula, spatial point when x refers at a time t in current configuration；For known solid phase displacement,For total loading force,For known pressure, n*It is border Γ outer normal vector；ΓpAnd ΓQLiquid phase volume domain Ω is represented respectivelyfBorder, And ΓTSolid phase volume domain Ω is represented respectivelysBorder, and border meets condition：
Step 3.2, construction shape function：
Define finite element tetrahedron element shape function be：
N=[Ni′·I3×3 Nj′·I3×3 Nm′·I3×3 Np′·I3×3] (44)
Wherein,
In formula, VeIt is the volume of tetrahedron element, ai′、bi′、ci′、di′For with Ni′Corresponding coefficient；
Nj′、Nm′、Np′With Ni′It is similar, have
Step 3.3, establish articular cartilage mechanical balance equation：
Cell type based on finite element solving is tetrahedron element, constructs interpolating function；Pass through the controlling party based on v-p variables Journey (21), (22) and boundary condition (23)-(26), articular cartilage mechanical balance side is established using the golden weighted residual method of gal the Liao Dynasty Journey；
Step 3.4, articular cartilage mechanical balance equation is solved, realize the emulation of articular cartilage tow phase model.
2. the method for building up according to claim 1 based on hyperelastic solids phase behaviour articular cartilage tow phase model, it is special Sign be step 3.3 establish articular cartilage mechanical balance equation detailed process it is as follows：
The cell type of finite element solving is tetrahedron element, selects shape function as the differential equation and the power letter of boundary condition Number, it is respectively N to make the weight function in Weighted Residual integral expressionJWith-NJ, Weighted Residual expression formula is obtained, passes through Gauss divergence Weighted Residual integral expression after the definition of theorem and gross tractive effort is simplified is：
Problem Areas is discretized into unit form, row interpolation is entered to speed and pressure：
V=Nve (30)
P=Npe (31)
Wherein, νeRepresent tetrahedron element speed, peRepresent tetrahedron element pressure；
Then the matrix form of the Weighted Residual expression formula on a finite element tetrahedron element grid is：
In formula, enFor a unit vector；
Order：
Obtain final articular cartilage mechanical balance equation：
Yv+M=f (34)
In formula, Y is a capacity matrix, and M is the nonlinear elasticity of related to elastic stress tensor nodal force in solid phase Vector, f are plus load force vectors.
3. the method for building up according to claim 2 based on hyperelastic solids phase behaviour articular cartilage tow phase model, it is special Sign is that the detailed process of the solution articular cartilage mechanical balance equation described in step 3.4 is as follows：
Numerical computations are carried out to articular cartilage mechanical balance equation (34) by finite difference calculus：
Method of addition is used first in time domain, turns to several time points by time t is discrete, i.e.,：
0=t0<t1<…<tn<tn+1<…<tN=t (35)
Incremental time is expressed as：
△tn+1=tn+1-tn (36)
Using complete Lagrangian method, initial time t is made0=0 configuration is as reference configuration, and the reference configuration is not at any time Between change and change；In tn+1Moment formula (34) is write as：
Yn+1vn+1+Mn+1=fn+1 (37)
By Newton-Raphson method, M linearised form is obtained by following recursive form：
In formula, Kn+1It is the tangent stiffness matrix of solid phase, Kn+1It isVariation；I represents iterations；
Iteration displacement is determined by trapezoidal methodRecurrence relation be：
In formula, u is modal displacement,For tn+1The modal displacement of moment ith iteration；ω is the time integral parameter specified, 0 ≤ ω≤1, and as ω >=1/2, is implicit integration, equation unconditional stability；
The incremental representation of speed is：
In the ith iteration of speed, equation (34) can be write as：
It can be obtained by (38) (39) (40) (41)：
In △ tn+1In time step, iteration is repeated to formula (42), until speed incrementMeet convergence criterion；One In individual time step, the convergence criterion of iteration is：
When current step inner iteration restrains generation, arrived obtained solution as initial value renewal in next time step, and And iteration is carried out until convergence occurs again；Iterate and carry out numerical computations, finally realize articular cartilage two-phase mould The emulation of type.
4. the method for building up according to claim 3 based on hyperelastic solids phase behaviour articular cartilage tow phase model, it is special Sign be iterating described in step 3.4 and carry out numerical computations iterative flow it is as follows：
(a) in △ t=tn+1-tnInitial velocity v known to definition on incrementnWith initial displacement un, and have
(b) iterations is set：I=1
(c) useThe matrix and remaining unknown quantity come in accounting equation (42), resolves this linear algebraic equation and obtains first Individual speed increment
(d) terminate first time iteration, update displacement and velocity vector：
(e) iterations is updated to i
(f) useMatrix and vector are updated, solution linear equation obtains speed increment
(g) current iteration, renewal speed and motion vector are terminated：
Return to step (e) and continue new iteration, until the increment of speedMeet convergence criterion.
5. the foundation side based on hyperelastic solids phase behaviour articular cartilage tow phase model according to claim 1,2,3 or 4 Method, it is characterised in that the V described in step 3.2e、ai′、bi′、ci′、di′Concrete form it is as follows
In formula, (x 'i′, y 'i′, z 'i′)、(x′j′, y 'j′, z 'j′)、(x′m′, y 'm′, z 'm′)、(x′p′, y 'p′, z 'p′) it is respectively four Four apex coordinates of face body unit.
6. the method for building up according to claim 5 based on hyperelastic solids phase behaviour articular cartilage tow phase model, it is special Sign is the elastic volume ratio J=1 described in step 2.1.
7. the method for building up according to claim 6 based on hyperelastic solids phase behaviour articular cartilage tow phase model, its feature exist Cartilage permeability described in step 2.1Wherein L is constant, subscript " 0 " represent it is not deformed before just Begin to configure.
8. the method for building up according to claim 7 based on hyperelastic solids phase behaviour articular cartilage tow phase model, its feature exist InIn φs0=0.2, κ0=2.519 × 10-15m4/ (Ns), L=0.0848.
9. the method for building up according to claim 8 based on hyperelastic solids phase behaviour articular cartilage tow phase model, its feature exist α described in step 2.20=0.1084Mpa, α1=0.592Mpa, α2=0.0846Mpa.
CN201610555178.1A 2016-07-14 2016-07-14 Method for building up based on hyperelastic solids phase behaviour articular cartilage tow phase model Expired - Fee Related CN106202738B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201610555178.1A CN106202738B (en) 2016-07-14 2016-07-14 Method for building up based on hyperelastic solids phase behaviour articular cartilage tow phase model

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201610555178.1A CN106202738B (en) 2016-07-14 2016-07-14 Method for building up based on hyperelastic solids phase behaviour articular cartilage tow phase model

Publications (2)

Publication Number Publication Date
CN106202738A CN106202738A (en) 2016-12-07
CN106202738B true CN106202738B (en) 2017-12-19

Family

ID=57476137

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201610555178.1A Expired - Fee Related CN106202738B (en) 2016-07-14 2016-07-14 Method for building up based on hyperelastic solids phase behaviour articular cartilage tow phase model

Country Status (1)

Country Link
CN (1) CN106202738B (en)

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN107808037A (en) * 2017-10-10 2018-03-16 哈尔滨理工大学 A kind of Modeling Calculation method of the articular cartilage based on machine direction
CN109033742B (en) * 2018-06-21 2019-06-21 哈尔滨理工大学 It is a kind of for simulating the method for building up of the shear-deformable hyperelastic model of soft tissue
CN109544673B (en) * 2018-10-19 2023-06-23 瑞梦德医药科技(北京)有限公司 Three-dimensional imaging method and system for cartilage segmentation of medical image

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8598144B1 (en) * 2005-11-10 2013-12-03 Reyn Pharma, Llc Method of administering hyaluronan formulation for the amelioration of osteoarthritis
WO2009052562A1 (en) * 2007-10-23 2009-04-30 Commonwealth Scientific And Industrial Research Organisation Automatic segmentation of articular cartilage in mr images
CN103310072B (en) * 2013-06-28 2015-12-23 哈尔滨理工大学 Based on the biomechanical properties finite element analysing system of force feedback

Also Published As

Publication number Publication date
CN106202738A (en) 2016-12-07

Similar Documents

Publication Publication Date Title
Fritzen et al. Two-stage data-driven homogenization for nonlinear solids using a reduced order model
Kamensky et al. An immersogeometric variational framework for fluid–structure interaction: Application to bioprosthetic heart valves
Teran et al. Finite volume methods for the simulation of skeletal muscle
Hartmann et al. Discontinuous Galerkin methods for computational aerodynamics—3D adaptive flow simulation with the DLR PADGE code
Li et al. Stable Orthotropic Materials.
CN110298105A (en) The CCPDI-IMPM method of saturated porous media analysis on Large Deformation
CN106202738B (en) Method for building up based on hyperelastic solids phase behaviour articular cartilage tow phase model
CN107133397B (en) A method of two-way wind-structure interaction is carried out to biovalve based on ALE methods
Furukawa et al. Accurate cyclic plastic analysis using a neural network material model
CN101496028A (en) Method of simulating deformable object using geometrically motivated model
Belnoue et al. A rapid multi-scale design tool for the prediction of wrinkle defect formation in composite components
CN106021977B (en) Subcutaneus adipose tissue biomethanics modeling method based on linear elasticity and hyperelastic model
Zhang et al. A total‐Lagrangian material point method for coupled growth and massive deformation of incompressible soft materials
Faisal et al. Computational study of the elastic properties of Rheum rhabarbarum tissues via surrogate models of tissue geometry
Koyanagi et al. Direct FE2 for simulating strain-rate dependent compressive failure of cylindrical CFRP
Wang From immersed boundary method to immersed continuum methods
CN106354954A (en) Three-dimensional mechanical modal simulation method based on hierarchical basis function
Bloxom Numerical simulation of the fluid-structure interaction of a surface effect ship bow seal
CN108710735A (en) A kind of mesh free soft tissue deformation analogy method of real-time, interactive
Zienkiewicz Computational mechanics today
Ge et al. Blending isogeometric and Lagrangian elements in three-dimensional analysis
Binesh et al. Elasto-plastic analysis of reinforced soils using mesh-free method
Zeng Numerical Simulations of the Two-phase flow and Fluid-Structure Interaction Problems with Adaptive Mesh Refinement
Wang A one-field fictitious domain method for fluid-structure interactions
Mazhar et al. A novel artificial neural network-based interface coupling approach for partitioned fluid–structure interaction problems

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant
CF01 Termination of patent right due to non-payment of annual fee
CF01 Termination of patent right due to non-payment of annual fee

Granted publication date: 20171219

Termination date: 20190714