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
|
# from https://github.com/seamplex/pipe-linearize
# convergence study of linearized stresses in an infinite pipe
# with respect to the number of elements in the pipe thickness
PROBLEM mechanical PC mumps
# PETSC_OPTIONS -mg_levels_pc_type sor
# t0 = wall_time() # initial wall time
# read mesh according to shape $1, order $2 and number of elements through thickness $2
DEFAULT_ARGUMENT_VALUE 2 2
DEFAULT_ARGUMENT_VALUE 3 2
READ_MESH pipe-$1-$2-$3.msh
# problem parameters for
# 12"-inch schedule 100
b = 323.8/2 # external radius [ mm ]
a = b-21.5 # internal radius [ mm ]
l = 2*(b-a) # axial length [ mm ]
E = 200e3 # Young modulus [ MPa ]
nu = 0.3 # Poisson's ratio [ non-dimensional ]
p = 10 # internal pressure [ MPa ]
# ------------------------------------------------------------------------
# definition of analytical solutions for comparison from
# <http://eprints.whiterose.ac.uk/110536/1/art%253A10.1007%252Fs00707-016-1762-7.pdf>
ur(x,y,z) = (p*a^2*sqrt(y^2+z^2))/(E*(b^2-a^2)) * ((1-2*nu)*(1+nu) + (1+nu)*b^2/(y^2+z^2))
sigmal(x,y,z) = 2*nu*p*a^2/(b^2-a^2)
sigmar(x,y,z) = p*a^2/(b^2-a^2) * (1 - b^2/(y^2+z^2))
sigmatheta(x,y,z) = p*a^2/(b^2-a^2) * (1 + b^2/(y^2+z^2))
# principal stresses along the radial coordinate (may be y or z)
sigma1_anal(r) = sigmatheta(0,0,r)
sigma2_anal(r) = sigmal(0,0,r)
sigma3_anal(r) = sigmar(0,0,r)
# computation of main membrane stresses
M1 = 1/(b-a)*integral(sigma1_anal(r), r, a, b)
M2 = 1/(b-a)*integral(sigma2_anal(r), r, a, b)
M3 = 1/(b-a)*integral(sigma3_anal(r), r, a, b)
# von mises and tresca membrane stresses
Mt_anal = max(abs(M1-M2), abs(M2-M3), abs(M3-M1))
Mv_anal = sqrt(((M1-M2)^2 + (M2-M3)^2 + (M3-M1)^2)/2)
# computation of membrane plus bending stresses
MB1 = M1 + 6/(b-a)^2*integral(sigma1_anal(r)*((a+b)/2-r), r, a, b)
MB2 = M2 + 6/(b-a)^2*integral(sigma2_anal(r)*((a+b)/2-r), r, a, b)
MB3 = M3 + 6/(b-a)^2*integral(sigma3_anal(r)*((a+b)/2-r), r, a, b)
# von mises and tresca
MBv_anal = sqrt(((MB1-MB2)^2 + (MB2-MB3)^2 + (MB3-MB1)^2)/2)
MBt_anal = max(abs(MB1-MB2), abs(MB2-MB3), abs(MB3-MB1))
# set boundary conditions
BC pressure pressure=p
BC front tangential radial
BC back tangential radial
SOLVE_PROBLEM
# write distribution of results in gmsh format (optional)
# WRITE_MESH pipe-out-$1-$2-$3.msh VECTOR u v w sigma1 sigma2 sigma3 sigmax sigmay sigmaz tauxy tauyz tauzx
# write same thing as ASCII data
# PRINT_FUNCTION v sigma1 sigma2 sigma3 MIN 0 a 0 MAX 0 b 0 NSTEPS 1 200 1 FILE pipe-dist-$1-$2-$3.dat
# compute linearized stresses for both von Mises and Tresca
LINEARIZE_STRESS FROM 0 a 0 TO 0 b 0 M Mv MB MBv Mt Mt MBt MBt
# compute L2 errors
volume = pi*(b^2-a^2)*l
h = (volume/cells)^(1/3)
ur_fea(x,y,z) = sqrt(v(x,y,z)^2+w(x,y,z)^2)
INTEGRATE (ur_fea(x,y,z)-ur(x,y,z))^2 RESULT e2ur
INTEGRATE (sigma1(x,y,z)-sigmatheta(x,y,z))^2 RESULT e2sigma1
INTEGRATE (sigma3(x,y,z)-sigmar(x,y,z))^2 RESULT e2sigma3
error_ur = sqrt(e2ur)/volume
error_sigma1 = sqrt(e2sigma1)/volume
error_sigma3 = sqrt(e2sigma3)/volume
error_Mv = abs(Mv-Mv_anal)/Mv_anal
error_MBv = abs(MBv-MBv_anal)/Mv_anal
error_Mt = abs(Mt-Mt_anal)/Mv_anal
error_MBt = abs(MBt-MBt_anal)/Mv_anal
PRINT error_ur+error_sigma1+error_sigma3+error_Mv+error_MBv+error_Mt+error_MBt
|
