Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F79664899
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
Tue, Aug 27, 15:34
Size
4 KB
Mime Type
text/x-python
Expires
Thu, Aug 29, 15:34 (1 d, 22 h)
Engine
blob
Format
Raw Data
Handle
20218025
Attached To
rTAMAAS tamaas
test_hertz_disp.py
View Options
# -*- coding: utf-8 -*-
# -----------------------------------------------------------------------------
# @author Lucas Frérot <lucas.frerot@epfl.ch>
# @author Valentine Rey <valentine.rey@epfl.ch>
#
# @section LICENSE
#
# Copyright (©) 2016 EPFL (Ecole Polytechnique Fédérale de
# Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des
# Solides)
#
# Tamaas is free software: you can redistribute it and/or modify it under the
# terms of the GNU Lesser General Public License as published by the Free
# Software Foundation, either version 3 of the License, or (at your option) any
# later version.
#
# Tamaas 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 Lesser General Public License for more
# details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
# -----------------------------------------------------------------------------
import
sys
from
tamaas
import
*
import
numpy
as
np
def
constructHertzProfile
(
size
,
curvature
):
radius
=
1.
/
curvature
x
=
np
.
linspace
(
-
0.5
,
0.5
,
size
)
y
=
np
.
linspace
(
-
0.5
,
0.5
,
size
)
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
main
():
grid_size
=
512
curvature
=
0.5
effective_modulus
=
1.
load
=
0.12
surface_energy
=
0.
rho
=
1.0e-3
surface
=
constructHertzProfile
(
grid_size
,
curvature
)
bem
=
BemGigipol
(
surface
)
bem
.
setDumpFreq
(
100
)
functional
=
ExponentialAdhesionFunctional
(
bem
)
functional
.
setParameter
(
'rho'
,
rho
)
functional
.
setParameter
(
'surface_energy'
,
surface_energy
)
bem
.
setEffectiveModulus
(
effective_modulus
)
bem
.
addFunctional
(
functional
)
bem
.
computeEquilibrium
(
1e-12
,
load
)
tractions
=
bem
.
getTractions
()
true_gap
=
bem
.
getGap
()
true_displacements
=
bem
.
getTrueDisplacements
()
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
,
true_gap
[
i
,
j
]
==
0
)
for
i
in
range
(
grid_size
)
for
j
in
range
(
grid_size
)]
contact_indexes
=
map
(
lambda
x
:
x
[
0
:
2
],
filter
(
lambda
x
:
x
[
2
],
contact_indexes
))
# 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
)
if
area_error
>
1e-2
or
pressure_error
>
1e-2
or
displacement_error
>
2e-3
:
return
1
return
0
if
__name__
==
"__main__"
:
sys
.
exit
(
main
())
Event Timeline
Log In to Comment