14.3.1.2. Dambreak with Elastic Obstacle Analysis using moving mesh¶
The source code is shown below, which can be downloaded
here
.Run the source code in your favorite Python program.
The ParaView is needed to view the results. To view the displaced shape of fluid, use the “Warp By Vector” filter with scale factor = 1.0.
1import os
2import openseespy.opensees as ops
3
4# ------------------------------
5# Start of model generation
6# -----------------------------
7
8# remove existing model
9ops.wipe()
10
11# set modelbuilder
12ops.model('basic', '-ndm', 2, '-ndf', 3)
13
14# geometric
15L = 0.146
16H = 2*L
17H2 = 0.3
18b = 0.012
19h = 0.005
20alpha = 1.4
21Hb = 20.0*b/3.0
22tw = 3*h
23
24# material
25rho = 1000.0
26mu = 0.0001
27b1 = 0.0
28b2 = -9.81
29thk = 0.012
30kappa = -1.0
31
32rhos = 2500.0
33A = thk*thk
34E = 1e6
35Iz = thk*thk*thk*thk/12.0
36bmass = A*Hb*rhos
37
38# analysis
39dtmax = 1e-3
40dtmin = 1e-6
41totaltime = 1.0
42
43filename = 'obstacle'
44
45# recorder
46if not os.path.exists(filename):
47 os.makedirs(filename)
48ops.recorder('PVD', filename, 'disp', 'vel', 'pressure')
49
50# nodes
51ops.node(1, 0.0, 0.0)
52ops.node(2, L, 0.0)
53ops.node(3, L, H, '-ndf', 2)
54ops.node(4, 0.0, H)
55ops.node(5, 0.0, H2)
56ops.node(6, 4*L, 0.0)
57ops.node(7, 4*L, H2)
58ops.node(8, -tw, H2)
59ops.node(9, -tw, -tw)
60ops.node(10, 4*L+tw, -tw)
61ops.node(11, 4*L+tw, H2)
62ops.node(12, 2*L, 0.0)
63ops.node(13, 2*L, Hb)
64
65# ids for meshing
66wall_id = 1
67beam_id = 2
68water_bound_id = -1
69water_body_id = -2
70
71# transformation
72transfTag = 1
73ops.geomTransf('Corotational', transfTag)
74
75# section
76secTag = 1
77ops.section('Elastic', secTag, E, A, Iz)
78
79# beam integration
80inteTag = 1
81numpts = 2
82ops.beamIntegration('Legendre', inteTag, secTag, numpts)
83
84# beam mesh
85beamTag = 6
86ndf = 3
87ops.mesh('line', beamTag, 2, 12, 13, beam_id, ndf, h, 'dispBeamColumn', transfTag, inteTag)
88
89ndmass = bmass/len(ops.getNodeTags('-mesh', beamTag))
90
91for nd in ops.getNodeTags('-mesh', beamTag):
92 ops.mass(nd, ndmass, ndmass, 0.0)
93
94# fluid mesh
95fluidTag = 4
96ndf = 2
97ops.mesh('line', 1, 10, 4,5,8,9,10,11,7,6,12,2, wall_id, ndf, h)
98ops.mesh('line', 2, 3, 2,1,4, wall_id, ndf, h)
99ops.mesh('line', 3, 3, 2,3,4, water_bound_id, ndf, h)
100
101eleArgs = ['PFEMElementBubble',rho,mu,b1,b2,thk,kappa]
102ops.mesh('tri', fluidTag, 2, 2,3, water_body_id, ndf, h, *eleArgs)
103
104# wall mesh
105wallTag = 5
106ops.mesh('tri', wallTag, 2, 1,2, wall_id, ndf, h)
107
108for nd in ops.getNodeTags('-mesh', wallTag):
109 ops.fix(nd, 1,1,1)
110
111# save the original modal
112ops.record()
113
114# create constraint object
115ops.constraints('Plain')
116
117# create numberer object
118ops.numberer('Plain')
119
120# create convergence test object
121ops.test('PFEM', 1e-5, 1e-5, 1e-5, 1e-5, 1e-15, 1e-15, 20, 3, 1, 2)
122
123# create algorithm object
124ops.algorithm('Newton')
125
126# create integrator object
127ops.integrator('PFEM')
128
129# create SOE object
130ops.system('PFEM')
131
132# create analysis object
133ops.analysis('PFEM', dtmax, dtmin, b2)
134
135# analysis
136while ops.getTime() < totaltime:
137
138 # analysis
139 if ops.analyze() < 0:
140 break
141
142 ops.remesh(alpha)