# Plastic Compression of a Pipe Assembly#

Problem Description:
• Two coaxial tubes, the inner one of 1020 CR steel and cross-sectional area $$A_{\mathrm{s}}$$, and the outer one of 2024-T4 aluminum alloy and of area $$A_{\mathrm{a}}$$, are compressed between heavy, flat end plates, as shown below. Determine the load-deflection curve of the assembly as it is compressed into the plastic region by an axial displacement. Assume that the end plates are so stiff that both tubes are shortened by exactly the same amount.

Reference:
• S. H. Crandall, N. C. Dahl, An Introduction to the Mechanics of Solids, McGraw-Hill Book Co., Inc., New York, NY, 1959, pg. 180, ex. 5.1.

Analysis Type(s):
• Static, Plastic Analysis (ANTYPE=0)

Element Type(s):
• Plastic Straight Pipe Element (PIPE288)

• 4-Node Finite Strain Shell (SHELL181)

• 3-D Structural Solid Elements (SOLID185)

Material Properties
• $$E_{\mathrm{s}} = 26875000\,psi$$

• $$\sigma_{\mathrm{(yp)s}} = 86000\,psi$$

• $$E_{\mathrm{a}} = 11000000\,psi$$

• $$\sigma_{\mathrm{(yp)a}} = 55000\,psi$$

• $$\nu = 0.3$$

Geometric Properties:
• $$l = 10\,in$$

• $$A_{\mathrm{s}} = 7\,in^2$$

• $$A_{\mathrm{a}} = 12\,in^2$$

• 1st Load Step: $$\delta = 0.032\,in$$

• 2nd Load Step: $$\delta = 0.050\,in$$

• 3rd Load Step: $$\delta = 0.100\,in$$

Analysis Assumptions and Modeling Notes:
• The following tube dimensions, which provide the desired cross-sectional areas, are arbitrarily chosen:

• Inner (steel) tube: inside radius = 1.9781692 in., wall thickness = 0.5 in.

• Outer (aluminum) tube: inside radius = 3.5697185 in., wall thickness = 0.5 in.

• The problem can be solved in three ways:

• using PIPE288 - the plastic straight pipe element

• using SOLID185 - the 3-D structural solid element

• using SHELL181 - the 4-Node Finite Strain Shell

• In the SOLID185 and SHELL181 cases, since the problem is axisymmetric, only a one element $$\theta$$ -sector is modeled. A small angle $$\theta = 6°$$ is arbitrarily chosen to reasonably approximate the circular boundary with straight sided elements. The nodes at the boundaries have the UX (radial) degree of freedom coupled. In the SHELL181 model, the nodes at the boundaries additionally have the ROTY degree of freedom coupled.

## Start MAPDL#

Start MAPDL and import Numpy and Pandas libraries.

# sphinx_gallery_thumbnail_path = '_static/vm7_setup.png'

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

from ansys.mapdl.core import launch_mapdl

# Start MAPDL.
mapdl = launch_mapdl()


## Pre-Processing#

Enter verification example mode and the pre-processing routine.

mapdl.clear()
mapdl.verify()
_ = mapdl.prep7()


## Parameterization#

# Angle of the model sector.
theta = 6

defl_ls1 = -0.032
defl_ls2 = -0.05
defl_ls3 = -0.1


## Define Element Type#

Set up the element types .

# Element type PIPE288.
mapdl.et(1, "PIPE288")

# Special Features are defined by keyoptions of pipe element.
# KEYOPT(4)(2)
# Hoop strain treatment:
# Thick pipe theory.
mapdl.keyopt(1, 4, 2)  # Cubic shape function

# Element type SOLID185.
mapdl.et(2, "SOLID185")

# Element type SHELL181.
mapdl.et(3, "SHELL181")  # FULL INTEGRATION

# Special Features are defined by keyoptions of shell element.
# KEYOPT(3)(2)
# Integration option:
# Full integration with incompatible modes.
mapdl.keyopt(3, 3, 2)

# Print
print(mapdl.etlist())


Out:

ELEMENT TYPE        1 IS PIPE288      3-D 2-NODE PIPE
KEYOPT( 1- 6)=        0      0      0        2      0      0
KEYOPT( 7-12)=        0      0      0        0      0      0
KEYOPT(13-18)=        0      0      0        0      0      0

ELEMENT TYPE        2 IS SOLID185     3-D 8-NODE STRUCTURAL SOLID
KEYOPT( 1- 6)=        0      0      0        0      0      0
KEYOPT( 7-12)=        0      0      0        0      0      0
KEYOPT(13-18)=        0      0      0        0      0      0

ELEMENT TYPE        3 IS SHELL181     4-NODE SHELL
KEYOPT( 1- 6)=        0      0      2        0      0      0
KEYOPT( 7-12)=        0      0      0        0      0      0
KEYOPT(13-18)=        0      0      0        0      0      0

CURRENT NODAL DOF SET IS  UX    UY    UZ    ROTX  ROTY  ROTZ
THREE-DIMENSIONAL MODEL


## Define Material#

Set up the material properties.

• Young Modulus of steel is: $$E_{\mathrm{s}} = 26875000\,psi$$,

• Yield strength of steel is: $$\sigma_{\mathrm{(yp)s}} = 86000\, psi$$,

• Young Modulus of aluminum is: $$E_{\mathrm{a}} = 11000000\,psi$$,

• Yield strength of aluminum is: $$\sigma_{\mathrm{(yp)a}} = 55000\,psi$$,

• Poisson’s ratio is: $$\nu = 0.3$$

# Steel material model.
# Define Young's moulus and Poisson ratio for Steel.
mapdl.mp("EX", 1, 26.875e6)
mapdl.mp("PRXY", 1, 0.3)

# Define non-linear material properties for Steel.
mapdl.tb("BKIN", 1, 1)
mapdl.tbtemp(0)
mapdl.tbdata(1, 86000, 0)

# Aluminum material model.
# Define Young's moulus and Poisson ratio for Aluminum.
mapdl.mp("EX", 2, 11e6)
mapdl.mp("PRXY", 2, 0.3)

# Define non-linear material properties for Aluminum.
mapdl.tb("BKIN", 2, 1)
mapdl.tbtemp(0)
mapdl.tbdata(1, 55000, 0)

# Print
print(mapdl.mplist())


Out:

MATERIAL NUMBER        1

TEMP        EX
0.2687500E+08

TEMP        PRXY
0.3000000

MATERIAL NUMBER        2

TEMP        EX
0.1100000E+08

TEMP        PRXY
0.3000000


## Plot Stress - Strain Curve#

Use Matplotlib library to plot material model curves of steel and aluminum.

# Define stress - strain properties of the steel.
steel = {"stress_s": [0, 86000, 86000, 86000], "strain_s": [0, 0.032, 0.1, 0.2]}

# Define yielding strength point of the steel on the curve.
xp = steel["strain_s"][1]
yp = steel["stress_s"][1]

# Set up the settings of the steel curve.
plt.plot(
steel["strain_s"],
steel["stress_s"],
label="1020 CR STEEL",
linewidth=2,
color="steelblue",
linestyle="-",
marker="o",
)
plt.plot(xp, yp, marker="o")

# Annotation settings
plt.annotate(
r"${(\sigma_{yp})_s}$",
xy=(xp, yp),
xytext=(0.05, 75000),
arrowprops=dict(facecolor="steelblue", shrink=0.05),
)

# Define stress - strain properties of the aluminum.
aluminum = {"stress_a": [0, 55000, 55000, 55000], "strain_a": [0, 0.05, 0.1, 0.2]}

# Define yielding strength point of the Aluminum on the curve.
xp = aluminum["strain_a"][1]
yp = aluminum["stress_a"][1]

# Set up the settings of the aluminum curve.
plt.plot(
aluminum["strain_a"],
aluminum["stress_a"],
label="2024-T4 Aluminum",
linewidth=2,
color="sandybrown",
linestyle="-",
marker="o",
)
plt.plot(xp, yp, marker="o")

# Annotation settings
plt.annotate(
r"${(\sigma_{yp})_a}$",
xy=(xp, yp),
xytext=(0.07, 45000),
arrowprops=dict(facecolor="sandybrown", shrink=0.05),
)

plt.grid(True)
plt.legend()
plt.title("Stress - Strain Curve", fontsize=18)
plt.show()


## Define Section#

Set up the cross-section properties for a shell and pipe elements.

# Shell cross-section for inside tube(steel).
mapdl.sectype(1, "SHELL")

# Thickness (SHELL181)
mapdl.secdata(0.5, 1, 0, 5)

# Shell cross-section for outside tube(aluminum).
mapdl.sectype(2, "SHELL")

# Thickness (SHELL181)
mapdl.secdata(0.5, 2, 0, 5)

# Define Pipe cross-section for inside tube(steel).
mapdl.sectype(3, "PIPE")

# Outside diameter and wall thickness settings for inside tube(PIPE288).
mapdl.secdata(4.9563384, 0.5)

# Pipe cross-section for outside tube(aluminum) .
mapdl.sectype(4, "PIPE")

# Outside diameter and wall thickness settings for outside tube (PIPE288).
mapdl.secdata(8.139437, 0.5)

# Print the section properties for all sections.
print(mapdl.slist())


Out:

*****ANSYS VERIFICATION RUN ONLY*****
DO NOT USE RESULTS FOR PRODUCTION

SECTION ID NUMBER:           1
SHELL SECTION TYPE:
SHELL SECTION NAME IS:
SHELL SECTION DATA SUMMARY:
Number of Layers    =      1
Total Thickness     =     0.500000

Layer      Thickness   MatID   Ori. Angle  Num Intg. Pts

1        0.5000     1        0.0000     5

Shell Section is offset to MID surface of Shell

Section Solution Controls
User Transverse Shear Stiffness (11)=  0.0000
(22)=  0.0000
(12)=  0.0000
Added Mass Per Unit Area            =  0.0000
Hourglass Scale Factor; Membrane    =  1.0000
Bending     =  1.0000
Drill Stiffness Scale Factor        =  1.0000

SECTION ID NUMBER:           2
SHELL SECTION TYPE:
SHELL SECTION NAME IS:
SHELL SECTION DATA SUMMARY:
Number of Layers    =      1
Total Thickness     =     0.500000

Layer      Thickness   MatID   Ori. Angle  Num Intg. Pts

1        0.5000     2        0.0000     5

Shell Section is offset to MID surface of Shell

Section Solution Controls
User Transverse Shear Stiffness (11)=  0.0000
(22)=  0.0000
(12)=  0.0000
Added Mass Per Unit Area            =  0.0000
Hourglass Scale Factor; Membrane    =  1.0000
Bending     =  1.0000
Drill Stiffness Scale Factor        =  1.0000

SECTION ID NUMBER:           3
PIPE SECTION NAME IS:
PIPE SECTION DATA SUMMARY:
Outside Diameter     =  4.9563
Thickness            = 0.50000
Area                 =  6.9946
Iyy                  =  17.559
Torsion Constant     =  35.118
Shear Correction-yy  = 0.50995

SECTION ID NUMBER:           4
PIPE SECTION NAME IS:
PIPE SECTION DATA SUMMARY:
Outside Diameter     =  8.1394
Thickness            = 0.50000
Area                 =  11.991
Iyy                  =  87.735
Torsion Constant     =  175.47
Shear Correction-yy  = 0.50305


## Define Geometry#

Set up the nodes and create the elements through the nodes.

# Generate nodes and elements for PIPE288.
mapdl.n(1, x=0, y=0, z=0)
mapdl.n(2, x=0, y=0, z=10)

# Create element for steel(inside) tube cross-section.
mapdl.mat(1)
mapdl.secnum(3)
mapdl.e(1, 2)

# Create element for aluminum(outside) tube cross-section.
mapdl.mat(2)
mapdl.secnum(4)
mapdl.e(1, 2)

# Activate the global cylindrical coordinate system.
mapdl.csys(1)

# Generate nodes and elements for SOLID185.
mapdl.n(node=101, x=1.9781692)
mapdl.n(node=101, x=1.9781692)
mapdl.n(node=102, x=2.4781692)
mapdl.n(node=103, x=3.5697185)
mapdl.n(node=104, x=4.0697185)
mapdl.n(node=105, x=1.9781692, z=10)
mapdl.n(node=106, x=2.4781692, z=10)
mapdl.n(node=107, x=3.5697185, z=10)
mapdl.n(node=108, x=4.0697185, z=10)

# Generate 2nd set of nodes to form a theta degree slice.
mapdl.ngen(itime=2, inc=10, node1=101, node2=108, dy=theta)

# Rotate nodal coordinate systems into the active system.
mapdl.nrotat(node1=101, node2=118, ninc=1)

# Create elements for Inside (Steel) tube.
mapdl.type(2)
mapdl.mat(1)
mapdl.e(101, 102, 112, 111, 105, 106, 116, 115)

# Create elements for Outside (Aluminum) tube
mapdl.mat(2)
mapdl.e(103, 104, 114, 113, 107, 108, 118, 117)

# Generate nodes.
mapdl.n(node=201, x=2.2281692)
mapdl.n(node=203, x=2.2281692, z=10)
mapdl.n(node=202, x=3.8197185)
mapdl.n(node=204, x=3.8197185, z=10)

# Generate nodes to form a theta degree slice
mapdl.ngen(itime=2, inc=4, node1=201, node2=204, dy=theta)

# Create element for Steel(Inside) tube cross-section.
mapdl.type(3)
mapdl.secnum(1)
mapdl.e(203, 201, 205, 207)

# Create element for Aluminum(Outside) tube cross-section.
mapdl.secnum(2)
mapdl.e(204, 202, 206, 208)

# Plot element model to demonstrate the axisymmetric element model.
cpos = [
(19.67899462804619, 17.856836088414664, 22.644135378046194),
(2.03485925, 0.21270071036846988, 5.0),
(0.0, 0.0, 1.0),
]
_ = mapdl.eplot(cpos=cpos)


## Define Boundary Conditions#

Application of boundary conditions (BC) for simplified axisymmetric model.

# Apply constraints to the PIPE288 model.
# Fix all DOFs for bottom end of PIPE288.
mapdl.d(node=1, lab="ALL")

# Allow only UZ DOF at top end of the PIPE288.
mapdl.d(node=2, lab="UX", lab2="UY", lab3="ROTX", lab4="ROTY", lab5="ROTZ")

# Apply constraints to SOLID185 and SHELL181 models"
# Couple nodes at boundary in RADIAL direction for SOLID185.
mapdl.cp(nset=1, lab="UX", node1=101, node2=111, node3=105, node4=115)
mapdl.cpsgen(itime=4, nset1=1)

# Couple nodes at boundary in RADIAL direction for the SHELL181.
mapdl.cp(5, lab="UX", node1=201, node2=205, node3=203, node4=20)
mapdl.cpsgen(itime=2, nset1=5)

# Couple nodes at boundary in ROTY dir for SHELL181.
mapdl.cp(7, lab="ROTY", node1=201, node2=205)
mapdl.cpsgen(itime=4, nset1=7)

# Select only nodes in SOLID185 and SHELL181 models.
mapdl.nsel(type_="S", item="NODE", vmin=101, vmax=212)

# Select only nodes at theta = 0 from the selected set.
mapdl.nsel("R", "LOC", "Y", 0)

# Apply symmetry boundary conditions.
mapdl.dsym("SYMM", "Y", 1)

# Select only nodes in SOLID185 and SHELL181 models.
mapdl.nsel(type_="S", item="NODE", vmin=101, vmax=212)

# elect nodes at theta from the selected set.
mapdl.nsel("R", "LOC", "Y", theta)

# Apply symmetry boundary conditions.
mapdl.dsym("SYMM", "Y", 1)

# Select all nodes and RESELECT only nodes at Z = 0.
mapdl.nsel("ALL")
mapdl.nsel("R", "LOC", "Z", 0)

# Constrain bottom nodes in Z direction.
mapdl.d("ALL", "UZ", 0)

# Select all nodes.
mapdl.nsel("ALL")
_ = mapdl.finish()


## Solve#

Enter solution mode and solve the system.

# Start solution procedure.
mapdl.slashsolu()

# Define solution function.
def solution(deflect):
mapdl.nsel("R", "LOC", "Z", 10)
mapdl.d(node="ALL", lab="UZ", value=deflect)
mapdl.nsel("ALL")
mapdl.solve()

# Run each load step to reproduce needed deflection subsequently.
solution(deflect=defl_ls1)

solution(deflect=defl_ls2)

solution(deflect=defl_ls3)
_ = mapdl.finish()


## Post-processing#

Enter post-processing.

# Enter the post-processing routine.
_ = mapdl.post1()


Set up the function to get load values of each load step of the simplified axisymmetric model and convert it to the full model.

def getload():

# Select the nodes in the PIPE288 element model.
mapdl.nsel(type_="S", item="NODE", vmin=1, vmax=2)
mapdl.nsel("R", "LOC", "Z", 0)

# Sum the nodal force contributions of elements.
mapdl.fsum()

# Extrapolation of the force results in the full 360 (deg) model.
load_288 = mapdl.get_value("FSUM", 0, "ITEM", "FZ")

# Select the nodes in the SOLID185 element model.
mapdl.nsel(type_="S", item="NODE", vmin=101, vmax=118)
mapdl.nsel("R", "LOC", "Z", 0)
mapdl.fsum()

# Get the force value of the simplified model.
load_185_theta = mapdl.get_value("FSUM", 0, "ITEM", "FZ")

# Extrapolation of the force results in the full 360 (deg) model.

# Select the nodes in the SHELL181 element model.
mapdl.nsel("S", "NODE", "", 201, 212)
mapdl.nsel("R", "LOC", "Z", 0)

# Sum the nodal force contributions of elements.
mapdl.fsum()

# Get the force value of the simplified model.
load_181_theta = mapdl.get_value("FSUM", 0, "ITEM", "FZ")

# Extrapolation of the force results in the full 360 (deg) model.

# Return load results of each element model.


Obtain the loads of the model using getload() function.

# Activate Load Step 1 and extract load data.
mapdl.set(1, 1)

mapdl.set(2, 1)

mapdl.set(3, 1)


## Check Results#

Finally we have the results of the loads for the simplified axisymmetric model, which can be compared with expected target values for models with PIPE288, SOLID185, and SHELL181 elements. Loads expected for each load step are:

• 1st Load Step with deflection $$\delta = 0.032 (in)$$ has $$load_1 = 1024400\,(lb)$$.

• 2nd Load Step with deflection $$\delta = 0.05 (in)$$ has $$load_2 = 1262000\,(lb)$$.

• 3rd Load Step with deflection $$\delta = 0.1 (in)$$ has $$load_3 = 1262000\,(lb)$$.

target_res = np.asarray(
[1024400, 1262000, 1262000, 1024400, 1262000, 1262000, 1024400, 1262000, 1262000]
)

simulation_res = np.asarray(
[
pipe288_ls1,
pipe288_ls2,
pipe288_ls2,
solid185_ls1,
solid185_ls2,
solid185_ls3,
shell181_ls1,
shell181_ls2,
shell181_ls3,
]
)

main_columns = {
"Target": target_res,
"Mechanical APDL": simulation_res,
"Ratio": list(np.divide(simulation_res, target_res)),
}

row_tuple = [
("PIPE288", "Load, lb for Deflection = 0.032 in"),
("PIPE288", "Load, lb for Deflection = 0.05 in"),
("PIPE288", "Load, lb for Deflection = 0.1 in"),
("SOLID185", "Load, lb for Deflection = 0.032 in"),
("SOLID185", "Load, lb for Deflection = 0.05 in"),
("SOLID185", "Load, lb for Deflection = 0.1 in"),
("SHELL181", "Load, lb for Deflection = 0.032 in"),
("SHELL181", "Load, lb for Deflection = 0.05 in"),
("SHELL181", "Load, lb for Deflection = 0.1 in"),
]

index_names = ["Element Type", "Load Step"]
row_indexing = pd.MultiIndex.from_tuples(row_tuple)
df = pd.DataFrame(main_columns, index=row_indexing)

df.style.set_caption("Results Comparison",).set_table_styles(
[
{
"props": [
("background-color", "#FFEFD5"),
("color", "black"),
("border", "0.5px solid black"),
("font-style", "italic"),
("text-align", "center"),
],
},
{
"props": [
("background-color", "#FFEFD5"),
("color", "black"),
("border", "0.5px solid black"),
("font-style", "italic"),
("text-align", "center"),
],
},
{"selector": "td:hover", "props": [("background-color", "#FFF8DC")]},
{"selector": "th", "props": [("max-width", "120px")]},
{"selector": "", "props": [("border", "0.5px solid black")]},
{
"selector": "caption",
"props": [
("color", "black"),
("font-style", "italic"),
("font-size", "24px"),
("text-align", "center"),
],
},
],
).set_properties(
**{
"background-color": "#FFFAFA",
"color": "black",
"border-color": "black",
"border-width": "0.5px",
"border-style": "solid",
"text-align": "center",
}
).format(
"{:.3f}"
)

Results Comparison
Target Mechanical APDL Ratio
PIPE288 Load, lb for Deflection = 0.032 in 1024400.000 1024400.000 1.000
Load, lb for Deflection = 0.05 in 1262000.000 1262000.000 1.000
Load, lb for Deflection = 0.1 in 1262000.000 1262000.000 1.000
SOLID185 Load, lb for Deflection = 0.032 in 1024400.000 1022529.000 0.998
Load, lb for Deflection = 0.05 in 1262000.000 1259695.000 0.998
Load, lb for Deflection = 0.1 in 1262000.000 1259695.000 0.998
SHELL181 Load, lb for Deflection = 0.032 in 1024400.000 1023932.000 1.000
Load, lb for Deflection = 0.05 in 1262000.000 1261654.000 1.000
Load, lb for Deflection = 0.1 in 1262000.000 1261423.000 1.000

stop mapdl

mapdl.exit()


Total running time of the script: ( 0 minutes 0.872 seconds)

Gallery generated by Sphinx-Gallery