Add basic force elements Bushing and SpringDamperPtP

- Add force element Bushing and test HarmonicOscillator.
- Add force element SpringDamperPtP and test BoxSpringDamperPtP.

Force elements do not work for non solid Object3Ds for now because these are not yet considered in the equations of motion.

Modia/src/Modia3D_Interface.jl needs to be synchronized with Modia3D/src/ModiaInterface/model3D.jl asap.

Author: GerhardHippmann

src/Composition/ForceElements/Bushing.jl

@@ -0,0 +1,46 @@
+
+"""
+force = Bushing(; obj1, obj2, stiffness, damping)
+
+Return a `force` acting as bushing between `obj1::`[`Object3D`](@ref) and
+`obj2::`[`Object3D`](@ref). Vectors `stiffness` and `damping` define the
+stiffness resp. damping values in x, y and z direction of `obj1`. The
+orientation of `obj2` does no influence the resulting forces.
+"""
+mutable struct Bushing <: Modia3D.AbstractForceElement
+
+    obj1::Modia3D.AbstractObject3D
+    obj2::Modia3D.AbstractObject3D
+
+    stiffness::SVector{3,Float64}
+    damping::SVector{3,Float64}
+
+    function Bushing(; obj1::Modia3D.AbstractObject3D,
+                       obj2::Modia3D.AbstractObject3D,
+                       stiffness::AbstractVector = Modia3D.ZeroVector3D,
+                       damping::AbstractVector   = Modia3D.ZeroVector3D)
+
+        stiff = Modia3D.convertAndStripUnit(SVector{3,Float64}, u"N/m" , stiffness)
+        damp  = Modia3D.convertAndStripUnit(SVector{3,Float64}, u"N*s/m", damping)
+
+        return new(obj1, obj2, stiff, damp)
+
+    end
+end
+
+
+function initializeForceElement(force::Bushing)
+    return nothing
+end
+
+function evaluateForceElement(time, force::Bushing)
+    r12 = measFramePosition(force.obj2; frameOrig=force.obj1, frameCoord=force.obj1)
+    v12 = measFrameTransVelocity(force.obj2; frameOrig=force.obj1, frameCoord=force.obj1, frameObsrv=force.obj1)
+    f12 = force.stiffness .* r12 + force.damping .* v12
+    applyFrameForcePair!(force.obj2, force.obj1, f12; frameCoord=force.obj1)
+    return nothing
+end
+
+function terminateForceElement(force::Bushing)
+    return nothing
+end

src/Composition/ForceElements/SpringDamperPtP.jl

@@ -0,0 +1,45 @@
+
+"""
+force = SpringDamperPtP(; obj1, obj2, stiffness, damping)
+
+Return a `force` acting as point-to-point parallel spring-damper between
+`obj1::`[`Object3D`](@ref) and `obj2::`[`Object3D`](@ref). Values
+`stiffness` and `damping` define the stiffness resp. damping coefficient.
+"""
+mutable struct SpringDamperPtP <: Modia3D.AbstractForceElement
+
+    obj1::Modia3D.AbstractObject3D
+    obj2::Modia3D.AbstractObject3D
+
+    stiffness::Float64
+    damping::Float64
+
+    function SpringDamperPtP(; obj1::Modia3D.AbstractObject3D,
+                               obj2::Modia3D.AbstractObject3D,
+                               stiffness::Float64 = 0.0,
+                               damping::Float64 = 0.0)
+
+        stiff = Modia3D.convertAndStripUnit(Float64, u"N/m"  , stiffness)
+        damp  = Modia3D.convertAndStripUnit(Float64, u"N*s/m", damping)
+
+        return new(obj1, obj2, stiff, damp)
+
+    end
+end
+
+
+function initializeForceElement(force::SpringDamperPtP)
+    return nothing
+end
+
+function evaluateForceElement(time, force::SpringDamperPtP)
+    (pos, norm) = measFrameDistance(force.obj2; frameOrig=force.obj1)
+    vel = measFrameDistVelocity(force.obj2; frameOrig=force.obj1)
+    frc = force.stiffness * pos + force.damping * vel
+    f12 = frc * norm
+    applyFrameForcePair!(force.obj2, force.obj1, f12; frameCoord=force.obj1)
+end
+
+function terminateForceElement(force::SpringDamperPtP)
+    return nothing
+end

src/Composition/_module.jl

@@ -32,6 +32,8 @@ export Revolute, setAngle!, connect
 export Prismatic, setDistance!
 export FreeMotion
 
+export Bushing, SpringDamperPtP
+
 export multibodyResiduals!, setModiaJointVariables!
 
 export distanceAndAngles, distance, planarRotationAngle
@@ -93,12 +95,15 @@ include("object3D.jl")
 include("supportPoints.jl")
 include("superObjects.jl")
 
-# Joints must be included before scene, because joint memory is stored in scene
+# Joints and Force Elements must be included before scene, because memory is stored in scene
 include(joinpath("joints", "FreeMotion.jl"))
 include(joinpath("joints", "Fix.jl"))
 include(joinpath("joints", "Revolute.jl"))
 include(joinpath("joints", "Prismatic.jl"))
 
+include(joinpath("ForceElements", "Bushing.jl"))
+include(joinpath("ForceElements", "SpringDamperPtP.jl"))
+
 include("contactPairs.jl")
 include(joinpath(Modia3D.path, "src", "contactDetection", "ContactDetectionMPR", "ContactDetectionMPR_handler.jl"))
 

src/ModiaInterface/_module.jl

@@ -17,6 +17,7 @@ export Fix
 export Revolute, RevoluteWithFlange
 export Prismatic, PrismaticWithFlange
 export J123, J132, J123or132, singularRem, FreeMotion, change_rotSequenceInNextIteration!
+export Bushing, SpringDamperPtP
 
 export buildModia3D
 

src/ModiaInterface/model3D.jl

@@ -21,12 +21,13 @@ Font(            ; kwargs...) = Par(; _constructor = :(Modia3D.Font)
 TextShape(       ; kwargs...) = Par(; _constructor = :(Modia3D.TextShape)             , kwargs...)
 ModelicaShape(   ; kwargs...) = Par(; _constructor = :(Modia3D.ModelicaShape)         , kwargs...)
 Fix(             ; kwargs...) = Par(; _constructor = :(Modia3D.Fix)                   , kwargs...)
+Bushing(         ; kwargs...) = Par(; _constructor = :(Modia3D.Bushing)               , kwargs...)
+SpringDamperPtP( ; kwargs...) = Par(; _constructor = :(Modia3D.SpringDamperPtP)       , kwargs...)
 
-MassPropertiesFromShape()  = Par(; _constructor = :(Modia3D.MassPropertiesFromShape))
-MassPropertiesFromShapeAndMass(;mass) = Par(; _constructor = :(Modia3D.MassPropertiesFromShapeAndMass), mass = mass)
+MassPropertiesFromShape() = Par(; _constructor = :(Modia3D.MassPropertiesFromShape))
+MassPropertiesFromShapeAndMass(; mass) = Par(; _constructor = :(Modia3D.MassPropertiesFromShapeAndMass), mass = mass)
 UniformGravityField(; kwargs...) = Par(; _constructor = :(Modia3D.UniformGravityField), kwargs...)
 
-
 RefPath(; kwargs...) = Modia3D.ReferencePath(; kwargs...)
 
 ptpJointSpace(; kwargs...) = Modia3D.ptpJointSpace(; kwargs...)

test/Basic/BoxSpringDamperPtP.jl

@@ -0,0 +1,34 @@
+module BoxSpringDamperPtP
+
+using ModiaLang
+
+import Modia3D
+using  Modia3D.ModiaInterface
+
+SpringDamper = Model(
+    Length = 0.1,
+    Mass = 1.0,
+    Stiffness = 100.0,
+    Damping = 2.0,
+    visualMaterial = VisualMaterial(color="IndianRed1", transparency=0.5),
+    world = Object3D(feature=Scene(gravityField=UniformGravityField(g=9.81, n=[0, 0, -1]), nominalLength=:Length)),
+    worldFrame = Object3D(parent=:world,
+                          feature=Visual(shape=CoordinateSystem(length=:Length))),
+    box = Object3D(feature=Solid(shape=Box(lengthX=:Length, lengthY=:Length, lengthZ=:Length),
+                                 massProperties=MassProperties(; mass=1.0, Ixx=0.1, Iyy=0.1, Izz=0.1),
+                                 visualMaterial=:(visualMaterial))),
+    joint = FreeMotion(obj1=:world, obj2=:box, r=Var(init=[0.1, 0.2, 0.3])),
+    force = SpringDamperPtP(obj1=:world, obj2=:box, stiffness=:Stiffness, damping=:Damping)
+)
+
+springDamper = @instantiateModel(buildModia3D(SpringDamper), aliasReduction=false, unitless=true)
+
+stopTime = 6.0
+dtmax = 0.1
+requiredFinalStates = [-0.026549009773552276, -0.05309801954749593, -0.10117097852565102, 0.0144068779028959, 0.028813755813570235, 0.02595802582858583, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+simulate!(springDamper, stopTime=stopTime, dtmax=dtmax, log=true, requiredFinalStates=requiredFinalStates)
+
+@usingModiaPlot
+plot(springDamper, ["joint.r", "joint.v"], figure=1)
+
+end

test/Basic/HarmonicOscillator.jl

@@ -0,0 +1,34 @@
+module HarmonicOscillator
+
+using ModiaLang
+
+import Modia3D
+using  Modia3D.ModiaInterface
+
+Oscillator = Model(
+    Length = 0.1,
+    Mass = 1.0,
+    Stiffness = 100.0,
+    Damping = 2.0,
+    visualMaterial = VisualMaterial(color="IndianRed1", transparency=0.5),
+    world = Object3D(feature=Scene(gravityField=UniformGravityField(g=9.81, n=[0, 0, -1]), nominalLength=:Length)),
+    worldFrame = Object3D(parent=:world,
+                          feature=Visual(shape=CoordinateSystem(length=:Length))),
+    oscillator = Object3D(feature=Solid(shape=Box(lengthX=:Length, lengthY=:Length, lengthZ=:Length),
+                                        massProperties=MassProperties(; mass=1.0, Ixx=0.1, Iyy=0.1, Izz=0.1),
+                                        visualMaterial=:(visualMaterial))),
+    joint = Prismatic(obj1=:world, obj2=:oscillator, axis=3, s=Var(init=0.0), v=Var(init=0.0)),
+    force = Bushing(obj1=:world, obj2=:oscillator, stiffness=:[0.0, 0.0, Stiffness], damping=:[0.0, 0.0, Damping])
+)
+
+oscillator = @instantiateModel(buildModia3D(Oscillator), aliasReduction=false, unitless=true)
+
+stopTime = 5.0
+dtmax = 0.1
+requiredFinalStates = [0.0032763240987720987, -0.09755848695453585]
+simulate!(oscillator, stopTime=stopTime, dtmax=dtmax, log=true, requiredFinalStates=requiredFinalStates)
+
+@usingModiaPlot
+plot(oscillator, ["joint.s", "joint.v"], figure=1)
+
+end

test/includeTests.jl

@@ -9,6 +9,8 @@ Test.@testset "Basic" begin
     include(joinpath("Basic", "PendulumWithParameterizedDamper.jl"))
     include(joinpath("Basic", "PendulumWithSpring.jl"))
     include(joinpath("Basic", "DoublePendulumWithDampers.jl"))
+    include(joinpath("Basic", "HarmonicOscillator.jl"))
+    include(joinpath("Basic", "BoxSpringDamperPtP.jl"))
     include(joinpath("Basic", "BoxPlanarMotion.jl"))
     include(joinpath("Basic", "ShaftFreeMotion.jl"))
     include(joinpath("Basic", "ShaftFreeMotionAdaptiveRotSequence.jl"))