This repository has been archived by the owner on Feb 13, 2024. It is now read-only.
-
-
Notifications
You must be signed in to change notification settings - Fork 56
/
Copy pathpipe-crush-elastic.py
102 lines (87 loc) · 3.06 KB
/
pipe-crush-elastic.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
#!/usr/bin/env python3
import sys
import pycalculix as pyc
# Model of a pipe being crushed, quarter-symmetrym plane strain
proj_name = 'pipe-crush-elastic'
model = pyc.FeaModel(proj_name)
model.set_units('m') # this sets dist units to meters
# set whether or not to show gui plots
show_gui = True
if '-nogui' in sys.argv:
show_gui = False
# set element shape
eshape = 'quad'
if '-tri' in sys.argv:
eshape = 'tri'
# Define variables we'll use to draw part geometry
rad_outer = 0.1143 # pipe outer radius
pipe_wall_th = .00887 # wall thickness
rad_inner = rad_outer - pipe_wall_th
pipe_length = 10*rad_outer
wall_elements = 4.0
e_size = pipe_wall_th/wall_elements
# Draw pipe, x, y = radial, axial
pipe = pyc.Part(model)
pipe.goto(0.0, rad_inner)
arc_inner = pipe.draw_arc_angle(90, 0.0, 0.0)[0]
pipe.draw_line_to(rad_outer, 0.0)
arc_outer = pipe.draw_arc_angle(-90, 0.0, 0.0)[0]
pipe.draw_line_to(0.0, rad_inner)
num_outer = arc_outer.length()/e_size
num_inner = arc_inner.length()/e_size
num_arc_eles = int(round(0.5*(num_outer + num_inner),0))
# Draw plate, x, y = radial, axial
plate = pyc.Part(model)
plate.goto(rad_outer, 0.0)
plate.draw_line_rad(pipe_wall_th)
plate.draw_line_ax(rad_outer)
plate.draw_line_rad(-pipe_wall_th)
plate_bottom = plate.draw_line_to(rad_outer, 0.0)[0]
num_plate_eles = int(round(plate_bottom.length()/e_size,0))
# view model
model.plot_geometry(proj_name+'_geometry', display=show_gui)
model.plot_parts(proj_name+'_parts', display=show_gui)
model.plot_areas(proj_name+'_areas', display=show_gui)
model.plot_lines(proj_name+'_lines', display=show_gui)
# set loads and constraints
model.set_constr('fix',pipe.left,'y')
model.set_constr('fix',pipe.bottom,'x')
model.set_constr('fix',plate.left,'y')
ux_total = -.050
percent_displ = 50
model.set_constr('displ',plate.top,'x',ux_total*percent_displ/100.0)
# set part material
mat = pyc.Material('steel')
youngs = 210*(10**9)
mat.set_mech_props(7800, youngs, 0.3)
model.set_matl(mat, pipe)
model.set_matl(mat, plate)
# set contact
factor = 5 # can be between 5 and 50
kval = youngs*factor
model.set_contact_linear(plate_bottom, arc_outer, kval)
# set the element type and mesh database
model.set_eshape(eshape, 2)
model.set_etype('plstrain', pipe, pipe_length)
model.set_etype('plstrain', plate, pipe_length)
# set element divisions
model.set_ediv(['L1','L3', 'L4', 'L6'], wall_elements)
model.set_ediv(['L0','L2'], num_arc_eles) # set element divisions
model.set_ediv(['L5','L7'], num_plate_eles)
model.mesh(1.0, 'gmsh') # mesh 1.0 fineness, smaller is finer
model.plot_elements(proj_name+'_elem', display=show_gui)
model.plot_constraints(proj_name+'_constr', display=show_gui)
# make and solve the model
prob = pyc.Problem(model, 'struct')
prob.solve()
# Plot results
fields = 'Sx,Sy,S1,S2,S3,Seqv,ux,uy,utot' # store the fields to plot
fields = fields.split(',')
for field in fields:
fname = proj_name+'_'+field
prob.rfile.nplot(field, fname, display=False)
model.view.select(pipe)
model.view.allsel_under('parts')
for field in fields:
fname = proj_name+'_PIPE_'+field
prob.rfile.nplot(field, fname, display=False)