Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F78460526
test_westergaard.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
Wed, Aug 21, 00:43
Size
4 KB
Mime Type
text/x-python
Expires
Fri, Aug 23, 00:43 (2 d)
Engine
blob
Format
Raw Data
Handle
20049716
Attached To
rTAMAAS tamaas
test_westergaard.py
View Options
#!/usr/bin/env python
# coding: utf-8
# -----------------------------------------------------------------------------
# @author Lucas Frérot <lucas.frerot@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
division
,
print_function
import
sys
import
numpy
as
np
from
numpy.linalg
import
norm
from
numpy
import
sin
,
cos
,
pi
,
sqrt
,
log
import
tamaas
as
tm
def
constructSinProfile
(
size
,
mode
,
amplitude
):
x
=
np
.
linspace
(
0
,
1
,
size
,
endpoint
=
False
)
y
=
np
.
linspace
(
0
,
1
,
size
,
endpoint
=
False
)
x
,
y
=
np
.
meshgrid
(
x
,
y
)
surface
=
amplitude
*
np
.
sin
(
2
*
np
.
pi
*
x
*
mode
)
return
surface
.
copy
()
def
main
():
grid_size
=
2187
# Power of 3 for better accuracy
lamda
=
1.
delta
=
0.1
effective_modulus
=
1.
p_star
=
np
.
pi
*
effective_modulus
*
delta
/
lamda
# Full contact load
load
=
0.1
*
p_star
surface
=
constructSinProfile
(
grid_size
,
1
/
lamda
,
delta
)
bem
=
tm
.
BemPolonski
(
surface
)
bem
.
setEffectiveModulus
(
effective_modulus
)
bem
.
computeEquilibrium
(
1e-12
,
load
)
tractions
=
bem
.
getTractions
()
displacements
=
bem
.
getDisplacements
()
# Testing contact area against Westergaard solution
contact_area
=
np
.
sum
(
tractions
>
0
)
/
grid_size
**
2
westergaard_contact_size
=
2
*
lamda
/
np
.
pi
*
np
.
arcsin
(
np
.
sqrt
(
load
/
p_star
))
area_error
=
np
.
abs
(
contact_area
-
westergaard_contact_size
/
lamda
)
/
(
westergaard_contact_size
/
lamda
)
print
(
"Area error: {}"
.
format
(
area_error
))
# Testing agains pressure from Westergaard solution
x
=
np
.
linspace
(
0
,
1
,
grid_size
,
endpoint
=
False
)
y
=
np
.
linspace
(
0
,
1
,
grid_size
,
endpoint
=
False
)
x
,
y
=
np
.
meshgrid
(
x
,
y
)
x
-=
1
/
4
a
=
westergaard_contact_size
/
2
pressure
=
np
.
zeros_like
(
tractions
)
pressure
=
2
*
load
*
np
.
cos
(
np
.
pi
*
x
/
lamda
)
/
np
.
sin
(
np
.
pi
*
a
/
lamda
)
**
2
*
\
np
.
sqrt
(
np
.
sin
(
np
.
pi
*
a
/
lamda
)
**
2
-
np
.
sin
(
np
.
pi
*
x
/
lamda
)
**
2
)
pressure
[
np
.
abs
(
x
)
>
a
]
=
0
pressure_error
=
norm
(
pressure
-
tractions
)
/
norm
(
pressure
)
print
(
"Pressure error: {}"
.
format
(
pressure_error
))
# Energetic error (displacement reference in Johnson)
displacements
-=
displacements
.
max
()
psi
=
np
.
pi
*
np
.
abs
(
x
)
/
lamda
psi_a
=
np
.
pi
*
a
/
lamda
disp_w
=
load
*
lamda
/
(
pi
*
effective_modulus
*
sin
(
psi_a
))
*
(
cos
(
2
*
psi
)
+
2
*
sin
(
psi
)
*
sqrt
(
sin
(
psi
)
**
2
-
sin
(
psi_a
)
**
2
)
-
2
*
sin
(
psi_a
)
**
2
*
log
((
sin
(
psi
)
+
sqrt
(
sin
(
psi
)
**
2
-
sin
(
psi_a
)
**
2
))
/
sin
(
psi_a
)))
disp_w
[
np
.
abs
(
x
)
<
a
]
=
load
*
lamda
/
(
pi
*
effective_modulus
*
sin
(
psi_a
))
*
\
cos
(
2
*
psi
[
np
.
abs
(
x
)
<
a
])
disp_w
-=
disp_w
.
max
()
energy_error
=
np
.
sum
((
tractions
-
pressure
)
*
(
displacements
-
disp_w
))
/
\
norm
(
pressure
*
disp_w
)
/
grid_size
**
2
print
(
"Energy error: {}"
.
format
(
energy_error
))
if
area_error
>
1e-4
or
pressure_error
>
5e-4
or
energy_error
>
2e-8
:
return
1
# import matplotlib.pyplot as plt
# plt.imshow(pressure, cmap='viridis')
# plt.figure()
# plt.imshow(tractions, cmap='viridis')
# plt.figure()
# plt.imshow(displacements, cmap='viridis')
# plt.figure()
# plt.imshow(disp_w, cmap='viridis')
# plt.show()
return
0
if
__name__
==
"__main__"
:
sys
.
exit
(
main
())
Event Timeline
Log In to Comment