-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathsimple_convergence_study.jou
151 lines (126 loc) · 4.64 KB
/
simple_convergence_study.jou
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
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
#!python
import numpy
import matplotlib.pyplot as plt
#!cubit
reset
#beam geometry
create brick x 10 y 1 z 0.2
create Cylinder height 2 radius 1
create Cylinder height 2 radius 0.5
create Cylinder height 2 radius 0.25
move volume 2 y -0.9
move volume 3 x 3 y 0.45
move volume 4 x -3 y -0.6
subtract volume 2 from volume 1
subtract volume 3 from volume 1
subtract volume 4 from volume 1
block 1 add volume 1
# node sets
nodeset 1 add surface 4
nodeset 1 name "x0"
nodeset 2 add surface 6
nodeset 2 name "xl"
# material
create material "Steel" property_group "CalculiX-FEA"
modify material "Steel" scalar_properties "CCX_ELASTIC_USE_CARD" 1
modify material "Steel" scalar_properties "CCX_ELASTIC_ISO_USE_CARD" 1
modify material "Steel" matrix_property "CCX_ELASTIC_ISO_MODULUS_VS_POISSON_VS_TEMPERATURE" 210000 0.3 0
modify material "Steel" scalar_properties "CCX_PLASTIC_ISO_USE_CARD" 1
modify material "Steel" scalar_properties "CCX_EXPANSION_ISO_USE_CARD" 1
modify material "Steel" scalar_properties "CCX_CONDUCTIVITY_ISO_USE_CARD" 1
# section
ccx create section solid block 1 material 1
# vertices for ref and rot node
create vertex location on surface 6 center
create vertex location on surface 6 center
# rigid body constraint
ccx create constraint rigid body nodeset 2 ref 27 rot 28
# load
create force on vertex 27 force value 0.5 direction 0 -1 0
# bc
create displacement on nodeset 1 dof 1 dof 2 dof 3 fix 0
# outputs
ccx create fieldoutput name "node_output" node
ccx modify fieldoutput 1 node
ccx modify fieldoutput 1 node key_on RF U
ccx modify fieldoutput 1 node key_off CP DEPF DEPT DTF HCRI KEQ MACH MAXU MF NT PNT POT PRF PS PSF PT PTF PU RFL SEN TS TSF TT TTF TURB V VF
ccx create fieldoutput name "element_output" element
ccx modify fieldoutput 2 element
ccx modify fieldoutput 2 element key_on E S ERR
ccx modify fieldoutput 2 element key_off CEEQ ECD EMFB EMFE ENER HER HFL HFLF MAXE MAXS ME PEEQ PHS SF SMID SNEG SPOS SVF SDV THE ZZS
ccx create historyoutput name "node_output" node
ccx create historyoutput name "element_output" element
ccx modify historyoutput 1 node nodeset 3
ccx modify historyoutput 1 node key_on U RF
ccx modify historyoutput 1 node key_off NT TSF TTF PN PSF PTF MACH CP VF DEPF TURB MF RFL
ccx modify historyoutput 2 element block 1
ccx modify historyoutput 2 element key_on S E
ccx modify historyoutput 2 element key_off SVF ME PEEQ CEEQ ENER SDV HFL HFLF COORD ELSE ELKE EVOL EMAS EBHE CENT
# steps
ccx create step name "Static" static
#ccx modify step 1 parameter nlgeom_yes
#ccx modify step 1 static totaltimeatstart 0 initialtimeincrement 0.1 timeperiodofstep 1 minimumtimeincrement 0.01 maximumtimeincrement 0.1
ccx step 1 add load force 1
ccx step 1 add bc displacement 1
ccx step 1 add fieldoutput 1 2
#ccx step 1 add historyoutput 1 2
# mesh
block 1 element type HEX20
curve 11 interval 1
curve 11 scheme equal
#remesh and run again until convergence is reached
#!python
abs_con = 1e-1
rel_con = 1e-3
i_values = []
node_values = []
i = 1
i_values.append(i)
# mesh
cubit.cmd(f"volume 1 size 1")
cubit.cmd(f"mesh vertex all")
cubit.cmd(f"mesh volume all")
# jobs
cubit.cmd(f"ccx create job name 'simple_convergence_study_{i}'")
cubit.cmd(f"ccx run job {i} no_conversion")
cubit.cmd(f"ccx wait job {i}")
node_values.append(ccx.frd_get_node_value(i,cubit.get_closest_node(0,0.1,0.1),1,"STRESS","MISES"))
i+=1
i_values.append(i)
# mesh
cubit.cmd(f"delete mesh")
cubit.cmd(f"volume 1 size 0.8")
cubit.cmd(f"mesh vertex all")
cubit.cmd(f"mesh volume all")
# jobs
cubit.cmd(f"ccx create job name 'simple_convergence_study_{i}'")
cubit.cmd(f"ccx run job {i} no_conversion")
cubit.cmd(f"ccx wait job {i}")
node_values.append(ccx.frd_get_node_value(i,cubit.get_closest_node(0,0.1,0.1),1,"STRESS","MISES"))
i+=1
while (abs(1-abs(node_values[i-2]/node_values[i-3])) > rel_con) and (abs(node_values[i-2]-node_values[i-3]) > abs_con):
i_values.append(i)
# mesh
cubit.cmd(f"delete mesh")
cubit.cmd(f"volume 1 size {6/(i*i)}")
cubit.cmd(f"mesh vertex all")
cubit.cmd(f"mesh volume all")
# jobs
cubit.cmd(f"ccx create job name 'simple_convergence_study_{i}'")
cubit.cmd(f"ccx run job {i} no_conversion")
cubit.cmd(f"ccx wait job {i}")
node_values.append(ccx.frd_get_node_value(i,cubit.get_closest_node(0,0.1,0.1),1,"STRESS","MISES"))
i+=1
if i==20:
break
# Data for plotting
x_data = numpy.asarray(i_values)
y_data = numpy.asarray(node_values)
fig, ax = plt.subplots()
ax.plot(x_data, y_data)
ax.set(xlabel='Iteration', ylabel='Mises',title='Mises over iteration for node at coordinate x 0 y 0.1 z 0.1')
ax.grid()
plt.show()
for ii in range(len(node_values)):
print(str(i_values[ii]) + " " + str(node_values[ii]))
print("finished!")