ExpressionBasedCoordinateForce quadratic spring sign issue

[Pasha van Bijlert]

Hello,

I'm looking to add a quadratic spring to my model, which can passively suspend a segment in horizontal equilibrium. I came across ExpressionBasedCoordinateForce in the forums, and was wondering if it is possible to set the expression using absolute signs? Otherwise I run into the following issue, assuming arbitrary stiffness of 30 & damping of 3, and a rest angle of 0.2:

Code: Select all

<ExpressionBasedCoordinateForce name="rotational_spring">
	<!--Coordinate (name) to apply force to.-->
	<coordinate>joint</coordinate>
	<!--Expression of the force magnitude as a function of the coordinate value (q)and its time derivative (qdot). Note, expression cannot have any whitespace separating characters-->
	<expression>-30*((q-0.2)^2)-3*qdot</expression>
</ExpressionBasedCoordinateForce>

The above expression will give a negative moment/torque in all cases, but if q goes below -0.2, this leads to an unstable situation (where the spring torque will lead to greater spring torques). Does EBCF allow for absolute signs, for instance using:

Code: Select all

<expression>-30*(q-0.2)*(|q-0.2|)-3*qdot</expression>

I was also wondering if there are any example snippets of code to add such a force to a model via the MATLAB API?

Best wishes,
Pasha

[Thomas Uchida]

The documentation for the setExpression() method (https://simtk.org/api_docs/opensim/api_ ... bb9e14ccc3) states that "Expressions with C-mathematical operations such as +,-,*,/ and common functions: exp, pow, sqrt, sin, cos, tan, and so on are acceptable." I doubt the vertical bar character would work, but you might try fabs().

[Pasha van Bijlert]

Hello!

While fabs() throws an error, abs() works, thank you! The following snippet of code seems to correctly add a damped, rotational quadratic spring to a model via the Matlab API:

Code: Select all

spring1=ExpressionBasedCoordinateForce(joint1_coorname,'-30*(q+0.2)*abs(q+0.2)-3*qdot')
    model.addForce(spring1);
    model.finalizeConnections();
    

Here, the stiffness is 30 Nm/rad, damping is 3 Nms/rad, and the rest angle of the spring is at -0.2 rad. joint1_coorname is the coordinate name of the joint (for instance, hip_angle).

I double checked the spring moment using:

Code: Select all

state = model.initSystem();
joint_angle = model.updCoordinateSet().get('joint1_coorname'); 
joint_angle.setValue(state,phi);
spring1.calcExpressionForce(state)

Where phi is any arbitrary angle at which you know how much moment you would like your spring to produce.

Thank you for your help!
Pasha