Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F70895072
test_hertz_adhesion.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
Mon, Jul 8, 03:43
Size
4 KB
Mime Type
text/x-python
Expires
Wed, Jul 10, 03:43 (2 d)
Engine
blob
Format
Raw Data
Handle
18883646
Attached To
rTAMAAS tamaas
test_hertz_adhesion.py
View Options
# -*- coding: utf-8 -*-
# -----------------------------------------------------------------------------
# @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/>.
# -----------------------------------------------------------------------------
from
__future__
import
print_function
,
division
import
sys
import
tamaas
as
tm
import
numpy
as
np
import
scipy.optimize
as
scio
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
=
-
(
x
**
2
+
y
**
2
)
/
(
2
*
radius
)
return
surface
.
copy
()
def
main
():
tm
.
initialize
()
grid_size
=
512
curvature
=
0.5
effective_modulus
=
1.
surface_energy
=
0.0001
rho
=
0.002
radius
=
1.
/
curvature
K
=
4
*
effective_modulus
/
3
KI
=
np
.
sqrt
(
surface_energy
*
2
*
effective_modulus
)
tabor
=
surface_energy
**
(
2
/
3
)
/
(
rho
*
effective_modulus
**
(
2
/
3
))
*
(
1
/
curvature
)
**
(
1
/
3
)
sig0
=
surface_energy
/
(
rho
)
lam
=
2
*
sig0
/
(
np
.
pi
*
surface_energy
*
K
**
2
/
radius
)
**
(
1
/
3
)
surface
=
constructHertzProfile
(
grid_size
,
curvature
)
bem
=
tm
.
BemGigipol
(
surface
)
bem
.
setDumpFreq
(
100
)
functional
=
tm
.
MaugisAdhesionFunctional
(
bem
)
functional
.
setParameter
(
'rho'
,
rho
)
functional
.
setParameter
(
'surface_energy'
,
surface_energy
)
bem
.
setEffectiveModulus
(
effective_modulus
)
bem
.
addFunctional
(
functional
)
bem
.
computeEquilibrium
(
1e-8
,
0.00116480372339
)
tractions
=
bem
.
getTractions
()
true_gap
=
bem
.
getGap
()
true_displacements
=
bem
.
getTrueDisplacements
()
adhesive_pressure
=
np
.
zeros
((
grid_size
,
grid_size
))
for
i
in
np
.
arange
(
0
,
grid_size
,
1
):
for
j
in
np
.
arange
(
0
,
grid_size
,
1
):
if
(
rho
>
true_gap
[
i
,
j
]):
adhesive_pressure
[
i
,
j
]
=
surface_energy
/
rho
# computed solution
p0
=
np
.
mean
(
tractions
)
p_max
=
np
.
max
(
tractions
)
contact_points
=
np
.
where
(
true_gap
==
0
)
contact_area
=
np
.
size
(
contact_points
)
/
2
/
(
float
(
grid_size
)
**
2
)
ag
=
(
contact_area
/
np
.
pi
)
**
0.5
adh_points
=
np
.
where
(
adhesive_pressure
!=
0.
)
adh_area
=
np
.
size
(
adh_points
)
/
2
/
(
float
(
grid_size
)
**
2
)
cg
=
(
adh_area
/
np
.
pi
)
**
0.5
delta
=
ag
**
2
/
radius
-
8
*
sig0
/
(
3
*
K
)
*
np
.
sqrt
(
cg
**
2
-
ag
**
2
)
# Exact solution
beta
=
(
1
-
np
.
exp
(
lam
/-
0.924
))
/
1.02
a_chap
=
1.54
+
((
2.28
*
lam
**
1.3
-
1
)
/
(
2.28
*
lam
**
1.3
+
1
))
*
0.279
L_chap
=
-
7.
/
4.
+
((
4.04
*
lam
**
1.4
-
1
)
/
(
4.04
*
lam
**
1.4
+
1
))
/
4.
a0
=
a_chap
*
(
np
.
pi
*
surface_energy
*
radius
**
2
/
K
)
**
(
1.
/
3.
)
pc
=
L_chap
*
np
.
pi
*
surface_energy
*
radius
P
=
-
0.00078
ac
=
a0
*
((
beta
+
np
.
sqrt
(
1
-
P
/
pc
))
/
(
1
+
beta
))
**
(
2.
/
3.
)
def
fonction1
(
m
):
return
np
.
sqrt
(
m
**
2
-
1
)
+
m
**
2
*
np
.
arctan
(
np
.
sqrt
(
m
**
2
-
1
))
-
np
.
pi
*
KI
/
(
np
.
sqrt
(
np
.
pi
*
ac
)
*
sig0
)
m
=
scio
.
brenth
(
fonction1
,
1.
,
2.
)
cc
=
m
*
ac
delta_cos
=
ac
**
2
/
radius
-
8
*
sig0
/
(
3
*
K
)
*
np
.
sqrt
(
cc
**
2
-
ac
**
2
)
x
=
np
.
linspace
(
-
0.5
,
0.5
,
grid_size
)
y
=
np
.
linspace
(
-
0.5
,
0.5
,
grid_size
)
x
,
y
=
np
.
meshgrid
(
x
,
y
)
a
=
ac
c
=
cc
# stress
sig_hertz
=
np
.
zeros
((
grid_size
,
grid_size
))
disp
=
np
.
zeros
((
grid_size
,
grid_size
))
for
i
in
np
.
arange
(
0
,
grid_size
,
1
):
for
j
in
np
.
arange
(
0
,
grid_size
,
1
):
if
a
**
2
>
(
x
[
i
,
j
]
**
2
+
y
[
i
,
j
]
**
2
):
disp
[
i
,
j
]
=
3
*
K
/
(
2
*
np
.
pi
*
radius
)
*
np
.
sqrt
(
a
**
2
-
(
x
[
i
,
j
]
**
2
+
y
[
i
,
j
]
**
2
))
-
2
*
sig0
/
np
.
pi
*
np
.
arctan
(
np
.
sqrt
((
c
**
2
-
a
**
2
)
/
(
a
**
2
-
(
x
[
i
,
j
]
**
2
+
y
[
i
,
j
]
**
2
))))
if
c
**
2
>
(
x
[
i
,
j
]
**
2
+
y
[
i
,
j
]
**
2
)
>
a
**
2
:
disp
[
i
,
j
]
=-
sig0
sig_hertz
=
disp
a_error
=
abs
((
ac
-
ag
)
/
ac
)
c_error
=
abs
((
cc
-
cg
)
/
cc
)
delta_error
=
abs
((
delta
-
delta_cos
)
/
delta
)
pressure_error
=
abs
((
p0
-
np
.
mean
(
sig_hertz
))
/
np
.
mean
(
sig_hertz
))
if
a_error
>
1.e-2
or
c_error
>
1.2e-2
or
delta_error
>
5.e-2
or
pressure_error
>
1.e-1
:
for
err
in
[
a_error
,
c_error
,
delta_error
,
pressure_error
]:
print
(
err
)
return
1
tm
.
finalize
()
return
0
if
__name__
==
"__main__"
:
sys
.
exit
(
main
())
Event Timeline
Log In to Comment