Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F79687969
hertz_3D.cc
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, 18:20
Size
5 KB
Mime Type
text/x-c
Expires
Thu, Aug 29, 18:20 (2 d)
Engine
blob
Format
Raw Data
Handle
20146348
Attached To
rAKA akantu
hertz_3D.cc
View Options
/**
* @file hertz_3D.cc
*
* @author Alejandro M. Aragón <alejandro.aragon@epfl.ch>
*
* @date creation: Sun Oct 19 2014
*
* @brief File used to obtain 3D implicit contact results and compare with
* Hertz theory
*
* @section LICENSE
*
* Copyright (©) 2015 EPFL (Ecole Polytechnique Fédérale de Lausanne) Laboratory
* (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu 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.
*
* Akantu 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 Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
#include <chrono>
#include "contact_impl.hh"
using
namespace
akantu
;
using
std
::
cout
;
using
std
::
endl
;
using
std
::
setw
;
using
std
::
setprecision
;
int
main
(
int
argc
,
char
*
argv
[])
{
// set dimension
static
const
UInt
dim
=
3
;
// type definitions
typedef
Point
<
dim
>
point_type
;
typedef
BoundingBox
<
dim
>
bbox_type
;
typedef
SolidMechanicsModel
model_type
;
typedef
Contact
<
dim
,
MasterAssignator
,
SelectResolution
<
_static
,
_augmented_lagrangian
>
>
contact_type
;
typedef
std
::
chrono
::
high_resolution_clock
clock
;
typedef
std
::
chrono
::
seconds
seconds
;
clock
::
time_point
t0
=
clock
::
now
();
initialize
(
"steel.dat"
,
argc
,
argv
);
// create mesh
Mesh
mesh
(
dim
);
// read mesh
mesh
.
read
(
"hertz_3D.msh"
);
// create model
model_type
model
(
mesh
);
SolidMechanicsModelOptions
opt
(
_static
);
// initialize material
model
.
initFull
(
opt
);
// create data structure that holds contact data
contact_type
cd
(
argc
,
argv
,
model
);
// optimal value of penalty multiplier
cd
[
Alpha
]
=
0.05
;
cd
[
Multiplier_tol
]
=
1.e-2
;
cd
[
Newton_tol
]
=
1.e-2
;
cd
[
Verbose
]
=
true
;
// set Paraview output resluts
model
.
setBaseName
(
"contact"
);
model
.
addDumpFieldVector
(
"displacement"
);
// call update current position to be able to call later
// the function to get current positions
model
.
updateCurrentPosition
();
// get physical names from Gmsh file
mesh
.
createGroupsFromMeshData
<
std
::
string
>
(
"physical_names"
);
// set-up bounding box to include slave nodes that lie inside it
Real
l1
=
1.
;
Real
l2
=
0.2
;
Real
l3
=
1.
;
point_type
c1
(
-
l1
/
2
,
-
l2
/
2
,
-
l3
/
2
);
point_type
c2
(
l1
/
2
,
l2
/
2
,
l3
/
2
);
bbox_type
bb
(
c1
,
c2
);
// get areas for the nodes of the circle
// this is done by applying a unit pressure to the contact surface elements
model
.
applyBC
(
BC
::
Neumann
::
FromHigherDim
(
Matrix
<
Real
>::
eye
(
3
,
1.
)),
"contact_surface"
);
Array
<
Real
>
&
areas
=
model
.
getForce
();
// loop over contact surface nodes and store node areas
ElementGroup
&
eg
=
mesh
.
getElementGroup
(
"contact_surface"
);
Array
<
Real
>
&
coords
=
mesh
.
getNodes
();
cout
<<
"- Adding areas to slave nodes. "
<<
endl
;
for
(
auto
nit
=
eg
.
node_begin
();
nit
!=
eg
.
node_end
();
++
nit
)
{
point_type
p
(
&
coords
(
*
nit
));
// ignore slave node if it doesn't lie within the bounding box
if
(
!
(
bb
&
p
))
continue
;
cd
.
addSlave
(
*
nit
);
// compute area contributing to the slave node
Real
a
=
0.
;
for
(
UInt
i
=
0
;
i
<
dim
;
++
i
)
a
+=
pow
(
areas
(
*
nit
,
i
),
2.
);
cd
.
addArea
(
*
nit
,
sqrt
(
a
));
}
// set force value to zero
areas
.
clear
();
// add master surface to find pairs
cd
.
searchSurface
(
"rigid_surface"
);
// apply boundary conditions for the rigid plane
model
.
applyBC
(
BC
::
Dirichlet
::
FixedValue
(
0.
,
_x
),
"bottom_body"
);
model
.
applyBC
(
BC
::
Dirichlet
::
FixedValue
(
0.
,
_y
),
"bottom_body"
);
model
.
applyBC
(
BC
::
Dirichlet
::
FixedValue
(
0.
,
_z
),
"bottom_body"
);
// block z-disp in extreme points of top surface
model
.
getBlockedDOFs
()(
1
,
2
)
=
true
;
model
.
getBlockedDOFs
()(
2
,
2
)
=
true
;
// block x-disp in extreme points of top surface
model
.
getBlockedDOFs
()(
3
,
0
)
=
true
;
model
.
getBlockedDOFs
()(
4
,
0
)
=
true
;
const
size_t
steps
=
30
;
Real
data
[
3
][
steps
];
// store results for printing
Real
step
=
0.001
;
// top displacement increment
size_t
k
=
0
;
for
(
Real
delta
=
0
;
delta
<=
step
*
steps
;
delta
+=
step
)
{
// apply displacement to the top surface of the half-sphere
model
.
applyBC
(
BC
::
Dirichlet
::
FixedValue
(
-
delta
,
_y
),
"top_surface"
);
// solve contact step, this function also dumps Paraview files
solveContactStep
<
_uzawa
>
(
cd
);
// solveContactStep<_generalized_newton>(cd);
data
[
0
][
k
]
=
delta
;
data
[
1
][
k
]
=
cd
.
getForce
();
data
[
2
][
k
]
=
cd
.
getMaxPressure
();
++
k
;
}
// print results
size_t
w
=
14
;
cout
<<
setprecision
(
4
);
cout
<<
endl
<<
setw
(
w
)
<<
"Disp."
<<
setw
(
w
)
<<
"Force"
<<
setw
(
w
)
<<
"Max pressure"
<<
endl
;
for
(
size_t
i
=
0
;
i
<
steps
;
++
i
)
cout
<<
setw
(
w
)
<<
data
[
0
][
i
]
<<
setw
(
w
)
<<
data
[
1
][
i
]
<<
setw
(
w
)
<<
data
[
2
][
i
]
<<
endl
;
clock
::
time_point
t1
=
clock
::
now
();
seconds
total_s
=
std
::
chrono
::
duration_cast
<
seconds
>
(
t1
-
t0
);
cout
<<
"*** INFO *** Simulation took "
<<
total_s
.
count
()
<<
" s"
<<
endl
;
// finalize simulation
finalize
();
return
EXIT_SUCCESS
;
}
Event Timeline
Log In to Comment