14.3.2.1. Dambreak Analysis using background 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 os.path
3import openseespy.opensees as ops
4
5# ------------------------------
6# Start of model generation
7# -----------------------------
8
9# wipe all previous objects
10ops.wipe()
11
12# create a model with fluid
13ops.model('basic', '-ndm', 2, '-ndf', 2)
14
15# geometric
16L = 0.146
17H = L * 2
18h = L / 40
19
20# number of particles per cell in each direction
21numx = 3.0
22numy = 3.0
23
24# material
25rho = 1000.0
26mu = 0.0001
27b1 = 0.0
28b2 = -9.81
29thk = 0.012
30kappa = -1.0
31
32# analysis
33dtmax = 1e-3
34dtmin = 1e-3
35totaltime = 1.0
36filename = 'dambreak-bg'
37
38# recorder
39ops.recorder('BgPVD', filename, 'disp', 'vel', 'pressure', '-dT', 1e-3)
40if not os.path.exists(filename):
41 os.makedirs(filename)
42
43# fluid particles
44ndf = 2
45
46# total number of particles in each direction
47nx = round(L / h * numx)
48ny = round(H / h * numy)
49
50# create particles
51eleArgs = ['PFEMElementBubble', rho, mu, b1, b2, thk, kappa]
52partArgs = ['quad', 0.0, 0.0, L, 0.0, L, H, 0.0, H, nx, ny]
53parttag = 1
54ops.mesh('part', parttag, *partArgs, *eleArgs, '-vel', 0.0, 0.0)
55
56print('num particles =', nx * ny)
57
58# wall
59ops.node(1, 0.0, H)
60ops.node(2, 0.0, 0.0)
61ops.node(3, 4 * L, 0.0)
62ops.node(4, 4 * L, H)
63
64walltag = 2
65wallid = 1
66ops.mesh('line', walltag, 4, 1, 2, 3, 4, wallid, ndf, h)
67
68wallnodes = ops.getNodeTags('-mesh', walltag)
69
70for nd in wallnodes:
71 ops.fix(nd, 1, 1)
72
73# background mesh
74lower = [-h, -h]
75upper = [4 * L + L, H + L]
76
77ops.mesh('bg', h, *lower, *upper,
78 '-structure', wallid, len(wallnodes), *wallnodes)
79
80# create constraint object
81ops.constraints('Plain')
82
83# create numberer object
84ops.numberer('Plain')
85
86# create convergence test object
87ops.test('PFEM', 1e-5, 1e-5, 1e-5, 1e-5, 1e-5, 1e-5, 10, 3, 1, 2)
88
89# create algorithm object
90ops.algorithm('Newton')
91
92# create integrator object
93ops.integrator('PFEM', 0.5, 0.25)
94
95# create SOE object
96ops.system('PFEM')
97# ops.system('PFEM', '-mumps) Linux version can use mumps
98
99# create analysis object
100ops.analysis('PFEM', dtmax, dtmin, b2)
101
102# analysis
103while ops.getTime() < totaltime:
104
105 # analysis
106 if ops.analyze() < 0:
107 break
108
109 ops.remesh()
110
111print("==========================================")