From d7939b7882b564d77387f4da57277148eb3b362e Mon Sep 17 00:00:00 2001 From: dsc Date: Thu, 7 Nov 2024 12:42:26 +0100 Subject: [PATCH 1/2] Add trunk cervical stiffness --- Body/AAUHuman/Trunk/CervicalDiscStiffness.any | 229 ++++++++++++++++++ Body/AAUHuman/Trunk/TrunkModel.root.any | 1 + 2 files changed, 230 insertions(+) create mode 100644 Body/AAUHuman/Trunk/CervicalDiscStiffness.any diff --git a/Body/AAUHuman/Trunk/CervicalDiscStiffness.any b/Body/AAUHuman/Trunk/CervicalDiscStiffness.any new file mode 100644 index 000000000..1ab47c1b0 --- /dev/null +++ b/Body/AAUHuman/Trunk/CervicalDiscStiffness.any @@ -0,0 +1,229 @@ + AnyFolder CervicalDiscStiffness = { + + #if (BM_TRUNK_CERVICAL_LIGAMENTS == ON) + AnyFloat coef = 0.0; + #else + AnyFloat coef = 1.0; + #endif + + // Cervical coefficients (Nm/deg) + /* + Moroney SP, Schultz AB, Miller JA, Andersson GB. + Load-displacement properties of lower cervical spine motion segments. + J Biomech. 1988;21(9):769-79. doi: 10.1016/0021-9290(88)90285-0. + */ + + // Stiffness (Nm/degree) for disc segments (IVD) including longitudinal ligaments + AnyFloat StiffnessScaleFactor ??= 1.0; + AnyFloat AxialRot ??= 0.42*StiffnessScaleFactor; + AnyFloat Flex ??= 0.21*StiffnessScaleFactor; + AnyFloat Ext ??= 0.32*StiffnessScaleFactor; + AnyFloat Latbend ??= 0.33*StiffnessScaleFactor; + + // Stiffness contribution (Nm/degree) of the posterior arches and ligaments to the cumulative IVD stiffness + // difference of intact segment and disc segment from Moroney et al. + AnyFloat ligAxialRot ??= coef*0.74*StiffnessScaleFactor; + AnyFloat ligFlex ??= coef*0.22*StiffnessScaleFactor; + AnyFloat ligExt ??= coef*0.41*StiffnessScaleFactor; + AnyFloat ligLatbend ??= coef*0.35*StiffnessScaleFactor; + + #if BM_TRUNK_CERVICAL_DISC_STIFFNESS == _DISC_STIFFNESS_NONE_ + AnyFunPolynomial Flexfun = { + PolyCoef={{0, 0}}; + }; + AnyFunPolynomial &Extfun = Flexfun; + AnyFunPolynomial &ARfun = Flexfun; + AnyFunPolynomial &LBfun = Flexfun; + AnyFunPolynomial &LBfun_R = Flexfun; + AnyFunPolynomial &LBfun_L = Flexfun; + #endif + + #if BM_TRUNK_CERVICAL_DISC_STIFFNESS == _DISC_STIFFNESS_LINEAR_ + + AnyFunPolynomial Flexfun = { + PolyCoef??=-1*{{0, .Flex+.ligFlex}}; + }; + AnyFunPolynomial Extfun = { + PolyCoef??=-1*{{0, .Ext+.ligExt}}; + }; + AnyFunPolynomial ARfun = { + PolyCoef??=-1*{{0, .AxialRot+.ligAxialRot}}; + }; + AnyFunPolynomial LBfun = { + PolyCoef??=-1*{{0, .Latbend+.ligLatbend}}; //linear behavior, 0-centered, symmetric with opposite sign for R/L + }; + AnyFunPolynomial &LBfun_R = LBfun; + AnyFunPolynomial &LBfun_L = LBfun; + + #endif + + #if BM_TRUNK_CERVICAL_DISC_STIFFNESS == _DISC_STIFFNESS_NONLINEAR_ + + AnyInt CervicalStiffnessNonLinearWarning = warn(0,strformat("\n"+ + "------------------------------------------------------------------------------------------------------------------------ \n"+ + " `BM_TRUNK_CERVICAL_DISC_STIFFNESS` is set to `_DISC_STIFFNESS_NONLINEAR_`\n"+ + " Nonlinear function is not available currently and linear function is used instead.\n"+ + "------------------------------------------------------------------------------------------------------------------------ ")); + + AnyFunPolynomial Flexfun = { + PolyCoef??=-1*{{0, .Flex+.ligFlex}}; + }; + AnyFunPolynomial Extfun = { + PolyCoef??=-1*{{0, .Ext+.ligExt}}; + }; + AnyFunPolynomial ARfun = { + PolyCoef??=-1*{{0, .AxialRot+.ligAxialRot}}; + }; + AnyFunPolynomial LBfun = { + PolyCoef??=-1*{{0, .Latbend+.ligLatbend}}; //linear behavior, 0-centered, symmetric with opposite sign for R/L + }; + AnyFunPolynomial &LBfun_R = LBfun; + AnyFunPolynomial &LBfun_L = LBfun; + + #endif + + //T1C7 + AnyForce T1C7DiscMoment = { + // find out where it is flexion or extension + AnyFloat A = iffun(gtfun(rot.Pos[0],0), 0, 1); + AnyFloat B = iffun(gtfun(rot.Pos[1],0), 0, 1); + F = { + (A*.Flexfun(Angle[0])[0] + (1-A)*.Extfun(Angle[0])[0]), + .ARfun( Angle[1])[0], + (B*.LBfun_L(Angle[2])[0] + (1-B)*.LBfun_R(Angle[2])[0]) + }; + AnyFloat Angle = rot.Pos*180/pi; // Z,Y,X (FlexExt,AxRot,LatBend) + AnyKinRotational rot = { + AnyRefNode &rf0 = ....Segments.T1Seg.T1C7JntNode; + AnyRefNode &rf1 = ....Segments.C7Seg.T1C7JntNode; + Type=RotAxesAngles; + }; + }; + + //C7C6 + AnyForce C7C6DiscMoment = { + // find out where it is flexion or extension + AnyFloat A = iffun(gtfun(rot.Pos[0],0), 0, 1); + AnyFloat B = iffun(gtfun(rot.Pos[1],0), 0, 1); + F = { + (A*.Flexfun(Angle[0])[0] + (1-A)*.Extfun(Angle[0])[0]), + .ARfun( Angle[1])[0], + (B*.LBfun_L(Angle[2])[0] + (1-B)*.LBfun_R(Angle[2])[0]) + }; + AnyFloat Angle = rot.Pos*180/pi; // Z,Y,X (FlexExt,AxRot,LatBend) + AnyKinRotational rot = { + AnyRefNode &rf0 = ....Segments.C7Seg.C7C6JntNode; + AnyRefNode &rf1 = ....Segments.C6Seg.C7C6JntNode; + Type=RotAxesAngles; + }; + }; + + //C6C5 + AnyForce C6C5DiscMoment = { + // find out where it is flexion or extension + AnyFloat A = iffun(gtfun(rot.Pos[0],0), 0, 1); + AnyFloat B = iffun(gtfun(rot.Pos[1],0), 0, 1); + F = { + (A*.Flexfun(Angle[0])[0] + (1-A)*.Extfun(Angle[0])[0]), + .ARfun( Angle[1])[0], + (B*.LBfun_L(Angle[2])[0] + (1-B)*.LBfun_R(Angle[2])[0]) + }; + AnyFloat Angle = rot.Pos*180/pi; // Z,Y,X (FlexExt,AxRot,LatBend) + AnyKinRotational rot = { + AnyRefNode &rf0 = ....Segments.C6Seg.C6C5JntNode; + AnyRefNode &rf1 = ....Segments.C5Seg.C6C5JntNode; + Type=RotAxesAngles; + }; + }; + + //C5C4 + AnyForce C5C4DiscMoment = { + // find out where it is flexion or extension + AnyFloat A = iffun(gtfun(rot.Pos[0],0), 0, 1); + AnyFloat B = iffun(gtfun(rot.Pos[1],0), 0, 1); + F = { + (A*.Flexfun(Angle[0])[0] + (1-A)*.Extfun(Angle[0])[0]), + .ARfun( Angle[1])[0], + (B*.LBfun_L(Angle[2])[0] + (1-B)*.LBfun_R(Angle[2])[0]) + }; + AnyFloat Angle = rot.Pos*180/pi; // Z,Y,X (FlexExt,AxRot,LatBend) + AnyKinRotational rot = { + AnyRefNode &rf0 = ....Segments.C5Seg.C5C4JntNode; + AnyRefNode &rf1 = ....Segments.C4Seg.C5C4JntNode; + Type=RotAxesAngles; + }; + }; + + //C4C3 + AnyForce C4C3DiscMoment = { + // find out where it is flexion or extension + AnyFloat A = iffun(gtfun(rot.Pos[0],0), 0, 1); + AnyFloat B = iffun(gtfun(rot.Pos[1],0), 0, 1); + F = { + (A*.Flexfun(Angle[0])[0] + (1-A)*.Extfun(Angle[0])[0]), + .ARfun( Angle[1])[0], + (B*.LBfun_L(Angle[2])[0] + (1-B)*.LBfun_R(Angle[2])[0]) + }; + AnyFloat Angle = rot.Pos*180/pi; // Z,Y,X (FlexExt,AxRot,LatBend) + AnyKinRotational rot = { + AnyRefNode &rf0 = ....Segments.C4Seg.C4C3JntNode; + AnyRefNode &rf1 = ....Segments.C3Seg.C4C3JntNode; + Type=RotAxesAngles; + }; + }; + + //C3C2 + AnyForce C3C2DiscMoment = { + // find out where it is flexion or extension + AnyFloat A = iffun(gtfun(rot.Pos[0],0), 0, 1); + AnyFloat B = iffun(gtfun(rot.Pos[1],0), 0, 1); + F = { + (A*.Flexfun(Angle[0])[0] + (1-A)*.Extfun(Angle[0])[0]), + .ARfun( Angle[1])[0], + (B*.LBfun_L(Angle[2])[0] + (1-B)*.LBfun_R(Angle[2])[0]) + }; + AnyFloat Angle = rot.Pos*180/pi; // Z,Y,X (FlexExt,AxRot,LatBend) + AnyKinRotational rot = { + AnyRefNode &rf0 = ....Segments.C3Seg.C3C2JntNode; + AnyRefNode &rf1 = ....Segments.C2Seg.C3C2JntNode; + Type=RotAxesAngles; + }; + }; + + //C2C1 - The specimens in Moroney et al. did not contain C2-C1. Therefore, it is excluded. +// AnyForce C2C1DiscMoment = { +// // find out where it is flexion or extension +// AnyFloat A = iffun(gtfun(rot.Pos[0],0), 0, 1); +// AnyFloat B = iffun(gtfun(rot.Pos[1],0), 0, 1); +// F = { +// (A*.Flexfun(Angle[0])[0] + (1-A)*.Extfun(Angle[0])[0]), +// .ARfun( Angle[1])[0], +// (B*.LBfun_L(Angle[2])[0] + (1-B)*.LBfun_R(Angle[2])[0]) +// }; +// AnyFloat Angle = rot.Pos*180/pi; // Z,Y,X (FlexExt,AxRot,LatBend) +// AnyKinRotational rot = { +// AnyRefNode &rf0 = ....Segments.C2Seg.C2C1JntNode; +// AnyRefNode &rf1 = ....Segments.C1Seg.C2C1JntNode; +// Type=RotAxesAngles; +// }; +// }; + + //C1C0 - The specimens in Moroney et al. did not contain C1-C0. Therefore, it is excluded. +// AnyForce C1C0DiscMoment = { +// // find out where it is flexion or extension +// AnyFloat A = iffun(gtfun(rot.Pos[0],0), 0, 1); +// AnyFloat B = iffun(gtfun(rot.Pos[1],0), 0, 1); +// F = { +// (A*.Flexfun(Angle[0])[0] + (1-A)*.Extfun(Angle[0])[0]), +// .ARfun( Angle[1])[0], +// (B*.LBfun_L(Angle[2])[0] + (1-B)*.LBfun_R(Angle[2])[0]) +// }; +// AnyFloat Angle = rot.Pos*180/pi; // Z,Y,X (FlexExt,AxRot,LatBend) +// AnyKinRotational rot = { +// AnyRefNode &rf0 = ....Segments.C1Seg.C1C0JntNode; +// AnyRefNode &rf1 = ....Segments.SkullSeg.C1C0JntNode; +// Type=RotAxesAngles; +// }; +// }; + + }; diff --git a/Body/AAUHuman/Trunk/TrunkModel.root.any b/Body/AAUHuman/Trunk/TrunkModel.root.any index 493e8287f..43946dae1 100644 --- a/Body/AAUHuman/Trunk/TrunkModel.root.any +++ b/Body/AAUHuman/Trunk/TrunkModel.root.any @@ -156,6 +156,7 @@ AnyFolder ElasticElements = { #include "DiscStiffness.any" + #include "CervicalDiscStiffness.any" }; /* Other sections */ From 2619f0c5830c586ddc5e4f1749c12bdd050e251a Mon Sep 17 00:00:00 2001 From: dsc Date: Thu, 7 Nov 2024 12:51:19 +0100 Subject: [PATCH 2/2] Update changelog --- CHANGELOG.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 4b7c7d8c9..ec8c7f1ab 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -50,6 +50,11 @@ The default pelvis model used in all models have changed. The pelvis morphology will distribute the mass to the different segments based on whether they are marked as being part of the distribution. +* The cervical model now has linear stiffness coefficients. These can be enabled with the switch {bm_statement}`BM_TRUNK_CERVICAL_DISC_STIFFNESS`: + ``` + #define BM_TRUNK_CERVICAL_DISC_STIFFNESS _DISC_STIFFNESS_LINEAR_ + ``` + Only linear stiffness function is available for cervical discs currently. **Changed:**