@@ -36,28 +36,24 @@ function resultantDampingCoefficient(cor, abs_vreln, vsmall, maximumContactDampi
3636 return d_res
3737end
3838
39-
4039function elasticContactPairCoefficients (obj1:: Object3D , obj2:: Object3D )
41- if typeof (obj1. feature. shape) <: Modia3D.Shapes.Sphere && typeof (obj2. feature. shape) <: Modia3D.Shapes.Sphere
42- r1 = obj1. feature. shape. diameter* 0.5
43- r2 = obj2. feature. shape. diameter* 0.5
44- mu_r_geo = r1* r2/ (r1 + r2)
45- n_geo = 1.5
46- c_geo = 4 / 3 * sqrt (mu_r_geo)
47- elseif typeof (obj1. feature. shape) <: Modia3D.Shapes.Sphere && typeof (obj2. feature. shape) != Modia3D. Shapes. Sphere
48- mu_r_geo = obj1. feature. shape. diameter* 0.5
49- n_geo = 1.5
50- c_geo = 4 / 3 * sqrt (mu_r_geo)
51- elseif typeof (obj1. feature. shape) != Modia3D. Shapes. Sphere && typeof (obj2. feature. shape) <: Modia3D.Shapes.Sphere
52- mu_r_geo = obj2. feature. shape. diameter* 0.5
40+ solid1:: Shapes.Solid = obj1. feature
41+ solid2:: Shapes.Solid = obj2. feature
42+
43+ if ! solid1. isFlat && solid2. isFlat
44+ mu_r_geo = solid1. contactSphereRadius
45+ elseif solid1. isFlat && ! solid2. isFlat
46+ mu_r_geo = solid2. contactSphereRadius
47+ else # (solid1.isFlat && solid2.isFlat) || (!solid1.isFlat && !solid2.isFlat)
48+ r1 = solid1. contactSphereRadius
49+ r2 = solid2. contactSphereRadius
50+ mu_r_geo = r1* r2/ (r1 + r2)
51+ end
52+
5353 n_geo = 1.5
5454 c_geo = 4 / 3 * sqrt (mu_r_geo)
55- else
56- mu_r_geo = 1.0
57- n_geo = 1.0
58- c_geo = 1.0
59- end
60- return (c_geo, n_geo, mu_r_geo)
55+
56+ return (c_geo, n_geo, mu_r_geo)
6157end
6258
6359
@@ -86,28 +82,28 @@ function contactStart(matPair::Shapes.ElasticContactPairMaterial,
8682 nu1 = mat1. PoissonsRatio
8783 nu2 = mat2. PoissonsRatio
8884 if E1 <= 0.0 || E2 <= 0.0 || nu1 <= 0.0 || nu1 >= 1.0 ||
89- nu2 <= 0.0 || nu2 >= 1.0 || elasticContactReductionFactor <= 0.0
90- responseMaterial = nothing
85+ nu2 <= 0.0 || nu2 >= 1.0 || elasticContactReductionFactor <= 0.0
86+ responseMaterial = nothing
9187 else
92- @assert (E1 > 0.0 )
93- @assert (E2 > 0.0 )
94- @assert (nu1 > 0.0 && nu1 < 1.0 )
95- @assert (nu2 > 0.0 && nu2 < 1.0 )
96- @assert (elasticContactReductionFactor > 0.0 )
97- c1 = E1/ (1 - nu1^ 2 )
98- c2 = E2/ (1 - nu2^ 2 )
99- c_res = elasticContactReductionFactor* c1* c2/ (c1 + c2)
100-
101- # Compute damping constant
102- delta_dot_start = normalRelativeVelocityAtContact (obj1, obj2, rContact, contactNormal)
103- d_res = Modia3D. resultantDampingCoefficient (matPair. coefficientOfRestitution, abs (delta_dot_start), matPair. vsmall, maximumContactDamping)
104-
105- # Determine other coefficients
106- (c_geo, n_geo, mu_r_geo) = elasticContactPairCoefficients (obj1,obj2)
107- responseMaterial = ElasticContactPairResponseMaterial (c_res, c_geo, n_geo, d_res,
108- matPair. slidingFrictionCoefficient,
109- matPair. rotationalResistanceCoefficient, mu_r_geo,
110- matPair. vsmall, matPair. wsmall)
88+ @assert (E1 > 0.0 )
89+ @assert (E2 > 0.0 )
90+ @assert (nu1 > 0.0 && nu1 < 1.0 )
91+ @assert (nu2 > 0.0 && nu2 < 1.0 )
92+ @assert (elasticContactReductionFactor > 0.0 )
93+ c1 = E1/ (1 - nu1^ 2 )
94+ c2 = E2/ (1 - nu2^ 2 )
95+ c_res = elasticContactReductionFactor* c1* c2/ (c1 + c2)
96+
97+ # Compute damping constant
98+ delta_dot_start = normalRelativeVelocityAtContact (obj1, obj2, rContact, contactNormal)
99+ d_res = Modia3D. resultantDampingCoefficient (matPair. coefficientOfRestitution, abs (delta_dot_start), matPair. vsmall, maximumContactDamping)
100+
101+ # Determine other coefficients
102+ (c_geo, n_geo, mu_r_geo) = elasticContactPairCoefficients (obj1,obj2)
103+ responseMaterial = ElasticContactPairResponseMaterial (c_res, c_geo, n_geo, d_res,
104+ matPair. slidingFrictionCoefficient,
105+ matPair. rotationalResistanceCoefficient, mu_r_geo,
106+ matPair. vsmall, matPair. wsmall)
111107 end
112108 return responseMaterial
113109end
@@ -163,11 +159,8 @@ function responseCalculation(material::ElasticContactPairResponseMaterial, obj1:
163159 e_t_reg = v_t/ Modia3D. regularize (norm (v_t),vsmall)
164160 delta = - s
165161
166- if c_geo != 1.0 && n_geo != 1.0
167- delta_comp = delta * sqrt (abs (delta))
168- else
169- delta_comp = delta
170- end
162+ delta_comp = delta * sqrt (abs (delta))
163+
171164 # fn = -max(F(0.0), c_res * c_geo * delta_comp * (1 - d_res*delta_dot) )
172165 fn = - c_res * c_geo * delta_comp * (1 - d_res* delta_dot)
173166 ft = - mu_k* fn* e_t_reg
0 commit comments