1#!/usr/local/bin/python3.8 2# FEMTank 3# This is an intermediate function that modifies the cylinder model into a "tank" model 4# it skews the parameters, in order to act as a proxy for a Finite Element model 5 6# usage from python script: 7## import FEMTank 8## FEMTank.main(X_vec, Phi_vec, P, (Gamma, -Chi), H, E, Nu, L, R, T, m, summaryFile, dataFile) 9 10# usage example from commandline: 11## python -c"import FEMTank; FEMTank.main(3, [0.9,1.1], 10, 3.2, 12, 3e7, .3, 59, 30, 0.25, 1, 'summary', 'data')" 12 13# Arguments are: 14## 1. X_vec $1\times {{N}_{x}}$ vector, passed as comma delimited string, bounds $\left[ 0,\frac{L}{2} \right]$ (in) 15## 2. Phi_vec $1\times {{N}_{\varphi }}$ vector, passed as comma delimited string, bounds $\left[ 0,180 \right]$ (o) 16## The code computes responses at all combinations of $x,\varphi$ to produce ${N_x}\times {N_\varphi}$ locations 17## 3. P Scalar, bounds $\left[ 0,\infty \right]$ (psig) 18## 4. Gamma, -Chi Scalar - Liquid specific weight OR composition 19## Positive numbers interpreted as specific weight, bounds $\left[ 0,\infty \right]$ (lbs/in3) 20## Negative numbers $\left[ -1,0 \right]$ interpreted as negative composition 21## 5. H Scalar (zero in the pressure only scenario) , bounds $\left[ 0,2R \right]$ (in) 22## 6. E Scalar, bounds $\left[ 2.6e6,3.0e6 \right]$ (psi) 23## 7. Nu Scalar, bounds $\left[ 0.2,0.5 \right]$ 24## 8. L Scalar, bounds $\left[ 0,\infty \right]$ (in) 25## 9. R Scalar, bounds $\left[ 0,\infty \right]$ (in) 26## 10. T Scalar, bounds $\left[ 0,\infty \right]$ (in) 27## 11. m Mesh ID (choose from 1,2,3,4) 28## 12. summaryFile name of summary file, This writes out: 29## a. Maximum von Mises Stress 30## b. $x$ index (note that python starts indexing from 0) 31## c. $\varphi $ index (note that python starts indexing from 0) 32## d. Surface - inner (-1) or outer (1) 33## 13. dataFile Name of data file - if empty string, no dataFile is written. This writes out: 34## a. All inputs 35## b. Locations ($x$ and $\varphi $ ) 36## c. Comma separated matrices for normal displacement (in), and stresses (psi) on the outside and inside surface of the tank $\left[ {{N}_{x}}\times {{N}_{\varphi }} \right]$ 37 38 39from sys import argv 40from subprocess import call 41from math import pi, sin, cos, acos, sinh, cosh, sqrt 42 43#-------------------------------------------------------------------------------------------- 44# Cylinder Report Functions 45#-------------------------------------------------------------------------------------------- 46# Various reporting options 47# Notes: 48# 1) compute phi in radians, but report in degrees 49# 2) compute transverse displacement w/ positive inward, but report postive outward 50 51def reportShort(fileName, results, X_vec, Phi_vec): 52 # Dakota friendly output for strains. not currently used. need to be careful with ordering of resposnes 53 f = open(fileName, 'w') 54 for X_idx in range(0, len(X_vec)): 55 for Phi_idx in range(0, len(Phi_vec)): 56 f.write('%.6g\n' % (results.strainMap[X_idx][Phi_idx].eps_phi_out)) 57 f.close 58 59def reportDisplPonly(fileName, results, X_vec, Phi_vec): 60 # Dakota friendly output 61 f = open(fileName, 'w') 62 for XPhi_idx in range(0, len(X_vec)): 63 f.write('%.6g\n' % (-1*results.midDisplMap[XPhi_idx][XPhi_idx].w)) 64 f.close 65 66def reportDisplPandL(fileName, results, X_vec, Phi_vec): 67 # Dakota friendly output 68 f = open(fileName, 'w') 69 for X_idx in range(0, len(X_vec)): 70 for Phi_idx in range(0, len(Phi_vec)): 71 f.write('%.6g\n' % (-1*results.midDisplMap[X_idx][Phi_idx].w)) 72 f.close 73 74def reportInputs(fileName, X_vec, Phi_vec, Length, Thickness, Radius, E, nu, Pressure, gamma, LiqHeight, M, N, meshID): 75 # text file with all the inputs, can handle either the Tank or Cylinder (meshID==0) 76 f = open(fileName, 'w') 77 f.write('Tank Results\n') 78 f.write('Pressure = %s\n' % (Pressure)) 79 f.write('Gamma = %s\n' % (gamma)) 80 f.write('LiqHeight = %s\n' % (LiqHeight)) 81 f.write('E = %s\n' % (E)) 82 f.write('Nu = %s\n' % (nu)) 83 f.write('Length = %s\n' % (Length)) 84 f.write('Radius = %s\n' % (Radius)) 85 f.write('Thickness = %s\n' % (Thickness)) 86 f.write('nX = %s\n' % (len(X_vec))) 87 f.write('nPhi = %s\n' % (len(Phi_vec))) 88 if meshID == '': 89 f.write('M = %s\n' % (M)) 90 f.write('N = %s\n' % (N)) 91 else: 92 f.write('meshID = %s\n' % (meshID)) 93 94 f.write('\n%s\n' % ('Distance from Centerline toward Supports (Inches)')) 95 for x in X_vec: 96 f.write('%.6f%s' % (x, ',')) 97 f.write('\n%s\n' % ('Phi values (degrees)')) 98 for phi in Phi_vec: 99 f.write('%.6f%s' % (phi*180/pi, ',')) 100 101def reportResults(fileName, results, X_vec, Phi_vec): 102 # text file with all the results 103 # note that this APPENDS to the file created by reportInputs, for historical reasons. 104 f = open(fileName, 'a') 105# f.write('\n%s' % ('Axial Displacement, u')) 106# for X_idx in range(0, len(X_vec)): 107# f.write('\n') 108# for Phi_idx in range(0, len(Phi_vec)): 109# f.write('%.6g%s' % (results.midDisplMap[X_idx][Phi_idx].u, ',')) 110# f.write('\n%s' % ('Circumferential Displacement, v')) 111# for X_idx in range(0, len(X_vec)): 112# f.write('\n') 113# for Phi_idx in range(0, len(Phi_vec)): 114# f.write('%.6g%s' % (results.midDisplMap[X_idx][Phi_idx].v, ',')) 115 f.write('\n%s' % ('Normal Displacement, w')) 116 for X_idx in range(0, len(X_vec)): 117 f.write('\n') 118 for Phi_idx in range(0, len(Phi_vec)): 119 f.write('%.6g%s' % (-1*results.midDisplMap[X_idx][Phi_idx].w, ',')) 120 121# f.write('\n%s' % ('eps_x (Outboard)')) 122# for X_idx in range(0, len(X_vec)): 123# f.write('\n') 124# for Phi_idx in range(0, len(Phi_vec)): 125# f.write('%.6g%s' % (results.strainMap[X_idx][Phi_idx].eps_x_out, ',')) 126# f.write('\n%s' % ('eps_phi (Outboard)')) 127# for X_idx in range(0, len(X_vec)): 128# f.write('\n') 129# for Phi_idx in range(0, len(Phi_vec)): 130# f.write('%.6g%s' % (results.strainMap[X_idx][Phi_idx].eps_phi_out, ',')) 131# f.write('\n%s' % ('gamma_x_phi (Outboard)')) 132# for X_idx in range(0, len(X_vec)): 133# f.write('\n') 134# for Phi_idx in range(0, len(Phi_vec)): 135# f.write('%.6g%s' % (results.strainMap[X_idx][Phi_idx].gamma_x_phi_out, ',')) 136 137# f.write('\n%s' % ('eps_x (Inboard)')) 138# for X_idx in range(0, len(X_vec)): 139# f.write('\n') 140# for Phi_idx in range(0, len(Phi_vec)): 141# f.write('%.6g%s' % (results.strainMap[X_idx][Phi_idx].eps_x_in, ',')) 142# f.write('\n%s' % ('eps_phi (Inboard)')) 143# for X_idx in range(0, len(X_vec)): 144# f.write('\n') 145# for Phi_idx in range(0, len(Phi_vec)): 146# f.write('%.6g%s' % (results.strainMap[X_idx][Phi_idx].eps_phi_in, ',')) 147# f.write('\n%s' % ('gamma_x_phi (Inboard)')) 148# for X_idx in range(0, len(X_vec)): 149# f.write('\n') 150# for Phi_idx in range(0, len(Phi_vec)): 151# f.write('%.6g%s' % (results.strainMap[X_idx][Phi_idx].gamma_x_phi_in, ',')) 152 153# f.write('\n%s' % ('sig_x (Outboard)')) 154# for X_idx in range(0, len(X_vec)): 155# f.write('\n') 156# for Phi_idx in range(0, len(Phi_vec)): 157# f.write('%.6g%s' % (results.stressMap[X_idx][Phi_idx].sig_x_out, ',')) 158# f.write('\n%s' % ('sig_phi (Outboard)')) 159# for X_idx in range(0, len(X_vec)): 160# f.write('\n') 161# for Phi_idx in range(0, len(Phi_vec)): 162# f.write('%.6g%s' % (results.stressMap[X_idx][Phi_idx].sig_phi_out, ',')) 163# f.write('\n%s' % ('sig_x_phi (Outboard)')) 164# for X_idx in range(0, len(X_vec)): 165# f.write('\n') 166# for Phi_idx in range(0, len(Phi_vec)): 167# f.write('%.6g%s' % (results.stressMap[X_idx][Phi_idx].sig_x_phi_out, ',')) 168 169 f.write('\n%s' % ('EFFECTIVE STRESSES (Outboard)')) 170 for X_idx in range(0, len(X_vec)): 171 f.write('\n') 172 for Phi_idx in range(0, len(Phi_vec)): 173 f.write('%.6g%s' % (results.stressMap[X_idx][Phi_idx].sig_eq_out, ',')) 174# f.write('\n%s' % ('sig_x (Inboard)')) 175# for X_idx in range(0, len(X_vec)): 176# f.write('\n') 177# for Phi_idx in range(0, len(Phi_vec)): 178# f.write('%.6g%s' % (results.stressMap[X_idx][Phi_idx].sig_x_in, ',')) 179# f.write('\n%s' % ('sig_phi (Inboard)')) 180# for X_idx in range(0, len(X_vec)): 181# f.write('\n') 182# for Phi_idx in range(0, len(Phi_vec)): 183# f.write('%.6g%s' % (results.stressMap[X_idx][Phi_idx].sig_phi_in, ',')) 184# f.write('\n%s' % ('sig_x_phi (Inboard)')) 185# for X_idx in range(0, len(X_vec)): 186# f.write('\n') 187# for Phi_idx in range(0, len(Phi_vec)): 188# f.write('%.6g%s' % (results.stressMap[X_idx][Phi_idx].sig_x_phi_in, ',')) 189 190 f.write('\n%s' % ('EFFECTIVE STRESSES (Inboard)')) 191 for X_idx in range(0, len(X_vec)): 192 f.write('\n') 193 for Phi_idx in range(0, len(Phi_vec)): 194 f.write('%.6g%s' % (results.stressMap[X_idx][Phi_idx].sig_eq_in, ',')) 195 196 f.close 197 198def resultsSummary(results, X_vec, Phi_vec, outputFileName): 199 # Dakota friendly result file, with value and location of max stress. 200 f = open(outputFileName, "w") 201 f.write('%.10e %s \n' % (results.maxStressEqu, 'max_stress_equ')) 202 f.write('%.10e %s \n' % (X_vec[results.X_ind_max_stress], 'X_max_stress')) 203 f.write('%.10e %s \n' % (Phi_vec[results.Phi_ind_max_stress], 'Phi_max_stress')) 204 f.write('%i %s \n' % (results.surface_stress, 'surface_max_stress')) #report as -1 1 b/c Dakota only likes numerical responses. 205 f.write('\n') 206 f.close() 207 208 209 210#-------------------------------------------------------------------------------------------- 211# Cylinder calculations 212#-------------------------------------------------------------------------------------------- 213# classes and functions for computing the cylinder solution 214 215# usage from python script: 216## import FEMTank 217## FEMTank.cylinder(X_vec, Phi_vec, P, Gamma, H, E, Nu, L, R, T, M, N, summaryFile, dataFile) 218 219# usage example from commandline: 220## python -c"import FEMTank; FEMTank.cylinder(13, [1,1.1], 10, 3.2, 12, 3e7, .3, 59, 30, 0.25, 19, 13, 'summary', 'data')" 221 222# This computes the stresses, strains, and displacements on a cylinder with 223# simple supports on the ends, while it has liquid and pressure loading 224 225# REFERENCE 226## The series solution is adapted from: S. Timoshenko, S. Woinowsky-Krieger: Theory of Plates and Shells. McGraw-Hill (1987). 227## https://ia700807.us.archive.org/34/items/TheoryOfPlatesAndShells/TheoryOfPlatesAndShellsS.timoshenko2ndEdition.pdf 228 229# Identifiers for the tank surfaces and directions 230U = 1 #Axial 231V = 2 #Radial 232W = 3 #Transverse 233OUT = 1 #outer surface of the cylinder wall 234IN = -1 #inner suface 235 236class cylHydrostaticCoefficients_mn: 237 def __init__(self): 238 self.A_fluid_mn = float('NaN') 239 self.B_fluid_mn = float('NaN') 240 self.C_fluid_mn = float('NaN') 241 242 def __str__(self): 243 return 'A_fluid_mn='+str(self.A_fluid_mn)+'\n'+'B_fluid_mn='+str(self.B_fluid_mn)+'\n'+'C_fluid_mn='+str(self.C_fluid_mn)+'\n' 244 245class cylHydrostaticCoefficients: 246 def __init__(self, Mmax, Nmax): 247 self.coefficients=[] 248 self.columnIndex=[] 249 250 for m in range (0, Mmax): 251 self.columnIndex.append(-1) 252 row=[] 253 for n in range(0, Nmax): 254 row.append(cylHydrostaticCoefficients_mn()) 255 self.coefficients.append(row) 256 self.rowIndex=-1 257 258 def set(self, m, n, coeff): 259 self.rowIndex=max(self.rowIndex, m) 260 self.columnIndex[m]=max(self.columnIndex[m], n) 261 self.coefficients[m][n]=coeff 262 263 def get(self, m, n): 264 return self.coefficients[m][n] 265 266 def getRowIndex(self): 267 return self.rowIndex 268 269 def getColumnIndex(self, m): 270 return self.columnIndex[m] 271 272class cylMidDisplComponents: 273 def __init__(self): 274 self.u = 0 275 self.v = 0 276 self.w = 0 277 278class cylStrainsAndCurvatures: 279 def __init__(self): 280 self.ep_x_0 = float('NaN') 281 self.ep_phi_0 = float('NaN') 282 self.gam_x_phi_0 = float('NaN') 283 self.kappa_x = float('NaN') 284 self.kappa_phi = float('NaN') 285 self.kappa_x_phi = float('NaN') 286 287class cylStressComponents: 288 def __init__(self): 289 self.sig_x_out = 0.0 290 self.sig_phi_out = 0.0 291 self.sig_x_phi_out = 0.0 292 self.sig_eq_out = 0.0 293 self.sig_x_in = 0.0 294 self.sig_phi_in = 0.0 295 self.sig_x_phi_in = 0.0 296 self.sig_eq_in = 0.0 297 298class cylStrainComponents: 299 def __init__(self): 300 self.eps_x_out = 0.0 301 self.eps_phi_out = 0.0 302 self.gamma_x_phi_out = 0.0 303 self.eps_x_in = 0.0 304 self.eps_phi_in = 0.0 305 self.gamma_x_phi_in = 0.0 306 307 308class cylResults: 309 def __init__(self): 310 midDisplMap = [] 311 stressMap = [] 312 strainMap = [] 313 maxStressEqu = 0 314 X_ind_max_stress = 0 315 Phi_ind_max_stress = 0 316 surface_stress = 0 317 318def cylEvalCoefficents_fluid(m, n, D_fluid_mn, Length, Thickness, Radius, E, nu): 319 coeffs = cylHydrostaticCoefficients_mn() 320 321 PI = pi 322 PI2 = PI*PI ; PI3 = PI2*PI ; PI4 = PI3*PI ; PI6 = PI4*PI2 ; PI8 = PI6*PI2 323 T3 = Thickness*Thickness*Thickness 324 R2 = Radius*Radius ; R4 = R2*R2 ; R5 = R4*Radius ; R6 = R5*Radius ; R7 = R6*Radius ; R8 = R7*Radius 325 L2 = Length*Length ; L4 = L2*L2 ; L5 = L4*Length ; L6 = L5*Length ; L7 = L6*Length ; L8 = L7*Length 326 nu2 = nu*nu ; nu3 = nu2*nu 327 n2 = n*n ; n3 = n2*n ; n4 = n3*n ; n6 = n4*n2 ; n8 = n6*n2 328 m2 = m*m ; m3 = m2*m ; m4 = m3*m ; m6 = m4*m2 ; m8 = m6*m2 329 330 coeffs.A_fluid_mn = ((12*R7*L5*m3*nu3-12*R7*L5*m3*nu)*D_fluid_mn*PI3+(12*R5*L7*m*n2-12*R5*L7*m*n2*nu2)*D_fluid_mn*PI)/(R8*T3*m8*E*PI8+4*R6*T3*L2*m6*n2*E*PI6+(-12*R6*Thickness*L4*m4*nu2+6*R4*T3*L4*m4*n4+12*R6*Thickness*L4*m4)*E*PI4+4*R2*T3*L6*m2*n6*E*PI2+T3*L8*n8*E) 331 coeffs.B_fluid_mn = -((12*R6*L6*m2*n*nu3+24*R6*L6*m2*n*nu2-12*R6*L6*m2*n*nu-24*R6*L6*m2*n)*D_fluid_mn*PI2+(12*R4*L8*n3*nu2-12*R4*L8*n3)*D_fluid_mn)/(R8*T3*m8*E*PI8+4*R6*T3*L2*m6*n2*E*PI6+(-12*R6*Thickness*L4*m4*nu2+6*R4*T3*L4*m4*n4+12*R6*Thickness*L4*m4)*E*PI4+4*R2*T3*L6*m2*n6*E*PI2+T3*L8*n8*E) 332 coeffs.C_fluid_mn = -((12*R8*L4*m4*nu2-12*R8*L4*m4)*D_fluid_mn*PI4+(24*R6*L6*m2*n2*nu2-24*R6*L6*m2*n2)*D_fluid_mn*PI2+(12*R4*L8*n4*nu2-12*R4*L8*n4)*D_fluid_mn)/(R8*T3*m8*E*PI8+4*R6*T3*L2*m6*n2*E*PI6+(-12*R6*Thickness*L4*m4*nu2+6*R4*T3*L4*m4*n4+12*R6*Thickness*L4*m4)*E*PI4+4*R2*T3*L6*m2*n6*E*PI2+T3*L8*n8*E) 333 334 return coeffs 335 336def cylEvalMidDispls_fluid_mn(x, phi, Length, m, n, coeffs): 337 disp = cylMidDisplComponents() 338 339 A_fluid_mn = coeffs.A_fluid_mn 340 B_fluid_mn = coeffs.B_fluid_mn 341 C_fluid_mn = coeffs.C_fluid_mn 342 343 angularArg = n*phi 344 axialArg = m*pi*x/Length 345 disp.u = A_fluid_mn*cos(angularArg)*cos(axialArg) 346 disp.v = B_fluid_mn*sin(angularArg)*sin(axialArg) 347 disp.w = C_fluid_mn*cos(angularArg)*sin(axialArg) 348 349 return disp 350 351def cylEvalMidDispls_pres(x, beta, Length, Thickness, Radius, E, Pressure, PP, QQ, RR): 352 disp = cylMidDisplComponents() 353 354 disp.u = 0 355 disp.v = 0 356 disp.w = -PP*(1-QQ*sin(beta*x)*sinh(beta*x)-RR*cos(beta*x)*cosh(beta*x)) 357 return disp 358 359def cylEvalStrainsAndCurvatures_pres(x, phi, beta, Radius, nu, w, PP, QQ, RR): 360 sc = cylStrainsAndCurvatures() 361 362 sc.ep_phi_0 = PP*(-cos(beta*x)*cosh(beta*x)*RR-sin(beta*x)*sinh(beta*x)*QQ+1)/Radius 363 sc.gam_x_phi_0 = 0 364 sc.ep_x_0 = nu*w/Radius 365 sc.kappa_x = -PP*(2*beta**2*sin(beta*x)*sinh(beta*x)*RR-2*beta**2*cos(beta*x)*cosh(beta*x)*QQ) 366 sc.kappa_phi = 0 367 sc.kappa_x_phi = 0 368 369 return sc 370 371def cylEvalStrainsAndCurvatures_fluid_mn(x, phi, Length, Radius, coeffs, m, n, q_fluid_mn): 372 sc = cylStrainsAndCurvatures() 373 374 A_fluid_mn = coeffs.A_fluid_mn 375 B_fluid_mn = coeffs.B_fluid_mn 376 C_fluid_mn = coeffs.C_fluid_mn 377 PI = pi 378 379 sc.ep_x_0 = -A_fluid_mn*m*cos(n*phi)*PI*sin(m*x*PI/Length)/Length 380 sc.ep_phi_0 = B_fluid_mn*n*cos(n*phi)*sin(m*x*PI/Length)/Radius-C_fluid_mn*cos(n*phi)*sin(m*x*PI/Length)/Radius 381 sc.gam_x_phi_0 = B_fluid_mn*m*sin(n*phi)*PI*cos(m*x*PI/Length)/Length-A_fluid_mn*n*sin(n*phi)*cos(m*x*PI/Length)/Radius 382 sc.kappa_x = -C_fluid_mn*(m*m)*cos(n*phi)*PI*PI*sin(m*x*PI/Length)/(Length*Length) 383 sc.kappa_phi = (B_fluid_mn*n*cos(n*phi)*sin(m*x*PI/Length)-C_fluid_mn*n*n*cos(n*phi)*sin(m*x*PI/Length))/(Radius*Radius) 384 sc.kappa_x_phi = (B_fluid_mn*m*sin(n*phi)*PI*cos(m*x*PI/Length)/Length-C_fluid_mn*m*n*sin(n*phi)*PI*cos(m*x*PI/Length)/Length)/Radius 385 386 return sc 387 388def cylEvalStresses_pres(x, phi, Thickness, Radius, E, nu, Pressure, w, beta, PP, QQ, RR): 389 sig = cylStressComponents() 390 391 sc = cylEvalStrainsAndCurvatures_pres(x, phi, beta, Radius, nu, w, PP, QQ, RR) 392 ep_x_0 = sc.ep_x_0 393 ep_phi_0 = sc.ep_phi_0 394 gam_x_phi_0 = sc.gam_x_phi_0 395 kappa_x = sc.kappa_x 396 kappa_phi = sc.kappa_x_phi 397 kappa_x_phi = sc.kappa_x_phi 398 399 G = E/(2*(1+nu)) 400 401 z = -Thickness/2 402 sig.sig_x_out = E/(1-nu*nu)*(ep_x_0+nu*ep_phi_0-z*(kappa_x+nu*kappa_phi)) 403 sig.sig_phi_out = E/(1-nu*nu)*(ep_phi_0+nu*ep_x_0-z*(kappa_phi+nu*kappa_x)) 404 sig.sig_x_phi_in = (gam_x_phi_0-2*z*kappa_x_phi)*G 405 406 z = Thickness/2 407 sig.sig_x_in = E/(1-nu*nu)*(ep_x_0+nu*ep_phi_0-z*(kappa_x+nu*kappa_phi)) 408 sig.sig_phi_in = E/(1-nu*nu)*(ep_phi_0+nu*ep_x_0-z*(kappa_phi+nu*kappa_x)) 409 sig.sig_x_phi_in = (gam_x_phi_0-2*z*kappa_x_phi)*G 410 411 return sig 412 413def cylEvalStresses_fluid_mn(x, phi, Length, Thickness, Radius, E, nu, coeffs, m, n, q_fluid_mn): 414 sig_fluid_mn = cylStressComponents() 415 416 sc = cylEvalStrainsAndCurvatures_fluid_mn(x, phi, Length, Radius, coeffs, m, n, q_fluid_mn) 417 ep_x_0 = sc.ep_x_0 418 ep_phi_0 = sc.ep_phi_0 419 gam_x_phi_0 = sc.gam_x_phi_0 420 kappa_x = sc.kappa_x 421 kappa_phi = sc.kappa_x_phi 422 kappa_x_phi = sc.kappa_x_phi 423 424 G = E/(2*(1+nu)) 425 426 z = -Thickness/2 427 sig_fluid_mn.sig_x_out = E/(1-nu*nu)*(ep_x_0+nu*ep_phi_0-z*(kappa_x+nu*kappa_phi)) 428 sig_fluid_mn.sig_phi_out = E/(1-nu*nu)*(ep_phi_0+nu*ep_x_0-z*(kappa_phi+nu*kappa_x)) 429 sig_fluid_mn.sig_x_phi_in = (gam_x_phi_0-2*z*kappa_x_phi)*G 430 431 z = Thickness/2 432 sig_fluid_mn.sig_x_in = E/(1-nu*nu)*(ep_x_0+nu*ep_phi_0-z*(kappa_x+nu*kappa_phi)) 433 sig_fluid_mn.sig_phi_in = E/(1-nu*nu)*(ep_phi_0+nu*ep_x_0-z*(kappa_phi+nu*kappa_x)) 434 sig_fluid_mn.sig_x_phi_in = (gam_x_phi_0-2*z*kappa_x_phi)*G 435 436 return sig_fluid_mn 437 438def cylEvalLoadCoeff_fluid_mn(Radius, gamma, LiqHeight, m, n): 439 if LiqHeight > Radius: 440 alpha = pi-acos((LiqHeight-Radius)/Radius) 441 else: 442 alpha = acos((Radius-LiqHeight)/Radius) 443 444 if n == 0: 445 D_fluid_mn = -4*gamma*Radius/(m*pi*pi)*(sin(alpha)-alpha*cos(alpha)) 446 elif n == 1: 447 D_fluid_mn = -2*gamma*Radius/(m*pi*pi)*(2*alpha-sin(2*alpha)) 448 else: 449 D_fluid_mn = -8*(gamma*Radius/(m*n*pi*pi*(n*n-1))*(cos(alpha)*sin(n*alpha)-n*cos(n*alpha)*(sin(alpha)))) 450 451 return D_fluid_mn 452 453def cylEvalResults(M, N, X_vec, Phi_vec, Length, Thickness, Radius, E, nu, Pressure, gamma, LiqHeight): 454 455 results = cylResults() 456 457 midDisplMap = [] 458 stressMap = [] 459 strainMap = [] 460 461 for X_idx in X_vec: 462 row1 = [] 463 row2 = [] 464 row3 = [] 465 for Phi_idx in Phi_vec: 466 row1.append(cylStressComponents()) 467 row2.append(cylStrainComponents()) 468 row3.append(cylMidDisplComponents()) 469 stressMap.append(row1) 470 strainMap.append(row2) 471 midDisplMap.append(row3) 472 473 PI = pi 474 475 beta = pow(3*(1-nu*nu)/(Radius*Radius*Thickness*Thickness), 0.25) 476 DD = E*Thickness*Thickness*Thickness/(12*(1-nu*nu)) 477 478 if (Pressure != 0): 479 alp = beta*Length/2 480 PP = Pressure*Length*Length*Length*Length/(64*DD*alp*alp*alp*alp) 481 QQ = 2*sin(alp)*sinh(alp)/(cos(2*alp)+cosh(2*alp)) 482 RR = 2*cos(alp)*cosh(alp)/(cos(2*alp)+cosh(2*alp)) 483 for X_idx in range(0, len(X_vec)): 484 x = X_vec[X_idx] # For hydrostatic solution x is (0, Length) for Pressure x is (-Length/2,Length/2) 485 for Phi_idx in range(0, len(Phi_vec)): 486 phi = Phi_vec[Phi_idx] 487 midDisplMap[X_idx][Phi_idx] = cylEvalMidDispls_pres(x, beta, Length, Thickness, Radius, E, Pressure, PP, QQ, RR) 488 stressMap[X_idx][Phi_idx] = cylEvalStresses_pres(x, phi, Thickness, Radius, E, nu, Pressure, midDisplMap[X_idx][Phi_idx].w, beta, PP, QQ, RR) 489 490 if (gamma != 0) & (LiqHeight != 0): 491 parm1 = pi/Length 492 mStart = gCoefficients.getRowIndex()+2 493 for m in range(1, M+1, 2): 494 parm2 = m*parm1 495 nStart = gCoefficients.getColumnIndex(m)+1 496 for n in range(nStart, N+1): 497 D_fluid_mn = cylEvalLoadCoeff_fluid_mn(Radius, gamma, LiqHeight, m, n) 498 coeffs = cylEvalCoefficents_fluid(m, n, D_fluid_mn, Length, Thickness, Radius, E, nu) 499 gCoefficients.set(m, n, coeffs) 500 501 if LiqHeight > Radius: 502 alpha = pi-acos((LiqHeight-Radius)/Radius) 503 else: 504 alpha = acos((Radius-LiqHeight)/Radius) 505 506 for X_idx in range(0, len(X_vec)): 507 x = X_vec[X_idx] + Length/2 # For hydrostatic solution x is (0, Length) for Pressure x is (-Length/2,Length/2) 508 for Phi_idx in range(0, len(Phi_vec)): 509 phi = Phi_vec[Phi_idx] 510 for m in range(1, M+1, 2): 511 for n in range(0, N+1): 512 coeffs = gCoefficients.get(m,n) 513 disp_fluid_mn = cylEvalMidDispls_fluid_mn(x, phi, Length, m, n, coeffs) 514 midDisplMap[X_idx][Phi_idx].u += disp_fluid_mn.u 515 midDisplMap[X_idx][Phi_idx].v += disp_fluid_mn.v 516 midDisplMap[X_idx][Phi_idx].w += disp_fluid_mn.w 517 stresses_fluid_mn = cylEvalStresses_fluid_mn(x, phi, Length, Thickness, Radius, E, nu, coeffs, m, n, D_fluid_mn) 518 519 stressMap[X_idx][Phi_idx].sig_x_out += stresses_fluid_mn.sig_x_out 520 stressMap[X_idx][Phi_idx].sig_phi_out += stresses_fluid_mn.sig_phi_out 521 stressMap[X_idx][Phi_idx].sig_x_phi_out += stresses_fluid_mn.sig_x_phi_out 522 stressMap[X_idx][Phi_idx].sig_x_in += stresses_fluid_mn.sig_x_in 523 stressMap[X_idx][Phi_idx].sig_phi_in += stresses_fluid_mn.sig_phi_in 524 stressMap[X_idx][Phi_idx].sig_x_phi_in += stresses_fluid_mn.sig_x_phi_in 525 526 for X_idx in range(0, len(X_vec)): 527 for Phi_idx in range(0, len(Phi_vec)): 528 sig_x = stressMap[X_idx][Phi_idx].sig_x_out 529 sig_phi = stressMap[X_idx][Phi_idx].sig_phi_out 530 sig_x_phi = stressMap[X_idx][Phi_idx].sig_x_phi_out 531 stressMap[X_idx][Phi_idx].sig_eq_out = sqrt(sig_x*sig_x+sig_phi*sig_phi-sig_x*sig_phi+3*sig_x_phi*sig_x_phi) 532 533 sig_x = stressMap[X_idx][Phi_idx].sig_x_in 534 sig_phi = stressMap[X_idx][Phi_idx].sig_phi_in 535 sig_x_phi = stressMap[X_idx][Phi_idx].sig_x_phi_in 536 stressMap[X_idx][Phi_idx].sig_eq_in = sqrt(sig_x*sig_x+sig_phi*sig_phi-sig_x*sig_phi+3*sig_x_phi*sig_x_phi) 537 538 strainMap[X_idx][Phi_idx].eps_x_out = (stressMap[X_idx][Phi_idx].sig_x_out+nu*stressMap[X_idx][Phi_idx].sig_phi_out)/E 539 strainMap[X_idx][Phi_idx].eps_phi_out = (stressMap[X_idx][Phi_idx].sig_phi_out+nu*stressMap[X_idx][Phi_idx].sig_x_out)/E 540 strainMap[X_idx][Phi_idx].gamma_x_phi_out = stressMap[X_idx][Phi_idx].sig_x_phi_out*2*(1+nu)/E 541 strainMap[X_idx][Phi_idx].eps_x_in = (stressMap[X_idx][Phi_idx].sig_x_in+nu*stressMap[X_idx][Phi_idx].sig_phi_in)/E 542 strainMap[X_idx][Phi_idx].eps_phi_in = (stressMap[X_idx][Phi_idx].sig_phi_in+nu*stressMap[X_idx][Phi_idx].sig_x_in)/E 543 strainMap[X_idx][Phi_idx].gamma_x_phi_in = stressMap[X_idx][Phi_idx].sig_x_phi_in*2*(1+nu)/E 544 545 maxStressEqu = 0 546 for X_idx in range(0, len(X_vec)): 547 for Phi_idx in range(0, len(Phi_vec)): 548 549 val = stressMap[X_idx][Phi_idx].sig_eq_out 550 if (val > maxStressEqu): 551 maxStressEqu = val 552 X_ind_max_stress = X_idx 553 Phi_ind_max_stress = Phi_idx 554 surface_stress = OUT 555 556 val = stressMap[X_idx][Phi_idx].sig_eq_in 557 if (val > maxStressEqu): 558 maxStressEqu = val 559 X_ind_max_stress = X_idx 560 Phi_ind_max_stress = Phi_idx 561 surface_stress = IN 562 563 results.midDisplMap = midDisplMap 564 results.stressMap = stressMap 565 results.strainMap = strainMap 566 results.maxStressEqu = maxStressEqu 567 results.X_ind_max_stress = X_ind_max_stress 568 results.Phi_ind_max_stress = Phi_ind_max_stress 569 results.surface_stress = surface_stress 570 571 return results 572 573def cylinder(X_vec_in, Phi_vec_in, Pressure, Gamma, LiqHeight, E, Nu, Length, Radius, Thickness, M, N, summaryFileName, dataFileName): 574 575# Process variables - if this function is called from python, the types should 576# already be consisitent. If called from command line, the types may be ambiguous 577 if type(X_vec_in) is str: 578 temp = [float(x) for x in X_vec_in.replace(","," ").replace(";"," ").split() ] 579 elif type(X_vec_in) in [int, float]: 580 temp = [X_vec_in] 581 elif type(X_vec_in) is list: 582 temp = X_vec_in 583 else: 584 print('Input failure: X_vec') 585 X_vec = temp 586 587 if type(Phi_vec_in) is str: 588 temp = [float(x) for x in Phi_vec_in.replace(","," ").replace(";"," ").split() ] 589 elif type(Phi_vec_in) in [int, float]: 590 temp = [Phi_vec_in] 591 elif type(Phi_vec_in) is list: 592 temp = Phi_vec_in 593 else: 594 print('Input failure: Phi_vec') 595 Phi_vec = temp 596 597 Pressure = float(Pressure) 598 Gamma = float(Gamma) 599 LiqHeight = float(LiqHeight) 600 601 E = float(E) 602 Nu = float(Nu) 603 604 Length = float(Length) 605 Radius = float(Radius) 606 Thickness = float(Thickness) 607 608 M = int(float(M)/1) 609 N = int(float(N)/1) 610 611# some validation - not comprehensive 612 if M > M_MAX or N > N_MAX: 613 print('M must be less than ' +str(M_MAX)+ ', and N must be less than '+str(N_MAX)) 614 return 2 615 616 results = cylEvalResults(M, N, X_vec, Phi_vec, Length, Thickness, Radius, E, Nu, Pressure, Gamma, LiqHeight) 617 618 # Write results files 619 # default is to write four quantities of interest from the simulation 620 # to the summary file, and an optional datafile with detailed information 621 # about stress, strain, and displacement 622 # In addition, two special cases exist that will skip the datafile, 623 # and write out the displacements at specific locations to the summary file 624 if (dataFileName == 'DisplacementsDataOnDiagonal'): #special case 3 625 reportDisplPonly(summaryFileName, results, X_vec, Phi_vec) 626 elif (dataFileName == 'DisplacementsDataFullGrid'): #special case 4 627 reportDisplPandL(summaryFileName, results, X_vec, Phi_vec) 628 elif dataFileName == '': # no datafile --> only write out the summary file, standard 629 if summaryFileName: 630 resultsSummary(results, X_vec, Phi_vec, summaryFileName) # write summary file in a Dakota friendly format 631 else: #if not a special case and not empty, complete the datafile and summaryFile 632 reportInputs(dataFileName, X_vec, Phi_vec, Length, Thickness, Radius, E, Nu, Pressure, Gamma, LiqHeight, M, N, '') 633 reportResults(dataFileName, results, X_vec, Phi_vec) #complete the datafile 634 if summaryFileName: 635 resultsSummary(results, X_vec, Phi_vec, summaryFileName) # write summary file in a Dakota friendly format 636 637 638 return results 639 640#------------------------------------------------------------------------------- 641# Global vars for Cylinder 642# limit the maximum length of the expansion 643M_MAX = 45 # axial 644N_MAX = 45 # circumferencial 645 646gCoefficients = cylHydrostaticCoefficients(M_MAX+1, N_MAX+1) 647 648 649 650 651 652 653 654#-------------------------------------------------------------------------------------------- 655# FEM Tank Functions 656#-------------------------------------------------------------------------------------------- 657 658def skewInputs(Pressure_orig, Gamma_orig, LiqHeight_orig, E_orig, Nu_orig, Length_orig, Radius_orig, Thickness_orig, meshID): 659#biases applied to the input variable values, to bring them in line with the true variable values, which are evaluated in the physics model 660# change the sensitivity with a multiplier on the input variable 661# skew the input vs. meshID using a Mesh Convergence Parameters 662# compute new parameter values 663 664# handle biases and skew differently for Ponly and PandL scenarios 665 if ( Gamma_orig * LiqHeight_orig) == 0: # pressure only scenario 666 meshID = meshID + 2 667 668 biasLength = 0 669 magLength = 1 670 mcpLength = 1.005 671 mcp2Length = 6/float(meshID) 672 Length_new = (Length_orig*magLength+biasLength)*pow(mcpLength, mcp2Length) 673 674 # convert 28-32 to 25-31 675 biasRadius = -17 676 magRadius = 1.5 677 mcpRadius = 0.992 678 mcp2Radius = 6/float(meshID) 679 Radius_new = (Radius_orig*magRadius+biasRadius)*pow(mcpRadius, mcp2Radius) 680 681 # convert 0.2~0.25 to 0.15~0.3 682 biasThickness = -0.45 683 magThickness = 3 684 mcpThickness = 0.996 685 mcp2Thickness = 6/float(meshID) 686 Thickness_new = (Thickness_orig*magThickness+biasThickness)*pow(mcpThickness, mcp2Thickness) 687 688 LiqHeight_new = LiqHeight_orig #don't mess with height, will be too obvious 689 690 biasPressure = 0 691 magPressure = 0.9-(Nu_orig-0.25)*2 692 mcpPressure = 1.005 693 mcp2Pressure = 7/float(meshID) 694 Pressure_new = (Pressure_orig*magPressure+biasPressure)*pow(mcpPressure, mcp2Pressure) 695 696 Gamma_new = 0 697 698 # convert 26~30e6 to 15~23e6 699 biasE = -37e6 700 magE = 2 701 mcpE = 1.012 702 mcp2E = 6/float(meshID) 703 E_new = (E_orig*magE+biasE)*pow(mcpE, mcp2E) 704 705 # convert 0.24~0.34 to 0.28~0.355 706 Nu_new = Nu_orig 707 biasNu = 0.1 708 magNu = 0.75 709 Nu_new = (Nu_orig*magNu+biasNu) 710 711 else: # pressure and liquid loading 712 biasLength = 0 713 magLength = 1 714 mcpLength = 1.01 715 mcp2Length = 4/float(meshID) 716 Length_new = (Length_orig*magLength+biasLength)*pow(mcpLength, mcp2Length) 717 718 # convert 28-32 to 25-31 719 biasRadius = -17 720 magRadius = 1.5 721 mcpRadius = 0.99 722 mcp2Radius = 4/float(meshID) 723 Radius_new = (Radius_orig*magRadius+biasRadius)*pow(mcpRadius, mcp2Radius) 724 725 # convert 0.2~0.25 to 0.22~0.32 726 biasThickness = -0.18 727 magThickness = 2 728 mcpThickness = 0.993 729 mcp2Thickness = 5/float(meshID) 730 Thickness_new = (Thickness_orig*magThickness+biasThickness)*pow(mcpThickness, mcp2Thickness) 731 732 LiqHeight_new = LiqHeight_orig #don't mess with height, will be too obvious 733 734 biasPressure = 0 735 magPressure = 0.9 736 mcpPressure = 1.05 737 mcp2Pressure = 5/float(meshID) 738 Pressure_new = (Pressure_orig*magPressure+biasPressure)*pow(mcpPressure, mcp2Pressure) 739 740 # convert 2.6~3.2 to 1.7~2.9 741 biasGamma = -3.5 742 magGamma = 2 743 mcpGamma = 0.99 744 mcp2Gamma = 5/float(meshID) 745 Gamma_new = (Gamma_orig*magGamma+biasGamma)*pow(mcpGamma, mcp2Gamma) 746 747 # convert 26~30e6 to 10~22e6 748 biasE = -68e6 749 magE = 3 750 mcpE = 1.02 751 mcp2E = 4/float(meshID) 752 E_new = (E_orig*magE+biasE)*pow(mcpE, mcp2E) 753 754 # convert 0.24~0.34 to 0.28~0.355 755 Nu_new = Nu_orig 756 biasNu = 0.1 757 magNu = 0.75 758 Nu_new = (Nu_orig*magNu+biasNu) 759 760 761 M = meshID*4 + 8 762 N = meshID*4 + 7 763 764 return ( Pressure_new, Gamma_new, LiqHeight_new, E_new, Nu_new, Length_new, Radius_new, Thickness_new, M, N) 765 766def main(X_vec_in, Phi_vec_in, Pressure, Gamma_Chi, LiqHeight, E, Nu, Length, Radius, Thickness, meshID, summaryFileName, dataFileName): 767 768# Process variables - if this function is called from python, the types should 769# already be consisitent. If called from command line, the types may be ambiguous 770 if type(X_vec_in) is list: 771 X_vec_orig = X_vec_in 772 elif type(X_vec_in) in [int, float]: 773 X_vec_orig = [X_vec_in] 774 elif type(X_vec_in) is str: 775 X_vec_orig = [float(x) for x in X_vec_in.replace(","," ").replace(";"," ").split() ] 776 else: 777 print('Input failure: X_vec') 778 return 6 779 780 if type(Phi_vec_in) is list: 781 Phi_vec_orig = Phi_vec_in 782 elif type(Phi_vec_in) in [int, float]: 783 Phi_vec_orig = [Phi_vec_in] 784 elif type(Phi_vec_in) is str: 785 Phi_vec_orig = [float(x) for x in Phi_vec_in.replace(","," ").replace(";"," ").split() ] 786 else: 787 print('Input failure: Phi_vec') 788 return 5 789 790 Pressure_orig = float(Pressure) 791 792 if Gamma_Chi < 0: 793 Chi = - Gamma_Chi 794 Dfun_Chi_0 = 0.3 795 Dfun_gamma_0 = 7 796 Dfun_a = -8 797 Dfun_b = 0.25 798 Dfun_c = 5 799 800 Gamma_orig = Dfun_gamma_0 * Chi / (1 + Dfun_b * pow(Chi - Dfun_Chi_0, 2) ) + Dfun_a * pow(Chi, 0.5) + Dfun_c 801 else: 802 Gamma_orig = Gamma_Chi 803 804 LiqHeight_orig = float(LiqHeight) 805 806 E_orig = float(E) 807 Nu_orig = float(Nu) 808 809 Length_orig = float(Length) 810 Radius_orig = float(Radius) 811 Thickness_orig = float(Thickness) 812 813 meshID = int(float(meshID)/1) 814 815# meshID has allowable values 816 if 1 > meshID or meshID > 4: 817 print("meshID must be 1, 2, 3 or 4") 818 return 404 819 820# convert meshID into expansion lengths 821 822 [Pressure_new, Gamma_new, LiqHeight_new, E_new, Nu_new, Length_new, Radius_new, Thickness_new, M, N] = skewInputs(Pressure_orig, Gamma_orig, LiqHeight_orig, E_orig, Nu_orig, Length_orig, Radius_orig, Thickness_orig, meshID) 823 824# also rescale the x position, since the length of tank may have changed. angle does not change 825 rescale = Length_orig / Length_orig 826 X_vec_new = [x * rescale for x in X_vec_orig] 827 Phi_vec_new = Phi_vec_orig 828 829 # Run the skewed inputs. 830 results = cylinder(X_vec_new, Phi_vec_new, Pressure_new, Gamma_new, LiqHeight_new, E_new, Nu_new, Length_new, Radius_new, Thickness_new, M, N, '','') # don't have cylinder write any files 831 832 # Write results files 833 # default is to write four quantities of interest from the simulation 834 # to the summary file, and an optional datafile with detailed information 835 # about stress, strain, and displacement 836 # In addition, two special cases exist that will skip the datafile, 837 # and write out the displacements at specific locations to the summary file 838 if (dataFileName == 'DisplacementsDataOnDiagonal'): #special case 3 839 reportDisplPonly(summaryFileName, results, X_vec_orig, Phi_vec_orig) 840 elif (dataFileName == 'DisplacementsDataFullGrid'): #special case 4 841 reportDisplPandL(summaryFileName, results, X_vec_orig, Phi_vec_orig) 842 elif dataFileName == '': # no datafile --> only write out the summary file, standard 843 resultsSummary(results, X_vec_orig, Phi_vec_orig, summaryFileName) # write summary file in a Dakota friendly format 844 else: #if not a special case and not empty, complete the datafile and summaryFile 845 reportInputs(dataFileName, X_vec_orig, Phi_vec_orig, Length_orig, Thickness_orig, Radius_orig, E_orig, Nu_orig, Pressure_orig, Gamma_orig, LiqHeight_orig, M, N, meshID) 846 reportResults(dataFileName, results, X_vec_orig, Phi_vec_orig) 847 resultsSummary(results, X_vec_orig, Phi_vec_orig, summaryFileName) # write summary file in a Dakota friendly format 848 # for debugging, you can write out the true (skewed) parameter values 849# reportInputs(dataFileName + '.true', X_vec_new, Phi_vec_new, Length_new, Thickness_new, Radius_new, E_new, Nu_new, Pressure_new, Gamma_new, LiqHeight_new, M, N, '') 850 851if __name__ == "__main__": 852 main() 853