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 3bf38910..7ea5799b 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, @@ -148,15 +148,58 @@ 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 [Shapes](@ref)). +If a user specific `contactSphereRadius` is defined in [Solid](@ref), it is taken instead of the default value. + +The default $r_{contact}$ 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 | + +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}$ | +|:-------|:-----------------|:------------------------| +|false | false | $\frac{r1 r2}{r1 + r2}$ | +|false | true | $r1$ | +|true | false | $r2$ | +|true | true | $\frac{r1 r2}{r1 + r2}$ | + + +$n_{geo} = 1.5$ +$c_{geo} = \frac{4}{3} \sqrt(\mu_{r,geo})$ -| 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 @@ -219,7 +262,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 @@ -262,3 +305,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. 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()