Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

feat: Trunk cervical stiffness #34

Merged
merged 2 commits into from
Nov 7, 2024
Merged

feat: Trunk cervical stiffness #34

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
d7939b7
Add trunk cervical stiffness
divyaksh-chander Nov 7, 2024
2619f0c
Update changelog
divyaksh-chander Nov 7, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
229 changes: 229 additions & 0 deletions Body/AAUHuman/Trunk/CervicalDiscStiffness.any
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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;
// };
// };

};
1 change: 1 addition & 0 deletions Body/AAUHuman/Trunk/TrunkModel.root.any
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -156,6 +156,7 @@ AnyFolder ElasticElements = {


#include "DiscStiffness.any"
#include "CervicalDiscStiffness.any"

};
/* Other sections */
Expand Down
5 changes: 5 additions & 0 deletions CHANGELOG.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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:**
Expand Down