Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F79247804
test_hertz_disp.py
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sun, Aug 25, 11:44
Size
4 KB
Mime Type
text/x-python
Expires
Tue, Aug 27, 11:44 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
20167218
Attached To
rTAMAAS tamaas
test_hertz_disp.py
View Options
# -*- coding: utf-8 -*-
# @file
# @section LICENSE
#
# Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne),
# Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program. If not, see <https://www.gnu.org/licenses/>.
from
__future__
import
print_function
,
division
import
tamaas
as
tm
import
numpy
as
np
def
constructHertzProfile
(
size
,
curvature
):
radius
=
1.
/
curvature
x
=
np
.
linspace
(
-
0.5
,
0.5
,
size
,
dtype
=
tm
.
dtype
)
y
=
np
.
linspace
(
-
0.5
,
0.5
,
size
,
dtype
=
tm
.
dtype
)
x
,
y
=
np
.
meshgrid
(
x
,
y
)
surface
=
np
.
sqrt
(
radius
**
2
-
x
**
2
-
y
**
2
)
surface
-=
surface
.
min
()
return
surface
.
copy
()
def
computeHertzDisplacement
(
e_star
,
contact_size
,
max_pressure
,
size
):
x
=
np
.
linspace
(
-
0.5
,
0.5
,
size
)
y
=
np
.
linspace
(
-
0.5
,
0.5
,
size
)
x
,
y
=
np
.
meshgrid
(
x
,
y
)
disp
=
np
.
pi
*
max_pressure
/
(
4
*
contact_size
*
e_star
)
\
*
(
2
*
contact_size
**
2
-
(
x
**
2
+
y
**
2
))
return
disp
.
copy
()
def
test_hertz_disp
():
tm
.
initialize
()
grid_size
=
512
curvature
=
0.5
effective_modulus
=
1.
load
=
0.12
surface
=
constructHertzProfile
(
grid_size
,
curvature
)
model
=
tm
.
ModelFactory
.
createModel
(
tm
.
model_type_basic_2d
,
[
1.
,
1.
],
[
grid_size
,
grid_size
])
model
.
setElasticity
(
1
,
0
)
solver
=
tm
.
PolonskyKeerRey
(
model
,
surface
,
1e-12
,
tm
.
PolonskyKeerRey
.
gap
,
tm
.
PolonskyKeerRey
.
gap
)
solver
.
solve
(
load
-
surface
.
mean
())
tractions
=
model
.
getTraction
()
true_displacements
=
model
.
getDisplacement
()
true_gap
=
true_displacements
-
surface
pressure
=
np
.
mean
(
tractions
)
contact_points
=
np
.
where
(
1.0e-10
>=
true_gap
)
contact_area
=
(
np
.
size
(
contact_points
)
/
2
)
/
float
(
grid_size
*
grid_size
)
print
(
'The contact area computed with the gap map is '
+
str
(
contact_area
))
hertz_contact_size
=
(
3
*
pressure
/
(
4
*
curvature
*
effective_modulus
))
**
(
1.
/
3.
)
print
(
'Exact Hertz contact radius is '
+
str
(
hertz_contact_size
))
hertz_area
=
np
.
pi
*
hertz_contact_size
**
2
print
(
'Exact Hertz contact area is '
+
str
(
hertz_area
))
area_error
=
np
.
abs
(
hertz_area
-
contact_area
)
/
hertz_area
print
(
"Area error: {}"
.
format
(
area_error
))
# Testing maximum pressure
max_pressure
=
tractions
.
max
()
hertz_max_pressure
=
(
6
*
pressure
*
effective_modulus
**
2
*
curvature
**
2
)
**
(
1.
/
3.
)
/
np
.
pi
pressure_error
=
np
.
abs
(
hertz_max_pressure
-
max_pressure
)
\
/
hertz_max_pressure
print
(
'Exact Hertz max pressure is '
+
str
(
hertz_max_pressure
))
print
(
'max pressure is '
+
str
(
np
.
max
(
tractions
)))
print
(
"Max pressure error: {}"
.
format
(
pressure_error
))
# Testing displacements
hertz_displacements
=
computeHertzDisplacement
(
effective_modulus
,
hertz_contact_size
,
hertz_max_pressure
,
grid_size
)
# Selecing only the points that are in contact
contact_indexes
=
[(
i
,
j
)
for
i
in
range
(
grid_size
)
for
j
in
range
(
grid_size
)
if
true_gap
[
i
,
j
]
<
1e-12
]
# Displacements of bem are centered around the mean of the whole surface
# and Hertz displacements are not centered, so we need to compute mean
# on the contact zone for both arrays
bem_mean
=
0.
hertz_mean
=
0.
for
index
in
contact_indexes
:
bem_mean
+=
true_displacements
[
index
]
hertz_mean
+=
hertz_displacements
[
index
]
bem_mean
/=
len
(
contact_indexes
)
hertz_mean
/=
len
(
contact_indexes
)
# Correction applied when computing error
correction
=
hertz_mean
-
bem_mean
# Computation of error of displacement in contact zone
error
=
0.
hertz_norm
=
0.
for
index
in
contact_indexes
:
error
+=
(
hertz_displacements
[
index
]
-
true_displacements
[
index
]
-
correction
)
**
2
hertz_norm
+=
(
hertz_displacements
[
index
]
-
hertz_mean
)
**
2
displacement_error
=
np
.
sqrt
(
error
/
hertz_norm
)
print
(
"Displacement error (in contact zone): {}"
.
format
(
displacement_error
))
assert
area_error
<
1e-2
\
and
pressure_error
<
1e-2
\
and
displacement_error
<
2e-3
tm
.
finalize
()
if
__name__
==
"__main__"
:
test_hertz_disp
()
Event Timeline
Log In to Comment