From a696c391e65bb932693c4b082ae066f9d6e9b87b Mon Sep 17 00:00:00 2001 From: "Martin.Otter@dlr.de" Date: Fri, 25 Feb 2022 01:25:38 +0100 Subject: [PATCH 1/5] Improve docu for contact force law --- docs/src/internal/ContactForceLaw.md | 45 +++++++++++++++++++++++----- 1 file changed, 37 insertions(+), 8 deletions(-) diff --git a/docs/src/internal/ContactForceLaw.md b/docs/src/internal/ContactForceLaw.md index 3bf38910..9b0da19e 100644 --- a/docs/src/internal/ContactForceLaw.md +++ b/docs/src/internal/ContactForceLaw.md @@ -148,15 +148,33 @@ In special cases (for example sphere rolling on a plane), the rotational coeffic can be interpreted as *rolling resistance coefficient*. Coefficients ``c_{geo}, n_{geo}, \mu_{r,geo}`` depend on the geometries of the objects -that are in contact. Only for spheres meaning values are provided based on Hertz' pressure, -because currently the collision handling in Modia3D does no provide enough information for other -geometries (``r_i`` is the radius of sphere ``i``): +that are in contact. The coefficients are computed approximately based on the contact theory +of Hertz [^5], [^6]: Here, it is assumed that each of the contacting surfaces can be described by a +quadratic polynomial in two variables that is basically defined by its principal curvatures +along two perpendicular directions at the point of contact. A characteristic feature is that +the contact volume increases nonlinearly with the penetration depth, so ``n_{geo} > 1`` (provided +the two contacting surfaces are not completely flat), and therefore the normal contact force +changes nonlinearly with the penetration depth. In the general case, elliptical integrals +have to be solved, as well as a nonlinear algebraic equation system to compute the normal +contact force as function of the penetration depth and the principal curvatures at the contact point. +An approximate *analytical* model is proposed in [^7]. + +In order that a numerical integration algorithm with step-size control +works reasonably, the contact force needs to be continuous and continuously differentiable with +respect to the penetration depth. This in turn means that the principal curvatures of the contacting +surfaces should also be continuous and continuously differentiable, which is usually not the case +(besides exceptional cases, such as a Sphere or an Ellipsoid). + +Since the determination of the principal curvatures of shapes is in general +complicated and the shapes have often areas with discontinuous curvatures, only a very rough approximation +is used in Modia3D: *The contact area of a shape is approximated by a quadratic polynomial +with constant mean principal curvature in all directions and on all points on the shape*. +In other words, a sphere with constant sphere radius ``r_{contact}`` is associated with every shape that +is used to compute coefficients ``c_{geo}, n_{geo}, \mu_{r,geo}``. A default value for ``r_{contact}`` +is determined based on the available data of the shape (see [shape data](https://modiasim.github.io/Modia3D.jl/stable/Components/Shapes.html)): + +xxx -| Object 1 | Object 2 | ``c_{geo}`` | ``n_{geo}`` | ``\mu_{r,geo}`` | -|:----------- |:--------- |:---------------------------------------- |:------------|:------------------- | -| Sphere | Sphere | ``\frac{4}{3} \sqrt{1/(1/r_1+1/r_2)}`` | ``1.5`` | ``1/(1/r_1+1/r_2)`` | -| Sphere | no Sphere | ``\frac{4}{3} \sqrt{r_1}`` | ``1.5`` | ``r_1`` | -| no Sphere | no Sphere | ``1`` | ``1.0`` | ``1.0`` | ## Regularized unit vectors @@ -262,3 +280,14 @@ similar responses: [^4]: Andrea Neumayr, Martin Otter (2019): [Collision Handling with Elastic Response Calculation and Zero-Crossing Functions](https://doi.org/10.1145/3365984.3365986). Proceedings of the 9th International Workshop on Equation-Based Object-Oriented Modeling Languages and Tools. EOOLT’19. ACM, pp. 57–65. + +[^5]: Hertz H. (1881): + [Über die Berührung fester elastischer Körper](https://home.uni-leipzig.de/pwm/web/download/Hertz1881.pdf). + Journal für die reine und angewandte Mathematik 92, S. 156-171. + +[^6]: Johnson K.L. (1985): + Contact Mechanics. Cambridge University Press. + +[^7]: Antoine J-F., Visa C., and Sauvey C. (2006): + [Approximate Analytical Model for Hertzian Elliptical Contact Problems](https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.1055.4455&rep=rep1&type=pdf). + Transactions of the ASME, Vol. 128. pp. 660-664. \ No newline at end of file From 195a84e6936615df49e8048c517f56dd0c1afbe5 Mon Sep 17 00:00:00 2001 From: "Martin.Otter@dlr.de" Date: Fri, 25 Feb 2022 08:22:17 +0100 Subject: [PATCH 2/5] Changed maximumContactDamping from 1000 to 2000 (text+figure), since this value used in Modia3D --- docs/resources/images/plot_damping1.svg | 1009 +++------ docs/resources/images/plot_damping2.svg | 2502 ----------------------- docs/src/internal/ContactForceLaw.md | 2 +- test/old/Plot_cor.jl | 46 +- 4 files changed, 316 insertions(+), 3243 deletions(-) delete mode 100644 docs/resources/images/plot_damping2.svg diff --git a/docs/resources/images/plot_damping1.svg b/docs/resources/images/plot_damping1.svg index 7fc254cb..03a79b4a 100644 --- a/docs/resources/images/plot_damping1.svg +++ b/docs/resources/images/plot_damping1.svg @@ -2,12 +2,12 @@ - + - 2022-01-19T17:53:14.793151 + 2022-02-25T08:10:47.721014 image/svg+xml @@ -22,42 +22,42 @@ - - - +" id="m50454c9a4c" style="stroke:#000000;stroke-width:0.8;"/> - + - + - - + - + - - + - + @@ -180,18 +180,18 @@ L 162.946987 38.5344 - - + - + - - + - + @@ -250,18 +250,18 @@ L 253.676338 38.5344 - - + - + - - + - + @@ -328,18 +328,18 @@ L 344.405688 38.5344 - - + - + - + - +" id="mb340fae211" style="stroke:#000000;stroke-width:0.8;"/> - + - + - - + - - + + - + - - + - - - + + + @@ -584,141 +584,143 @@ L 405.648 184.644 - - + - - + + - +" id="DejaVuSans-55"/> - - + + - - + - - - - - - + + + + - - + - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - + - + - + - - + - - + - - + - - + - - - - - - - + - - + + - + - - + - - + + - - - - + @@ -1540,15 +1119,15 @@ z - - + - - + + - + @@ -1573,15 +1152,15 @@ L 259.148 80.989087 - - + - - + + - + @@ -1606,15 +1185,15 @@ L 259.148 95.667213 - - + - - + + - + @@ -1643,8 +1222,8 @@ L 259.148 110.345337 - - + + diff --git a/docs/resources/images/plot_damping2.svg b/docs/resources/images/plot_damping2.svg deleted file mode 100644 index 4b993f98..00000000 --- a/docs/resources/images/plot_damping2.svg +++ /dev/null @@ -1,2502 +0,0 @@ - - - - - - - - image/svg+xml - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/docs/src/internal/ContactForceLaw.md b/docs/src/internal/ContactForceLaw.md index 9b0da19e..07e46b41 100644 --- a/docs/src/internal/ContactForceLaw.md +++ b/docs/src/internal/ContactForceLaw.md @@ -237,7 +237,7 @@ The damping coefficient ``d`` is basically computed with the formulation from [^ because a response calculation with impulses gives similar results for some experiments as shown in [^3]. However, (a) instead of ``cor``, the current coefficient of restitution ``cor_{cur}`` is used and (b) the damping coefficient is limited to ``d_{max}`` (this parameter is set via `maximumContactDamping` -in object `Scene` and has a default value of ``1000 ~Ns/m``) to avoid an unphysical +in object `Scene` and has a default value of ``2000 ~Ns/m``) to avoid an unphysical strong creeping effect for collisions with small ``cor_{cur}`` values: ```math diff --git a/test/old/Plot_cor.jl b/test/old/Plot_cor.jl index 08ba95d3..6223014c 100644 --- a/test/old/Plot_cor.jl +++ b/test/old/Plot_cor.jl @@ -13,10 +13,10 @@ cor2 = 0.3 cor3 = 0.1 cor4 = 0.01 -vrela = collect(range(0.5*vsmall, 4*vsmall, length=100)) -vrelb = collect(range(0 , vsmall, length=500)) -vrelc = collect(range(0 , 4*vsmall, length=100)) -vreld = collect(range(0 , 2*vsmall, length=100)) +vrela = collect(range(0.5*vsmall, 4*vsmall, length=1000)) +vrelb = collect(range(0 , vsmall, length=5000)) +vrelc = collect(range(0 , 4*vsmall, length=1000)) +vreld = collect(range(0 , 2*vsmall, length=1000)) d_res0a = zeros( length(vrela) ) d_res1a = zeros( length(vrela) ) @@ -39,36 +39,35 @@ cor_res4 = zeros( length(vrelc) ) reg = zeros( length(vreld) ) w = 0.0 for i in 1:length(vrela) - d_res0a[i] = Modia3D.resultantDampingCoefficient(cor0,vrela[i],vsmall,1000) - d_res1a[i] = Modia3D.resultantDampingCoefficient(cor1,vrela[i],vsmall,1000) - d_res2a[i] = Modia3D.resultantDampingCoefficient(cor2,vrela[i],vsmall,1000) - d_res3a[i] = Modia3D.resultantDampingCoefficient(cor3,vrela[i],vsmall,1000) - d_res4a[i] = Modia3D.resultantDampingCoefficient(cor4,vrela[i],vsmall,1000) + d_res0a[i] = Modia3D.resultantDampingCoefficient(cor0,vrela[i],vsmall,2000) + d_res1a[i] = Modia3D.resultantDampingCoefficient(cor1,vrela[i],vsmall,2000) + d_res2a[i] = Modia3D.resultantDampingCoefficient(cor2,vrela[i],vsmall,2000) + d_res3a[i] = Modia3D.resultantDampingCoefficient(cor3,vrela[i],vsmall,2000) + d_res4a[i] = Modia3D.resultantDampingCoefficient(cor4,vrela[i],vsmall,2000) end for i in 1:length(vrelb) - d_res0b[i] = Modia3D.resultantDampingCoefficient(cor0,vrelb[i],vsmall,1000) - d_res1b[i] = Modia3D.resultantDampingCoefficient(cor1,vrelb[i],vsmall,1000) - d_res2b[i] = Modia3D.resultantDampingCoefficient(cor2,vrelb[i],vsmall,1000) - d_res3b[i] = Modia3D.resultantDampingCoefficient(cor3,vrelb[i],vsmall,1000) - d_res4b[i] = Modia3D.resultantDampingCoefficient(cor4,vrelb[i],vsmall,1000) + d_res0b[i] = Modia3D.resultantDampingCoefficient(cor0,vrelb[i],vsmall,2000) + d_res1b[i] = Modia3D.resultantDampingCoefficient(cor1,vrelb[i],vsmall,2000) + d_res2b[i] = Modia3D.resultantDampingCoefficient(cor2,vrelb[i],vsmall,2000) + d_res3b[i] = Modia3D.resultantDampingCoefficient(cor3,vrelb[i],vsmall,2000) + d_res4b[i] = Modia3D.resultantDampingCoefficient(cor4,vrelb[i],vsmall,2000) end -#= +resultantCoefficientOfRestitution(cor, abs_vreln, vsmall) = abs_vreln > vsmall ? cor : 0.0 + for i in 1:length(vrelc) - cor_res0[i] = Modia3D.resultantCoefficientOfRestitution(cor0,vrelc[i],vsmall) - cor_res1[i] = Modia3D.resultantCoefficientOfRestitution(cor1,vrelc[i],vsmall) - cor_res2[i] = Modia3D.resultantCoefficientOfRestitution(cor2,vrelc[i],vsmall) - cor_res3[i] = Modia3D.resultantCoefficientOfRestitution(cor3,vrelc[i],vsmall) - cor_res4[i] = Modia3D.resultantCoefficientOfRestitution(cor4,vrelc[i],vsmall) + cor_res0[i] = resultantCoefficientOfRestitution(cor0,vrelc[i],vsmall) + cor_res1[i] = resultantCoefficientOfRestitution(cor1,vrelc[i],vsmall) + cor_res2[i] = resultantCoefficientOfRestitution(cor2,vrelc[i],vsmall) + cor_res3[i] = resultantCoefficientOfRestitution(cor3,vrelc[i],vsmall) + cor_res4[i] = resultantCoefficientOfRestitution(cor4,vrelc[i],vsmall) end -=# for i in 1:length(vreld) reg[i] = Modia3D.regularize(vreld[i],vsmall) end -#= PyPlot.figure(1) PyPlot.clf() PyPlot.plot(vrelc, cor_res0, vrelc, cor_res1, vrelc, cor_res2, vrelc, cor_res3, vrelc, cor_res4) @@ -80,7 +79,6 @@ PyPlot.legend(["\$cor = 1.0, v_{small}=0.01 \\; m/s\$", "\$cor = 0.3, v_{small}=0.01 \\; m/s\$", "\$cor = 0.1, v_{small}=0.01 \\; m/s\$", "\$cor = 0.0, v_{small}=0.01 \\; m/s\$"],loc="upper right") -=# PyPlot.figure(2) PyPlot.clf() @@ -94,7 +92,6 @@ PyPlot.legend(["\$cor = 1.0, v_{small}=0.01 \\; m/s\$", "\$cor = 0.1, v_{small}=0.01 \\; m/s\$", "\$cor = 0.0, v_{small}=0.01 \\; m/s\$"],loc="upper right") -#= PyPlot.figure(3) PyPlot.clf() PyPlot.plot(vrelb, d_res0b, vrelb, d_res1b, vrelb, d_res2b, vrelb, d_res3b, vrelb, d_res4b) @@ -106,7 +103,6 @@ PyPlot.legend(["\$cor = 1.0, v_{small}=0.01 \\; m/s\$", "\$cor = 0.3, v_{small}=0.01 \\; m/s\$", "\$cor = 0.1, v_{small}=0.01 \\; m/s\$", "\$cor = 0.0, v_{small}=0.01 \\; m/s\$"],loc="upper right") -=# PyPlot.figure(4) PyPlot.clf() From 2aab48dc872fc94bc370f4fbcd5882ae69f8603c Mon Sep 17 00:00:00 2001 From: "Martin.Otter@dlr.de" Date: Fri, 25 Feb 2022 08:30:04 +0100 Subject: [PATCH 3/5] Correct formulation --- docs/src/internal/ContactForceLaw.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/src/internal/ContactForceLaw.md b/docs/src/internal/ContactForceLaw.md index 07e46b41..35cc8915 100644 --- a/docs/src/internal/ContactForceLaw.md +++ b/docs/src/internal/ContactForceLaw.md @@ -10,7 +10,7 @@ penetration depth and the relative motion of the objects in contact. It is planned to optionally also support impulsive response calculation in the future. -This section is based on [^4] with some minor improvements as used in the current +This section is based on [^4] with some improvements as used in the current version of Modia3D. The current approach has several **limitations** that a user must know, From ed7cd7b6747c9cfd795ba204b6197912015215fe Mon Sep 17 00:00:00 2001 From: Andrea Neumayr Date: Fri, 25 Feb 2022 13:32:06 +0100 Subject: [PATCH 4/5] add table contactSphereRadius + Herz law --- docs/src/internal/ContactForceLaw.md | 47 +++++++++++++++++++++------- 1 file changed, 36 insertions(+), 11 deletions(-) diff --git a/docs/src/internal/ContactForceLaw.md b/docs/src/internal/ContactForceLaw.md index 35cc8915..7de6f2ae 100644 --- a/docs/src/internal/ContactForceLaw.md +++ b/docs/src/internal/ContactForceLaw.md @@ -150,14 +150,14 @@ can be interpreted as *rolling resistance coefficient*. Coefficients ``c_{geo}, n_{geo}, \mu_{r,geo}`` depend on the geometries of the objects that are in contact. The coefficients are computed approximately based on the contact theory of Hertz [^5], [^6]: Here, it is assumed that each of the contacting surfaces can be described by a -quadratic polynomial in two variables that is basically defined by its principal curvatures -along two perpendicular directions at the point of contact. A characteristic feature is that +quadratic polynomial in two variables that is basically defined by its principal curvatures +along two perpendicular directions at the point of contact. A characteristic feature is that the contact volume increases nonlinearly with the penetration depth, so ``n_{geo} > 1`` (provided the two contacting surfaces are not completely flat), and therefore the normal contact force changes nonlinearly with the penetration depth. In the general case, elliptical integrals have to be solved, as well as a nonlinear algebraic equation system to compute the normal -contact force as function of the penetration depth and the principal curvatures at the contact point. -An approximate *analytical* model is proposed in [^7]. +contact force as function of the penetration depth and the principal curvatures at the contact point. +An approximate *analytical* model is proposed in [^7]. In order that a numerical integration algorithm with step-size control works reasonably, the contact force needs to be continuous and continuously differentiable with @@ -168,12 +168,37 @@ surfaces should also be continuous and continuously differentiable, which is usu Since the determination of the principal curvatures of shapes is in general complicated and the shapes have often areas with discontinuous curvatures, only a very rough approximation is used in Modia3D: *The contact area of a shape is approximated by a quadratic polynomial -with constant mean principal curvature in all directions and on all points on the shape*. -In other words, a sphere with constant sphere radius ``r_{contact}`` is associated with every shape that -is used to compute coefficients ``c_{geo}, n_{geo}, \mu_{r,geo}``. A default value for ``r_{contact}`` -is determined based on the available data of the shape (see [shape data](https://modiasim.github.io/Modia3D.jl/stable/Components/Shapes.html)): +with constant mean principal curvature in all directions and on all points on the shape*. +In other words, a sphere with constant sphere radius ``r_{contact}`` is associated with every shape that is used to compute coefficients ``c_{geo}, n_{geo}, \mu_{r,geo}``. +A default value for ``r_{contact}`` is determined based on the available data of the shape (see [Shapes](@ref)): -xxx +The default values for each shape are: + +| Shape | $r_{contact}$ | isFlat | +|:----------------|:--------------------------------|:-------| +|[Sphere](@ref) | diameter/2 | false | +|[Ellipsoid](@ref)| min(lengthX, lengthY, lengthZ)/2| false | +|[Box](@ref) | min(lengthX, lengthY, lengthZ)/2| true | +|[Cylinder](@ref) | min(diameter, length)/2 | false | +|[Cone](@ref) | (diameter + topDiameter)/4 | false | +|[Capsule](@ref) | diameter/2 | false | +|[Beam](@ref) | min(length, width, thickness)/2 | true | +|[FileMesh](@ref) | shortestEdge/2 | false | + +It is possible to define a user specific `contactSphereRadius` in [Solid](@ref). +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. + + +| isFlat | isFlat | $\mu_{r,geo}$ | +|:-------|:-----------------|:------------------------| +|true | true | $\frac{r1 r2}{r1 + r2}$ | +|false | false | $\frac{r1 r2}{r1 + r2}$ | +|true | false | $r1$ | +|false | true | $r2$ | + + +$n_{geo} = 1.5$ +$c_{geo} = \frac{4}{3} \sqrt(\mu_{r,geo})$ @@ -288,6 +313,6 @@ similar responses: [^6]: Johnson K.L. (1985): Contact Mechanics. Cambridge University Press. -[^7]: Antoine J-F., Visa C., and Sauvey C. (2006): +[^7]: Antoine J-F., Visa C., and Sauvey C. (2006): [Approximate Analytical Model for Hertzian Elliptical Contact Problems](https://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.1055.4455&rep=rep1&type=pdf). - Transactions of the ASME, Vol. 128. pp. 660-664. \ No newline at end of file + Transactions of the ASME, Vol. 128. pp. 660-664. From dc2877758d6d3465cb4bceb58dc95a62fc8dee37 Mon Sep 17 00:00:00 2001 From: Andrea Neumayr Date: Mon, 28 Feb 2022 10:43:51 +0100 Subject: [PATCH 5/5] little update ContactForceLaw.md --- docs/src/internal/ContactForceLaw.md | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/docs/src/internal/ContactForceLaw.md b/docs/src/internal/ContactForceLaw.md index 7de6f2ae..7ea5799b 100644 --- a/docs/src/internal/ContactForceLaw.md +++ b/docs/src/internal/ContactForceLaw.md @@ -170,9 +170,10 @@ complicated and the shapes have often areas with discontinuous curvatures, only is used in Modia3D: *The contact area of a shape is approximated by a quadratic polynomial with constant mean principal curvature in all directions and on all points on the shape*. In other words, a sphere with constant sphere radius ``r_{contact}`` is associated with every shape that is used to compute coefficients ``c_{geo}, n_{geo}, \mu_{r,geo}``. -A default value for ``r_{contact}`` is determined based on the available data of the shape (see [Shapes](@ref)): +A default value for ``r_{contact}`` is determined based on the available data of the shape (see [Shapes](@ref)). +If a user specific `contactSphereRadius` is defined in [Solid](@ref), it is taken instead of the default value. -The default values for each shape are: +The default $r_{contact}$ values for each shape are: | Shape | $r_{contact}$ | isFlat | |:----------------|:--------------------------------|:-------| @@ -185,16 +186,15 @@ The default values for each shape are: |[Beam](@ref) | min(length, width, thickness)/2 | true | |[FileMesh](@ref) | shortestEdge/2 | false | -It is possible to define a user specific `contactSphereRadius` in [Solid](@ref). -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. +For flat shapes, [Box](@ref) and [Beam](@ref), no $r_{contact}$ is taken. For Herz' pressure it is needed only if two flat shapes are colliding ($r_i$ is the contact sphere radius $r_{contact}$ of shape $i$):). | isFlat | isFlat | $\mu_{r,geo}$ | |:-------|:-----------------|:------------------------| -|true | true | $\frac{r1 r2}{r1 + r2}$ | |false | false | $\frac{r1 r2}{r1 + r2}$ | -|true | false | $r1$ | -|false | true | $r2$ | +|false | true | $r1$ | +|true | false | $r2$ | +|true | true | $\frac{r1 r2}{r1 + r2}$ | $n_{geo} = 1.5$