Note
Go to the end to download the full example code
Stresses in a Long Cylinder#
- Problem description:
A long thick-walled cylinder is initially subjected to an internal pressure p. Determine the radial displacement \(\delta_r\) at the inner surface, the radial stress \(\sigma_r\) , and tangential stress \(\sigma_t\) , at the inner and outer surfaces and at the middle wall thickness. Internal pressure is then removed and the cylinder is subjected to a rotation ω about its center line. Determine the radial \(\sigma_r\) and tangential \(\sigma_t\) stresses at the inner wall and at an interior point located at r = Xi.
- Reference:
S. Timoshenko, Strength of Materials, Part II, Elementary Theory and Problems, 3rd Edition, D. Van Nostrand Co., Inc., New York, NY, 1956, pg. 213, problem 1 and pg. 213, article 42.
- Analysis type(s):
Static Analysis
ANTYPE=0
- Element type(s):
2-D 8-Node Structural Solid Elements (PLANE183)
- Material properties:
\(E = 30 \cdot 10^6 psi\)
\(\mu = 0.3\)
\(\rho = 0.00073 lb-sec^2/in^4\)
- Geometric properties:
\(a = 4 inches\)
\(b = 8 inches\)
\(X_i = 5.43 inches\)
- Loading:
\(p = 30,000 psi\)
\(\Omega = 1000 rad/sec\)
- Analysis Assumptions and Modeling Notes:
The axial length is arbitrarily selected. Elements are oriented such that surface stresses may be obtained at the inner and outer cylinder surfaces. POST1 is used to display linearized stresses through the thickness of the cylinder when it is subjected to an internal pressure.
# sphinx_gallery_thumbnail_path = '_static/vm25_setup.png'
# Importing the `launch_mapdl` function from the `ansys.mapdl.core` module
from ansys.mapdl.core import launch_mapdl
import numpy as np
import pandas as pd
# Launching MAPDL with specified settings
mapdl = launch_mapdl(loglevel="WARNING", print_com=True, remove_temp_dir_on_exit=True)
# Clearing the MAPDL database
mapdl.clear()
# Run the FINISH command to exists normally from a processor
mapdl.finish()
# Set the ANSYS version
mapdl.com("ANSYS MEDIA REL. 2022R2 (05/13/2022) REF. VERIF. MANUAL: REL. 2022R2")
# Run the /VERIFY command for VM25
mapdl.run("/VERIFY,VM25")
# Set the title of the analysis
mapdl.title("VM25 Stresses in a Long Cylinder")
# Enter the model creation /Prep7 preprocessor
mapdl.prep7()
# Deactivate automatic (smart) element sizing
mapdl.smrtsize(sizlvl="OFF")
DEACTIVATE SMART SIZING
Define element type and properties#
Use 2-D 8-Node or 6-Node Structural Solid and specify Axisymmetric element behavior via setting Keyopt(3)=1.
mapdl.et(1, "PLANE183", "", "", 1)
1
Define material#
Set up the material and its type (a single material), Young’s modulus of 30e6, Density, rho = 0.00073 and Poisson’s ratio NUXY of 0.3 is specified.
mapdl.mp("EX", 1, 30e6)
mapdl.mp("DENS", 1, 0.00073)
mapdl.mp("NUXY", 1, 0.3)
MATERIAL 1 NUXY = 0.3000000
Define geometry#
Set up the nodes and elements. This creates a mesh just like in the problem setup.
# Defining node 1 with coordinates (4)
mapdl.n(1, 4)
# Defining node 2 with coordinates (4+4/14)
mapdl.n(2, "4+4/14")
# Defining node 3 with coordinates (4+4/14, 0.5)
mapdl.n(3, "4+4/14", 0.5)
# Defining node 4 with coordinates (4, 0.5)
mapdl.n(4, 4, 0.5)
# Defining node 5 with coordinates (4+4/28)
mapdl.n(5, "4+4/28")
# Defining node 6 with coordinates (4+4/14, 0.25)
mapdl.n(6, "4+4/14", 0.25)
# Defining node 7 with coordinates (4+4/28, 0.5)
mapdl.n(7, "4+4/28", 0.5)
# Defining node 8 with coordinates (4, 0.25)
mapdl.n(8, 4, 0.25)
# Creating element 1 with nodes 2, 3, 4, 5, 6, 7, and 8
mapdl.e(1, 2, 3, 4, 5, 6, 7, 8)
# Generating additional nodes along element 14 using the specified
# parameters (from an existing pattern)
mapdl.egen(14, 8, 1, "", "", "", "", "", "", "", "4/14")
# Merging nodes based on their coordinates
mapdl.nummrg("NODE")
# Generate a mesh using EGEN command
mapdl.egen(2, 111, 1, 14, "", "", "", "", "", "", "", 0.5)
# Merge the nodes that share the same coordinates
mapdl.nummrg("NODE")
MERGE COINCIDENT NODES WITHIN TOLERANCE OF 0.10000E-03
NODE 4 USED FOR NODES 112
NODE 7 USED FOR NODES 116
NODE 3 USED FOR NODES 113
NODE 15 USED FOR NODES 124
NODE 11 USED FOR NODES 121
NODE 23 USED FOR NODES 132
NODE 19 USED FOR NODES 129
NODE 31 USED FOR NODES 140
NODE 27 USED FOR NODES 137
NODE 39 USED FOR NODES 148
NODE 35 USED FOR NODES 145
NODE 47 USED FOR NODES 156
NODE 43 USED FOR NODES 153
NODE 55 USED FOR NODES 164
NODE 51 USED FOR NODES 161
NODE 63 USED FOR NODES 172
NODE 59 USED FOR NODES 169
NODE 71 USED FOR NODES 180
NODE 67 USED FOR NODES 177
NODE 79 USED FOR NODES 188
NODE 75 USED FOR NODES 185
NODE 87 USED FOR NODES 196
NODE 83 USED FOR NODES 193
NODE 95 USED FOR NODES 204
NODE 91 USED FOR NODES 201
NODE 103 USED FOR NODES 212
NODE 99 USED FOR NODES 209
NODE 111 USED FOR NODES 220
NODE 107 USED FOR NODES 217
A total of 29 nodes were merged at 29 locations.
Define coupling and boundary conditions#
Apply a displacement boundary condition in the vertical direction (UY) to all nodes. Couple the axial displacements at the unconstrained Y-dir. Apply internal pressure of 30000 psi is applied on nodes. Also, apply dummy pressure for surface printout. Then exit prep7 processor.
mapdl.nsel("S", "LOC", "Y", 0) # Select nodes located on the Y-axis
# Apply a displacement boundary condition
mapdl.d("ALL", "UY")
mapdl.nsel("S", "LOC", "Y", 1) # Select nodes located on the positive Y-axis
# Couple the axial displacements at the unconstrained Y-dir
mapdl.cp(1, "UY", "ALL")
mapdl.nsel(
"S", "LOC", "X", 4
) # Select nodes located on the X-axis at a specific coordinate
# Apply internal pressure on nodes
mapdl.sf("", "PRES", 30000)
mapdl.nsel(
"S", "LOC", "X", 8
) # Select nodes located on the X-axis at a different coordinate
# Apply dummy pressure for surface printout
mapdl.sf("", "PRES", 1e-10)
# Selects all entities
mapdl.allsel()
# Element plot
mapdl.eplot(background="w")
# Finish the pre-processing processor
mapdl.finish()
# Save the finite element model
mapdl.save("MODEL")
ALL CURRENT MAPDL DATA WRITTEN TO FILE NAME=
FOR POSSIBLE RESUME FROM THIS POINT
Solve#
Enter solution mode and solve the system.
mapdl.slashsolu()
# Set the analysis type to STATIC
mapdl.antype("STATIC")
# Output results for all nodes
mapdl.outpr("", "ALL")
# Perform the solution for the load step 1, which is the internal pressure
mapdl.solve()
# exists solution processor for load case 1
mapdl.finish()
FINISH SOLUTION PROCESSING
***** ROUTINE COMPLETED ***** CP = 0.000
Post-processing#
Enter post-processing. Compute displacement and stress components.
mapdl.post1()
# Set the load step 1 and substep to 1
mapdl.set(1, 1)
USE LOAD STEP 1 SUBSTEP 1 FOR LOAD CASE 0
SET COMMAND GOT LOAD STEP= 1 SUBSTEP= 1 CUMULATIVE ITERATION= 1
TIME/FREQUENCY= 1.0000
TITLE= VM25 Stresses in a Long Cylinder
Inline functions in PyMAPDL to query node#
Retrieve nodal deflection and section stresses#
# Retrieves the displacement "DEF_4" of nodes associated to
# "LFT_NODE" in the X direction
def_4 = mapdl.get("DEF_4", "NODE", LFT_NODE, "U", "X")
# Retrieves the stress and store it "parm" of nodes associated to a
# variables 'LFT_NODE','MID_NODE' and 'RT_NODE'
rst_4_c1 = mapdl.get("RST_4_C1", "NODE", LFT_NODE, "S", "X")
rst_6_c1 = mapdl.get("RST_6_C1", "NODE", MID_NODE, "S", "X")
rst_8_c1 = mapdl.get("RST_8_C1", "NODE", RT_NODE, "S", "X")
tst_4_c1 = mapdl.get("TST_4_C1", "NODE", LFT_NODE, "S", "Z")
tst_6_c1 = mapdl.get("TST_6_C1", "NODE", MID_NODE, "S", "Z")
tst_8_c1 = mapdl.get("TST_8_C1", "NODE", RT_NODE, "S", "Z")
# Print the nodal stress solution (COMP means all stress components)
mapdl.prnsol("S", "COMP")
# Define a path with the name 'STRESS' and ID 2, no limits specified
mapdl.path("STRESS", 2, "", 48)
mapdl.ppath(1, LFT_NODE) # Define the path points using the variable 'LFT_NODE'
mapdl.ppath(2, RT_NODE) # Define the path points using the variable 'RT_NODE'
mapdl.plsect("S", "Z", -1) # Display the SZ stresses in a sectional plot
mapdl.plsect("S", "X", -1) # Display the SX stresses in a sectional plot
mapdl.prsect(-1) # Print linearized stresses
'PRINT LINEARIZED STRESS THROUGH A SECTION DEFINED BY PATH= STRESS DSYS= 0\n *****MAPDL VERIFICATION RUN ONLY*****\n DO NOT USE RESULTS FOR PRODUCTION\n\n ***** POST1 LINEARIZED STRESS LISTING *****\n INSIDE NODE = 1 OUTSIDE NODE = 106\n\n LOAD STEP 1 SUBSTEP= 1\n TIME= 1.0000 LOAD CASE= 0\n\n ** AXISYMMETRIC OPTION ** RHO = 0.47392E+12\n THE FOLLOWING X,Y,Z STRESSES ARE IN SECTION COORDINATES. \n\n ** MEMBRANE **\n SX SY SZ SXY SYZ SXZ\n -0.1000E+05 0.1255E-08 0.3000E+05 0.8005E-10 0.000 0.000 \n S1 S2 S3 SINT SEQV\n 0.3000E+05 0.000 -0.1000E+05 0.4000E+05 0.3606E+05\n\n ** BENDING ** I=INSIDE C=CENTER O=OUTSIDE\n SX SY SZ SXY SYZ SXZ\n I -0.1991E+05 0.1083E-08 0.1364E+05 0.000 0.000 0.000 \n C -4948. 0.1083E-09 0.1919E-07 0.000 0.000 0.000 \n O 0.1001E+05 -0.8666E-09 -0.1364E+05 0.000 0.000 0.000 \n S1 S2 S3 SINT SEQV\n I 0.1364E+05 0.000 -0.1991E+05 0.3355E+05 0.2922E+05\n C 0.000 0.000 -4948. 4948. 4948. \n O 0.1001E+05 0.000 -0.1364E+05 0.2365E+05 0.2056E+05\n\n ** MEMBRANE PLUS BENDING ** I=INSIDE C=CENTER O=OUTSIDE\n SX SY SZ SXY SYZ SXZ\n I -0.2991E+05 0.2338E-08 0.4364E+05 0.8005E-10 0.000 0.000 \n C -0.1495E+05 0.1363E-08 0.3000E+05 0.8005E-10 0.000 0.000 \n O 6.734 0.3879E-09 0.1636E+05 0.8005E-10 0.000 0.000 \n S1 S2 S3 SINT SEQV\n I 0.4364E+05 0.000 -0.2991E+05 0.7355E+05 0.6407E+05\n C 0.3000E+05 0.000 -0.1495E+05 0.4495E+05 0.3965E+05\n O 0.1636E+05 6.734 0.000 0.1636E+05 0.1636E+05\n\n ** PEAK ** I=INSIDE C=CENTER O=OUTSIDE\n SX SY SZ SXY SYZ SXZ\n I -0.3638E-11 0.2863E-08 6264. 0.1697E-07 0.000 0.000 \n C 7193. -0.9225E-09 -2245. 0.4677E-09 0.000 0.000 \n O 0.7061E-12 0.2176E-08 3632. 0.2536E-09 0.000 0.000 \n S1 S2 S3 SINT SEQV\n I 6264. 0.000 0.000 6264. 6264. \n C 7193. 0.000 -2245. 9438. 8540. \n O 3632. 0.000 0.000 3632. 3632. \n\n ** TOTAL ** I=INSIDE C=CENTER O=OUTSIDE \n SX SY SZ SXY SYZ SXZ\n I -0.2991E+05 0.5201E-08 0.4991E+05 0.1705E-07 0.000 0.000 \n C -7758. 0.4404E-09 0.2776E+05 0.5477E-09 0.000 0.000 \n O 6.734 0.2564E-08 0.1999E+05 0.3336E-09 0.000 0.000 \n S1 S2 S3 SINT SEQV TEMP\n I 0.4991E+05 0.000 -0.2991E+05 0.7982E+05 0.6984E+05 0.000 \n C 0.2776E+05 0.000 -7758. 0.3552E+05 0.3234E+05\n O 0.1999E+05 6.734 0.000 0.1999E+05 0.1999E+05 0.000'
Verify the results.#
# Set target values
target_def = 0.0078666
target_strss = [-30000, -7778, 0]
target_tst_strss = [50000, 27778, 20000]
# Fill result values
res_def = def_4
res_strss = [rst_4_c1, rst_6_c1, rst_8_c1]
res_tst_strss = [tst_4_c1, tst_6_c1, tst_8_c1]
title = f"""
------------------- VM25 RESULTS COMPARISON ---------------------
RESULTS FOR CASE p = 30,000 psi:
--------------------------------
"""
col_headers = ["TARGET", "Mechanical APDL", "RATIO"]
row_headers = ["Displacement, in (r = 4 in)"]
data = [
[target_def, res_def, abs(target_def / res_def)],
]
print(title)
print(pd.DataFrame(data, row_headers, col_headers))
# Radial stress results comparison
row_headers = [
"Stress_r, psi (r = 4 in)",
"Stress_r, psi (r = 6 in)",
"Stress_r, psi (r = 8 in)",
]
data = [target_strss, res_strss, np.abs(target_strss) / np.abs(res_strss)]
print(pd.DataFrame(np.transpose(data), row_headers, col_headers))
# Tangential stress results comparison
row_headers = [
"Stress_t, psi (r = 4 in)",
"Stress_t, psi (r = 6 in)",
"Stress_t, psi (r = 8 in)",
]
data = [
target_tst_strss,
res_tst_strss,
np.abs(target_tst_strss) / np.abs(res_tst_strss),
]
print(pd.DataFrame(np.transpose(data), row_headers, col_headers))
------------------- VM25 RESULTS COMPARISON ---------------------
RESULTS FOR CASE p = 30,000 psi:
--------------------------------
TARGET Mechanical APDL RATIO
Displacement, in (r = 4 in) 0.007867 0.007867 0.999992
TARGET Mechanical APDL RATIO
Stress_r, psi (r = 4 in) -30000.0 -29908.046900 1.003075
Stress_r, psi (r = 6 in) -7778.0 -7757.541500 1.002637
Stress_r, psi (r = 8 in) 0.0 6.734007 0.000000
TARGET Mechanical APDL RATIO
Stress_t, psi (r = 4 in) 50000.0 49908.0469 1.001842
Stress_t, psi (r = 6 in) 27778.0 27757.5410 1.000737
Stress_t, psi (r = 8 in) 20000.0 19993.2656 1.000337
Finish the post-processing processor.#
mapdl.finish()
EXIT THE MAPDL POST1 DATABASE PROCESSOR
***** ROUTINE COMPLETED ***** CP = 0.000
Set a new title for the analysis#
mapdl.title("VM25 Stresses in a Long Cylinder - Rotation About Axis")
# Resume the Finite Element (FE) "MODEL" save previously
mapdl.resume("MODEL")
RESUME ANSYS DATA FROM FILE NAME=
*** MAPDL GLOBAL STATUS ***
TITLE = VM25 Stresses in a Long Cylinder
1 ELEM TYPES DEFINED MAX ELEM TYPE NUMBER = 1
28 ELEMENTS DEFINED MAX ELEMENT NUMBER = 28
117 NODES DEFINED MAX NODE NUMBER = 222
1 MATERIALS DEFINED MAX MATERIAL NUMBER = 1
0 REAL CONSTS DEFINED MAX REAL CONST NUMBER = 0
0 SECTIONS DEFINED MAX SECTION NUMBER = 0
0 COORD SYS DEFINED MAX COORD SYS NUMBER = 0
1 COUPLING EQNS DEFINED MAX COUPLING EQN NUMBER = 1
ACTIVE COORDINATE SYSTEM = 0 (CARTESIAN)
NUMBER OF DEFINED NODAL CONSTRAINTS = 29
NUMBER OF DEFINED NODAL LOADS = 0
NUMBER OF DEFINED ELEM SURFACE LOADS = 4
NUMBER OF DEFINED ELEM BODY LOADS = 0
NUMBER OF DEFINED NODE BODY FORCES = 0
INITIAL JOBNAME =
CURRENT JOBNAME =
Solve#
Enter solution mode and solve the system.
mapdl.slashsolu()
# Print all results
mapdl.outpr("", "ALL")
PRINT BASI ITEMS WITH A FREQUENCY OF ALL
FOR ALL APPLICABLE ENTITIES
Define boundary conditions#
Apply a displacement boundary condition in the vertical direction (UY) to all nodes. Rotate the cylinder with an angular velocity of 1000 rad/sec. Also, apply dummy pressure for surface printout. Then exit solution processor.
mapdl.nsel(
"S", "LOC", "Y", 0
) # Select nodes located at Y=0 to prevent rigid body motion
mapdl.nsel("R", "LOC", "X", 4) # Select nodes located at X=4
# Displace all nodes in the Y-direction
mapdl.d("ALL", "UY")
mapdl.nsel("S", "LOC", "X", 4) # Select nodes located at X=4
# Apply a small pressure to allow stress printout
mapdl.sf("", "PRES", 1e-10)
mapdl.nsel("ALL") # Select all nodes
# Rotate the cylinder with an angular velocity of 1000 RAD/SEC
mapdl.omega("", 1000)
# Solve the problem in load step 2 - centrifugal loading
mapdl.solve()
# exists solution processor
mapdl.finish()
FINISH SOLUTION PROCESSING
***** ROUTINE COMPLETED ***** CP = 0.000
Post-processing#
Enter post-processing. Compute displacement and stress components.
mapdl.post1()
*****MAPDL VERIFICATION RUN ONLY*****
DO NOT USE RESULTS FOR PRODUCTION
***** MAPDL RESULTS INTERPRETATION (POST1) *****
USE LAST SUBSTEP ON RESULT FILE FOR LOAD CASE 0
SET COMMAND GOT LOAD STEP= 1 SUBSTEP= 1 CUMULATIVE ITERATION= 1
TIME/FREQUENCY= 1.0000
TITLE= VM25 Stresses in a Long Cylinder
Inline functions in PyMAPDL to query node#
Retrieve nodal deflection and section stresses#
Verify the results.#
# Set target values
target_strss = [0, 4753]
target_tst_strss = [40588, 29436]
# Fill result values
res_strss = [rst_4_c2, rst_x_c2]
res_tst_strss = [tst_4_c2, tst_x_c2]
title = f"""
RESULTS FOR CASE Rotation = 1000 rad/sec:
-----------------------------------------
"""
col_headers = ["TARGET", "Mechanical APDL", "RATIO"]
# Radial stress results comparison
row_headers = ["Stress_r, psi (r = 4 in)", "Stress_r, psi (r = 5.43 in)"]
data = [target_strss, res_strss, np.abs(target_strss) / np.abs(res_strss)]
print(title)
print(pd.DataFrame(np.transpose(data), row_headers, col_headers))
# Tangential stress results comparison
row_headers = ["Stress_t, psi (r = 4 in)", "Stress_t, psi (r = 5.43 in)"]
data = [
target_tst_strss,
res_tst_strss,
np.abs(target_tst_strss) / np.abs(res_tst_strss),
]
print(pd.DataFrame(np.transpose(data), row_headers, col_headers))
RESULTS FOR CASE Rotation = 1000 rad/sec:
-----------------------------------------
TARGET Mechanical APDL RATIO
Stress_r, psi (r = 4 in) 0.0 50.366142 0.000000
Stress_r, psi (r = 5.43 in) 4753.0 4957.113530 0.958824
TARGET Mechanical APDL RATIO
Stress_t, psi (r = 4 in) 40588.0 41671.1641 0.974007
Stress_t, psi (r = 5.43 in) 29436.0 29740.9404 0.989747
Finish the post-processing processor.#
mapdl.finish()
EXIT THE MAPDL POST1 DATABASE PROCESSOR
***** ROUTINE COMPLETED ***** CP = 0.000
Stop MAPDL.#
mapdl.exit()
Total running time of the script: (0 minutes 0.808 seconds)