Merge pull request #94 from ModiaSim/an_elasticResponse

An elastic response

Author: AndreaNeumayr

Project.toml

@@ -1,7 +1,7 @@
 name = "Modia3D"
 uuid = "07f2c1e0-90b0-56cf-bda7-b44b56e34eed"
 authors = ["Andrea Neumayr <[email protected]>", "Martin Otter <[email protected]>", "Gerhard Hippmann <[email protected]>"]
-version = "0.8.0"
+version = "0.8.1"
 
 [deps]
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"

docs/src/index.md

@@ -48,12 +48,21 @@ Download and install the free DLR SimVis Community Edition, e.g. with [https://v
 
 ## Release Notes
 
+### Version 0.8.1
+
+- `collisionSmoothingRadius`: edges are smoothed with `collisionSmoothingRadius`
+    - its default value is `0.001`, if it is set to `0.0` no smoothing takes place
+    - the edges of Box, Cylinder, Cone, and Beam are smoothened with `collisionSmoothingRadius`
+- For using Herz' pressure `contactSphereRadius` is introduced for each shape
+    - it can be user set
+    - otherwise, it is calculated from shape geometry
+
 ### Version 0.8.0
 
 - Require ModiaLang 0.10.0 (e.g. require DifferentialEquations 7)
-- Support animation export for FloatType = Measurements.XXX (animation is performed for the nominal value).
-- Remove PathPlanning struct und functions (and instead use the functionality from ModiaLang 0.9.1). 
-- Reactivate test model test/old/Test_PathPlanning.jl
+- Support animation export for `FloatType = Measurements.XXX` (animation is performed for the nominal value).
+- Remove PathPlanning struct und functions (and instead use the functionality from ModiaLang 0.9.1).
+- Reactivate test model `test/old/Test_PathPlanning.jl`
 - Remove unused files: .codecov.yml, .travis.yml, appveyor.yml
 - Fix import/using of some tests.
 

src/Composition/responseCalculation/elasticCollisionResponse.jl

@@ -36,28 +36,24 @@ function resultantDampingCoefficient(cor, abs_vreln, vsmall, maximumContactDampi
     return d_res
 end
 
-
 function elasticContactPairCoefficients(obj1::Object3D, obj2::Object3D)
-  if typeof(obj1.feature.shape) <: Modia3D.Shapes.Sphere && typeof(obj2.feature.shape) <: Modia3D.Shapes.Sphere
-    r1 = obj1.feature.shape.diameter*0.5
-    r2 = obj2.feature.shape.diameter*0.5
-    mu_r_geo = r1*r2/(r1 + r2)
-    n_geo = 1.5
-    c_geo = 4/3*sqrt(mu_r_geo)
-  elseif typeof(obj1.feature.shape) <: Modia3D.Shapes.Sphere && typeof(obj2.feature.shape) != Modia3D.Shapes.Sphere
-    mu_r_geo = obj1.feature.shape.diameter*0.5
-    n_geo = 1.5
-    c_geo = 4/3*sqrt(mu_r_geo)
-  elseif typeof(obj1.feature.shape) != Modia3D.Shapes.Sphere && typeof(obj2.feature.shape) <: Modia3D.Shapes.Sphere
-    mu_r_geo = obj2.feature.shape.diameter*0.5
+    solid1::Shapes.Solid = obj1.feature
+    solid2::Shapes.Solid = obj2.feature
+
+    if !solid1.isFlat && solid2.isFlat
+        mu_r_geo = solid1.contactSphereRadius
+    elseif solid1.isFlat && !solid2.isFlat
+        mu_r_geo = solid2.contactSphereRadius
+    else # (solid1.isFlat && solid2.isFlat) || (!solid1.isFlat && !solid2.isFlat)
+        r1 = solid1.contactSphereRadius
+        r2 = solid2.contactSphereRadius
+        mu_r_geo = r1*r2/(r1 + r2)
+    end
+
     n_geo = 1.5
     c_geo = 4/3*sqrt(mu_r_geo)
-  else
-    mu_r_geo = 1.0
-    n_geo = 1.0
-    c_geo = 1.0
-  end
-  return (c_geo, n_geo, mu_r_geo)
+
+    return (c_geo, n_geo, mu_r_geo)
 end
 
 
@@ -86,28 +82,28 @@ function contactStart(matPair::Shapes.ElasticContactPairMaterial,
     nu1   = mat1.PoissonsRatio
     nu2   = mat2.PoissonsRatio
     if E1 <= 0.0 || E2 <= 0.0 || nu1 <= 0.0 || nu1 >= 1.0 ||
-       nu2 <= 0.0 || nu2 >= 1.0 || elasticContactReductionFactor <= 0.0
-      responseMaterial = nothing
+        nu2 <= 0.0 || nu2 >= 1.0 || elasticContactReductionFactor <= 0.0
+        responseMaterial = nothing
     else
-      @assert(E1 > 0.0)
-      @assert(E2 > 0.0)
-      @assert(nu1 > 0.0 && nu1 < 1.0)
-      @assert(nu2 > 0.0 && nu2 < 1.0)
-      @assert(elasticContactReductionFactor > 0.0)
-      c1 = E1/(1 - nu1^2)
-      c2 = E2/(1 - nu2^2)
-      c_res = elasticContactReductionFactor*c1*c2/(c1 + c2)
-
-      # Compute damping constant
-      delta_dot_start = normalRelativeVelocityAtContact(obj1, obj2, rContact, contactNormal)
-      d_res = Modia3D.resultantDampingCoefficient(matPair.coefficientOfRestitution, abs(delta_dot_start), matPair.vsmall, maximumContactDamping)
-
-      # Determine other coefficients
-      (c_geo, n_geo, mu_r_geo) = elasticContactPairCoefficients(obj1,obj2)
-      responseMaterial = ElasticContactPairResponseMaterial(c_res, c_geo, n_geo, d_res,
-                              matPair.slidingFrictionCoefficient,
-                              matPair.rotationalResistanceCoefficient, mu_r_geo,
-                              matPair.vsmall, matPair.wsmall)
+        @assert(E1 > 0.0)
+        @assert(E2 > 0.0)
+        @assert(nu1 > 0.0 && nu1 < 1.0)
+        @assert(nu2 > 0.0 && nu2 < 1.0)
+        @assert(elasticContactReductionFactor > 0.0)
+        c1 = E1/(1 - nu1^2)
+        c2 = E2/(1 - nu2^2)
+        c_res = elasticContactReductionFactor*c1*c2/(c1 + c2)
+
+        # Compute damping constant
+        delta_dot_start = normalRelativeVelocityAtContact(obj1, obj2, rContact, contactNormal)
+        d_res = Modia3D.resultantDampingCoefficient(matPair.coefficientOfRestitution, abs(delta_dot_start), matPair.vsmall, maximumContactDamping)
+
+        # Determine other coefficients
+        (c_geo, n_geo, mu_r_geo) = elasticContactPairCoefficients(obj1,obj2)
+        responseMaterial = ElasticContactPairResponseMaterial(c_res, c_geo, n_geo, d_res,
+                                matPair.slidingFrictionCoefficient,
+                                matPair.rotationalResistanceCoefficient, mu_r_geo,
+                                matPair.vsmall, matPair.wsmall)
     end
     return responseMaterial
 end
@@ -163,11 +159,8 @@ function responseCalculation(material::ElasticContactPairResponseMaterial, obj1:
     e_t_reg   = v_t/Modia3D.regularize(norm(v_t),vsmall)
     delta     = -s
 
-    if c_geo != 1.0 && n_geo != 1.0
-      delta_comp = delta * sqrt(abs(delta))
-    else
-      delta_comp = delta
-    end
+    delta_comp = delta * sqrt(abs(delta))
+
     #fn = -max(F(0.0), c_res * c_geo * delta_comp * (1 - d_res*delta_dot) )
     fn = -c_res * c_geo * delta_comp * (1 - d_res*delta_dot)
     ft = -mu_k*fn*e_t_reg

src/Modia3D.jl

@@ -4,8 +4,8 @@
 module Modia3D
 
 const path = dirname(dirname(@__FILE__))   # Absolute path of package directory
-const Version = "0.8.0"
-const Date = "2022-02-07"
+const Version = "0.8.1"
+const Date = "2022-02-11"
 
 
 # println("\nImporting Modia3D Version $Version ($Date)")

src/Shapes/geometry.jl

@@ -235,6 +235,7 @@ mutable struct FileMesh <: Modia3D.AbstractGeometry
 
     # for solids only
     centroid::SVector{3,Float64}
+    shortestEdge::Float64
     longestEdge::Float64
     objPoints::Vector{SVector{3,Float64}}
     facesIndizes::Vector{SVector{3,Int64}}

src/Shapes/setCollisionSmoothingRadius.jl

@@ -30,3 +30,37 @@ function setCollisionSmoothingRadius(shape::Beam{F}, collisionSmoothingRadius::F
     collisionSmoothingRadius2 = min(collisionSmoothingRadius, F(0.1)*min(shape.length, shape.width, shape.thickness)) # at most 10% of the smallest edge length
     return collisionSmoothingRadius2
 end
+
+
+setContactSphereRadius(shape::Nothing, contactSphereRadius::Nothing, ::Type{F}) where F <: Modia3D.VarFloatType = (false, F(NaN))
+
+setContactSphereRadius(shape::Sphere{F}, contactSphereRadius::Nothing, ::Type{F}) where F <: Modia3D.VarFloatType = (false, F(shape.diameter * 0.5) )
+setContactSphereRadius(shape::Sphere{F}, contactSphereRadius::F, ::Type{F}) where F <: Modia3D.VarFloatType = (false, contactSphereRadius)
+
+
+setContactSphereRadius(shape::Ellipsoid{F}, contactSphereRadius::Nothing, ::Type{F}) where F <: Modia3D.VarFloatType = (false, F(min(shape.lengthX, shape.lengthY, shape.lengthZ)*0.5) )
+setContactSphereRadius(shape::Ellipsoid{F}, contactSphereRadius::F, ::Type{F}) where F <: Modia3D.VarFloatType = (false, contactSphereRadius )
+
+
+setContactSphereRadius(shape::Cylinder{F}, contactSphereRadius::Nothing, ::Type{F}) where F <: Modia3D.VarFloatType = (false, F(min(shape.diameter, shape.length)*0.5) )
+setContactSphereRadius(shape::Cylinder{F}, contactSphereRadius::F, ::Type{F}) where F <: Modia3D.VarFloatType = (false, contactSphereRadius )
+
+
+setContactSphereRadius(shape::Capsule{F}, contactSphereRadius::Nothing, ::Type{F}) where F <: Modia3D.VarFloatType = (false, F(shape.diameter*0.5) )
+setContactSphereRadius(shape::Capsule{F}, contactSphereRadius::F, ::Type{F}) where F <: Modia3D.VarFloatType = (false, contactSphereRadius )
+
+
+setContactSphereRadius(shape::Cone{F}, contactSphereRadius::Nothing, ::Type{F}) where F <: Modia3D.VarFloatType = (false, F((shape.diameter + shape.topDiameter)*0.25) )
+setContactSphereRadius(shape::Cone{F}, contactSphereRadius::F, ::Type{F}) where F <: Modia3D.VarFloatType = (false, contactSphereRadius )
+
+
+setContactSphereRadius(shape::Box{F}, contactSphereRadius::Nothing, ::Type{F}) where F <: Modia3D.VarFloatType = (true, F(min(shape.lengthX, shape.lengthY, shape.lengthZ)*0.5) )
+setContactSphereRadius(shape::Box{F}, contactSphereRadius::F, ::Type{F}) where F <: Modia3D.VarFloatType = (true, contactSphereRadius )
+
+
+setContactSphereRadius(shape::Beam{F}, contactSphereRadius::Nothing, ::Type{F}) where F <: Modia3D.VarFloatType = (true, F(min(shape.length, shape.width, shape.thickness)*0.5) )
+setContactSphereRadius(shape::Beam{F}, contactSphereRadius::F, ::Type{F}) where F <: Modia3D.VarFloatType = (true, contactSphereRadius )
+
+
+setContactSphereRadius(shape::FileMesh, contactSphereRadius::Nothing, ::Type{F}) where F <: Modia3D.VarFloatType = (false, F(shape.shortestEdge * 0.5) )
+setContactSphereRadius(shape::FileMesh, contactSphereRadius::F, ::Type{F}) where F <: Modia3D.VarFloatType = (false, contactSphereRadius )

src/Shapes/solid.jl

@@ -4,7 +4,8 @@
             massProperties = nothing,
             collision = false,
             contactMaterial = "",
-            collisionSmoothingRadius = 0.0,
+            collisionSmoothingRadius = 0.001,
+            contactSphereRadius = nothing,
             visualMaterial = VisualMaterial())
 
 Generate a [Solid](@ref) with physical behavior of a rigid body with mass, visualization and collision properties.
@@ -35,6 +36,22 @@ Generate a [Solid](@ref) with physical behavior of a rigid body with mass, visua
 
 - `collisionSmoothingRadius`: Defines a collision smoothing radius for surface edges, its default value is `0.001`. It takes the minimum value of your collision smoothing radius and 10% of the smallest shape length, like `min(collisionSmoothingRadius, 0.1 min(shape dimensions))`. If it is set to `0.0` no `collisionSmoothingRadius` is used. A `collisionSmoothingRadius` is used for `Box`, `Cylinder`, `Cone`, and `Beam`.
 
+- `contactSphereRadius`: for each shape a `contactSphereRadius` is defined. So that Herz' pressure is used in [Response calculation](@ref) not only for spheres. You can define your own `contactSphereRadius`, otherwise it is computed from shape geometry (sketched in the following table).
+
+
+| Shape                     | `contactSphereRadius` from shape |
+|:--------------------------|:---------------------------------|
+| [Sphere](@ref)            | `diameter/2` |
+| [Ellipsoid](@ref)         | `min(lengthX, lengthY, lengthZ)/2` |
+| [Box](@ref)               | `min(lengthX, lengthY, lengthZ)/2` |
+| [Cylinder](@ref)          | `min(diameter, length)/2` |
+| [Cone](@ref)              | `(diameter + topDiameter)/4` |
+| [Capsule](@ref)           | `diameter/2` |
+| [Beam](@ref)              | `min(length, width, thickness)/2` |
+| [FileMesh](@ref)          | `shortestEdge/2` |
+
+For flat shapes, [Box](@ref) and [Beam](@ref), no `contactSphereRadius` is taken.  For Herz' pressure it is needed only if two flat shapes are colliding.
+
 - `visualMaterial`: Defines the material of the solid used for visualization. A pre-defined [Visual material](@ref)
     from palettes/visualMaterials.json (e.g. `"RedTransparent"`) or a user-defined [Visual material](@ref) (e.g.
     `VisualMaterial(color="DeepSkyBlue4", transparency=0.75)`) can be used.
@@ -59,6 +76,8 @@ struct Solid{F <: Modia3D.VarFloatType} <: Modia3D.AbstractObject3DFeature
     contactMaterial::Union{String,Modia3D.AbstractContactMaterial,Nothing}
     collisionSmoothingRadius::F
     visualMaterial::Union{Shapes.VisualMaterial,Nothing}
+    isFlat::Bool
+    contactSphereRadius::F
 
     function Solid{F}(;
         shape::Union{Modia3D.AbstractGeometry,Nothing} = nothing,
@@ -68,7 +87,8 @@ struct Solid{F <: Modia3D.VarFloatType} <: Modia3D.AbstractObject3DFeature
         contactMaterial::AbstractString = "",
         collisionSmoothingRadius=F(0.001),
         visualMaterial::Union{Shapes.VisualMaterial,AbstractString,Nothing} = Shapes.VisualMaterial(),
-        visualMaterialConvexDecomposition::Union{Shapes.VisualMaterial,AbstractString,Nothing} = Shapes.VisualMaterial()) where F <: Modia3D.VarFloatType
+        visualMaterialConvexDecomposition::Union{Shapes.VisualMaterial,AbstractString,Nothing} = Shapes.VisualMaterial(),
+        contactSphereRadius::Union{Nothing, F} = nothing) where F <: Modia3D.VarFloatType
 
         if collision && isnothing(shape)
             error("For collision/gripping simulations, a shape must be defined.")
@@ -94,12 +114,13 @@ struct Solid{F <: Modia3D.VarFloatType} <: Modia3D.AbstractObject3DFeature
         end
 
         if typeof(shape) == FileMesh
-            (shape.centroid, shape.longestEdge, shape.objPoints, shape.facesIndizes) = getMeshInfos(shape.filename, shape.scaleFactor)
+            (shape.centroid, shape.shortestEdge, shape.longestEdge, shape.objPoints, shape.facesIndizes) = getMeshInfos(shape.filename, shape.scaleFactor)
             (shape.volume, shape.centroidAlgo, shape.inertia) = computeMassProperties(shape.objPoints, shape.facesIndizes; bodyCoords=false)
         end
 
         massProperties = createMassProperties(F, massProperties, shape, solidMaterial)
-        new(shape, solidMaterial, massProperties, collision, contactMaterial, setCollisionSmoothingRadius(shape, F(collisionSmoothingRadius)), visualMaterial)
+        (isFlat, contactSphereRadius) = setContactSphereRadius(shape, contactSphereRadius, F)
+        new(shape, solidMaterial, massProperties, collision, contactMaterial, setCollisionSmoothingRadius(shape, F(collisionSmoothingRadius)), visualMaterial, isFlat, contactSphereRadius)
     end
 end
 Solid(args...; kwargs...) = Solid{Float64}(args...; kwargs...)

src/Shapes/utilities.jl

@@ -36,7 +36,7 @@ end
 
 
 # Utility functions for FileMesh.jl
-# getMeshInfos(): read file and return centroid, longestEdge and mesh data
+# getMeshInfos(): read file and return centroid, shortestEdge, longestEdge and mesh data
 #
 # VHACD(): decomposes a concave file mesh in convex sub parts
 #          see: https://github.com/kmammou/v-hacd
@@ -47,6 +47,7 @@ function getMeshInfos(filename::AbstractString, scaleFactor::SVector{3,Float64})
     vertices = Vector{SVector{3,Float64}}()
     faces = Vector{SVector{3,Int64}}()
     centroid::SVector{3,Float64} = Modia3D.ZeroVector3D(Float64)
+    shortestEdge::Float64 = 0.0
     longestEdge::Float64 = 0.0
 
     mesh = FileIO.load(filename, pointtype=MeshIO.Point3{Float64}, facetype=MeshIO.TriangleFace{MeshIO.OneIndex{Int64}})
@@ -86,10 +87,11 @@ function getMeshInfos(filename::AbstractString, scaleFactor::SVector{3,Float64})
             sum += vertices[i]
         end
         centroid = sum / length(vertices)
+        shortestEdge = min((x_max-x_min), (y_max-y_min), (z_max-z_min))
         longestEdge = max((x_max-x_min), (y_max-y_min), (z_max-z_min))
     end
 
-    return (centroid, longestEdge, vertices, faces)
+    return (centroid, shortestEdge, longestEdge, vertices, faces)
 end
 
 getMeshInfos(filename::AbstractString, scaleFactor::AbstractVector) =

test/Collision/BouncingBeams.jl

@@ -77,7 +77,8 @@ bouncingBeams = @instantiateModel(buildModia3D(BouncingBeams), unitless=true, lo
 
 stopTime = 1.2
 tolerance = 1e-8
-requiredFinalStates = [-0.8527742342910635, 1.6578017942923853, 0.04577174265530922, -0.7904903307190998, 0.3784209895248903, -0.048541407443668916, 4.511726495865782, 1.371932459718465, -1.510024797043685, -7.155505879183188, -0.6205909597968676, 0.6202303171080142, 0.04600259422130081, -0.8603712066366127, 1.6591872146181508, -0.05012483182230522, -0.8502027335229523, 0.39228915465700764, 3.3084796917299557, 0.1357058823264399, -1.6011728205254234, 0.5652834075330013, -7.593556822874329, -0.55350211235096, 1.6580899516035195, 0.04573950567885933, -0.8547165339765472, 0.38182577225339254, -0.048923116390466426, -0.8057409130731289, -1.559903565016631, 0.140455404024548, 3.3294946971080854, -0.6029565330774422, 0.6048964881953146, -7.269193867332278]
+requiredFinalStates = [-0.8104568745056254, 1.7973874611151963,
+-0.04182468212042101, -0.6843977605721305, 1.1626826154110008, -0.4579740595633442, 4.49795019668019, -1.2036633128916912, -1.4661864515779968, 7.374259609209814, 0.25467343019105243, 0.09496510178329644, -0.04316792788452747, -0.8118074155564607, 1.79878733239027, -0.45837679978147955, -0.7109008517331268, 1.1662784938747761, 5.9164930138891965, -0.3057751186242098, -1.6453594624487566, 0.09509462328382953, 7.372754384939963, 0.2531821792141959, 1.8050541127713922, -0.04525571158341206, -0.8223444493426333, 1.1809620524502702, -0.4635597953934979, -0.8311350825736231, -1.527376219825907, 0.31839879289639594, 5.985690312517948, 0.3021665793506058, 0.14598734167247776, 7.745079948909601]
 simulate!(bouncingBeams, stopTime=stopTime, tolerance=tolerance, log=true, logStates=false, logEvents=false, requiredFinalStates=requiredFinalStates)
 
 @usingModiaPlot

test/Collision/BouncingCapsules.jl

@@ -66,9 +66,9 @@ bouncingCapsules = @instantiateModel(buildModia3D(BouncingCapsules), unitless=tr
 #@show bouncingCapsules.parameterExpressions
 #@show bouncingCapsules.parameters
 
-stopTime = 1.7
+stopTime = 1.4
 tolerance = 1e-8
-requiredFinalStates = [-0.7965194324923187, 1.0792671324200518, 0.22983676627347815, -0.1039097850692818, 0.4814176273397363, -0.07448832005070478, 2.944368266298714, -0.23188092804556582, 1.5637806069044045, 2.4254147963752337, 0.1958593686673586, -0.049049991414040636, 0.22982554355367993, -0.7965131155676861, 1.0792581962112746, -0.07973508965627857, -0.11418691979366435, 0.48039300578570787, 3.377894872964063, 0.19183612343753909, -1.6238015814380535, -0.0420638591264601, 2.4255988934027917, 0.1666656948426664, 1.0793274621567706, 0.23012495483972084, -0.7967193671718964, 0.4805756438722489, -0.07905877747249435, -0.12952010214850182, -1.5655006683501074, 0.19873974775783254, 3.3726826053720087, 0.17806826424566943, -0.04410873189732083, 2.426030228460745]
+requiredFinalStates = [-0.7673604039864523, 0.9516988414284928, 0.24982620243724654, 0.014180822268884124, 0.5036323273267752, -0.18340022474242668, 1.5220127656990843, 1.1890544892808237, 1.4103120626088645, 2.622092435351728, -0.22352345880886096, -0.10544477069267612, 0.2498438864904243, -0.7676017038880056, 0.951656101955695, -0.18458496829981588, 0.013038547274916143, 0.503371330630657, 2.7566146861930125, 0.20558772188486563, -1.552314769388157, -0.10825787299295013, 2.6224532099355535, -0.23084017453690012, 0.9516460970951072, 0.249610396714428, -0.767900681642965, 0.5031767718168185, -0.1856868932825187, 0.017287282775231304, -1.510292025990061, 0.19727749473740247, 2.753891529551364, -0.23790362799901707, -0.11054266322392534, 2.6231546726843544]
 simulate!(bouncingCapsules, stopTime=stopTime, tolerance=tolerance, log=true, logStates=false, logEvents=false, requiredFinalStates=requiredFinalStates)
 
 @usingModiaPlot