-
Notifications
You must be signed in to change notification settings - Fork 107
/
timoshenko.py
110 lines (91 loc) · 2.73 KB
/
timoshenko.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
__doc__ = """Timoshenko beam validation case, for detailed explanation refer to
Gazzola et. al. R. Soc. 2018 section 3.4.3 """
import numpy as np
import sys
# FIXME without appending sys.path make it more generic
sys.path.append("../../")
from elastica import *
from examples.TimoshenkoBeamCase.timoshenko_postprocessing import plot_timoshenko
class TimoshenkoBeamSimulator(BaseSystemCollection, Constraints, Forcing):
pass
timoshenko_sim = TimoshenkoBeamSimulator()
final_time = 5000
# Options
PLOT_FIGURE = True
SAVE_FIGURE = False
SAVE_RESULTS = False
ADD_UNSHEARABLE_ROD = True
# setting up test params
n_elem = 100
start = np.zeros((3,))
direction = np.array([0.0, 0.0, 1.0])
normal = np.array([0.0, 1.0, 0.0])
base_length = 3.0
base_radius = 0.25
base_area = np.pi * base_radius ** 2
density = 5000
nu = 0.1
E = 1e6
# For shear modulus of 1e4, nu is 99!
poisson_ratio = 99
shear_modulus = E / (poisson_ratio + 1.0)
shearable_rod = CosseratRod.straight_rod(
n_elem,
start,
direction,
normal,
base_length,
base_radius,
density,
nu,
E,
shear_modulus=shear_modulus,
)
timoshenko_sim.append(shearable_rod)
timoshenko_sim.constrain(shearable_rod).using(
OneEndFixedBC, constrained_position_idx=(0,), constrained_director_idx=(0,)
)
end_force = np.array([-15.0, 0.0, 0.0])
timoshenko_sim.add_forcing_to(shearable_rod).using(
EndpointForces, 0.0 * end_force, end_force, ramp_up_time=final_time / 2.0
)
if ADD_UNSHEARABLE_ROD:
# Start into the plane
unshearable_start = np.array([0.0, -1.0, 0.0])
shear_modulus = E / (-0.7 + 1.0)
unshearable_rod = CosseratRod.straight_rod(
n_elem,
unshearable_start,
direction,
normal,
base_length,
base_radius,
density,
nu,
E,
# Unshearable rod needs G -> inf, which is achievable with -ve poisson ratio
shear_modulus=shear_modulus,
)
timoshenko_sim.append(unshearable_rod)
timoshenko_sim.constrain(unshearable_rod).using(
OneEndFixedBC, constrained_position_idx=(0,), constrained_director_idx=(0,)
)
timoshenko_sim.add_forcing_to(unshearable_rod).using(
EndpointForces, 0.0 * end_force, end_force, ramp_up_time=final_time / 2.0
)
timoshenko_sim.finalize()
timestepper = PositionVerlet()
# timestepper = PEFRL()
dl = base_length / n_elem
dt = 0.01 * dl
total_steps = int(final_time / dt)
print("Total steps", total_steps)
integrate(timestepper, timoshenko_sim, final_time, total_steps)
if PLOT_FIGURE:
plot_timoshenko(shearable_rod, end_force, SAVE_FIGURE, ADD_UNSHEARABLE_ROD)
if SAVE_RESULTS:
import pickle
filename = "Timoshenko_beam_data.dat"
file = open(filename, "wb")
pickle.dump(shearable_rod, file)
file.close()