# 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