|
| 1 | +""" |
| 2 | +.. _pressure_vessel_example: |
| 3 | +
|
| 4 | +Pressure Vessel |
| 5 | +--------------- |
| 6 | +This example demonstrates how to create a basic pressure vessel and |
| 7 | +apply a pressure to it. |
| 8 | +
|
| 9 | +Also shown here: |
| 10 | +- Various ways of accessing stress results from MAPDL. |
| 11 | +- Comparison between PRNSOL, *VGET (efficient wrapping), and the legacy reader. |
| 12 | +- Notes regarding FULL vs. POWER graphics when using PRNSOL. |
| 13 | +
|
| 14 | +""" |
| 15 | +import numpy as np |
| 16 | +from ansys.mapdl.core import launch_mapdl |
| 17 | + |
| 18 | +# start mapdl, enter the preprocessor, and set the units |
| 19 | +mapdl = launch_mapdl() |
| 20 | + |
| 21 | +mapdl.clear() |
| 22 | +mapdl.prep7() |
| 23 | + |
| 24 | +# US Customary system using inches (in, lbf*s2/in, s, °F). |
| 25 | +mapdl.units('BIN') |
| 26 | + |
| 27 | + |
| 28 | +############################################################################### |
| 29 | +# Set the materials and element type |
| 30 | + |
| 31 | +mapdl.et(1, "SOLID285") |
| 32 | +mapdl.mp("EX", 1, 10e6) |
| 33 | +mapdl.mp("PRXY", 1, 0.3) |
| 34 | +mapdl.mp("DENS", 1, 0.1) |
| 35 | +print(mapdl.mplist()) |
| 36 | + |
| 37 | + |
| 38 | +############################################################################### |
| 39 | +# Create the Geometry |
| 40 | + |
| 41 | +# area generation |
| 42 | +height = 10 |
| 43 | +inner_width = 2.5 |
| 44 | +outer_width = 3.5 |
| 45 | +mapdl.rectng(inner_width, outer_width, 0, height) |
| 46 | +mapdl.cyl4(0, height, inner_width, 0, outer_width, 90) |
| 47 | + |
| 48 | +# combine areas |
| 49 | +a_comb = mapdl.aadd(1, 2) |
| 50 | +mapdl.aplot(color='grey', background='w', show_area_numbering=True) |
| 51 | + |
| 52 | +# Generate a cylindrical volume by rotating an area pattern about an axis |
| 53 | +mapdl.vrotat(a_comb, pax1=6, arc=90) |
| 54 | +mapdl.vplot(color='grey', background='w', show_area_numbering=True, cpos='zy') |
| 55 | + |
| 56 | + |
| 57 | +############################################################################### |
| 58 | +# Create the mesh |
| 59 | +mapdl.smrtsize(1) |
| 60 | +mapdl.esize(0.25, 0) |
| 61 | +mapdl.mshape(1, "3D") |
| 62 | +mapdl.mshkey(0) |
| 63 | +mapdl.vmesh("ALL") |
| 64 | +mapdl.eplot(color='grey', background='w') |
| 65 | + |
| 66 | + |
| 67 | +############################################################################### |
| 68 | +# Solve |
| 69 | + |
| 70 | +# boundary condition selection |
| 71 | +mapdl.geometry.area_select([3, 5, 7]) |
| 72 | +mapdl.da("ALL", "SYMM") |
| 73 | +mapdl.allsel() |
| 74 | + |
| 75 | +# apply pressure |
| 76 | +mapdl.geometry.area_select([1, 6]) |
| 77 | +mapdl.sfa("ALL", 1, "PRES", 1000) |
| 78 | +mapdl.allsel() |
| 79 | + |
| 80 | +# solver |
| 81 | +mapdl.run("/SOL") |
| 82 | +mapdl.antype(0) |
| 83 | +mapdl.outres("ALL", "ALL") |
| 84 | +mapdl.run("/STATUS,SOLU") |
| 85 | +sol_output = mapdl.solve() |
| 86 | +mapdl.finish() |
| 87 | + |
| 88 | + |
| 89 | +############################################################################### |
| 90 | +# Post-Processing |
| 91 | +# ~~~~~~~~~~~~~~~ |
| 92 | +# Enter the MAPDL post-postprocessing routine (/POST1) and obtain the |
| 93 | +# von-mises stress for the single static solution. Here, we use MAPDL |
| 94 | +# directly to obtain the results using a wrapper around the VGET |
| 95 | +# method to efficiently obtain results without outputting to disk. |
| 96 | + |
| 97 | +# enter the postprocessing routine |
| 98 | +mapdl.post1() |
| 99 | +mapdl.set(1, 1) |
| 100 | + |
| 101 | +# results directly from MAPDL's VGET command |
| 102 | +# *VGET, __VAR__, NODE, , S, EQV |
| 103 | +nnum = mapdl.mesh.nnum |
| 104 | +von_mises_mapdl = mapdl.post_processing.nodal_eqv_stress |
| 105 | + |
| 106 | +# we could print out the solution for each node with: |
| 107 | + |
| 108 | +print(f'\nNode Stress (psi)') |
| 109 | +for node_num, stress_value in zip(nnum[:5], von_mises_mapdl[:5]): |
| 110 | + print(f'{node_num:<5d} {stress_value:.3f}') |
| 111 | +print('...') |
| 112 | + |
| 113 | +# or simply get the maximum stress value and corresponding node |
| 114 | +idx = np.argmax(von_mises_mapdl) |
| 115 | +node_num = nnum[idx] |
| 116 | +stress_value = von_mises_mapdl[idx] |
| 117 | +print(f'\nMaximum Stress') |
| 118 | +print(f'Node Stress (psi)') |
| 119 | +print(f'{node_num:<5d} {stress_value:.3f}') |
| 120 | + |
| 121 | +############################################################################### |
| 122 | +# Plot the results |
| 123 | + |
| 124 | +mapdl.post_processing.plot_nodal_eqv_stress(cpos='zy') |
| 125 | + |
| 126 | + |
| 127 | +############################################################################### |
| 128 | +# We could, alternatively, get the exact same results by directly |
| 129 | +# accessing the result file using the legacy file reader |
| 130 | +# `ansys-mapdl-reader <https://github.com/pyansys/pymapdl-reader>`_. |
| 131 | + |
| 132 | +# access the result |
| 133 | +result = mapdl.result |
| 134 | + |
| 135 | +# Get the von mises stess and show that this is equivalent to the |
| 136 | +# stress obtained from MAPDL. |
| 137 | +nnum, stress = result.principal_nodal_stress(0) |
| 138 | +von_mises = stress[:, -1] # von-Mises stress is the right most column |
| 139 | +min_von_mises, max_von_mises = np.min(von_mises), np.max(von_mises) |
| 140 | +print('All close:', np.allclose(von_mises, von_mises_mapdl)) |
| 141 | + |
| 142 | +############################################################################### |
| 143 | +# That these results are equivalent to results from PRNSOL. |
| 144 | +# |
| 145 | +# .. note:: |
| 146 | +# Enabling POWER GRAPHICS with ``mapdl.graphics('POWER') will |
| 147 | +# change the averaging scheme. |
| 148 | + |
| 149 | +mapdl.header('OFF', 'OFF', 'OFF', 'OFF', 'OFF', 'OFF') |
| 150 | +table = mapdl.prnsol('S', 'PRIN').splitlines()[1:] |
| 151 | +prnsol_eqv = np.genfromtxt(table)[:, -1] # eqv is the last column |
| 152 | + |
| 153 | +# show these are equivalent (RTOL due to rounding within PRNSOL) |
| 154 | +print('All close:', np.allclose(von_mises, prnsol_eqv, rtol=1E-4)) |
| 155 | + |
| 156 | +print(f'LEGACY Reader and MAPDL VGET Min: {min_von_mises}') |
| 157 | +print(f'PRNSOL MAPDL Min: {prnsol_eqv.min()}') |
| 158 | +print() |
| 159 | +print(f'LEGACY Reader and MAPDL VGET Min: {max_von_mises}') |
| 160 | +print(f'PRNSOL MAPDL Min: {prnsol_eqv.max()}') |
0 commit comments