COMSOL 6.4 - Compression of an Elastoplastic Pipe

In offshore applications, it is sometimes necessary to quickly seal a pipe as part of the prevention of a blowout. This example shows a simulation in which a circular pipe is squeezed by a tool. The main target of the analysis is to find the remaining cross-section area through which a fluid can pass.

The tutorial serves as an example of an analysis involving very large plastic strains and contact. Moreover, the computation is performed by two distinct formulations using either a multiplicative decomposition of the deformation gradient or an additive decomposition of the logarithmic strain.

Model Definition

The pipe has an external radius, R0, of 200 mm and a wall thickness of 25 mm. The pipe is compressed between the two sides of a tool that can be considered as rigid. The geometry is shown in Figure 1. The original position of the tool is 0.1 mm from the outer pipe wall. During the compression of the pipe, the width of the tool is decreased until the maximum available force, 6 MN, is reached. The tool is then retracted, while the pipe expands as an effect of elastic springback. Due to the symmetries, only one quarter of the geometry needs to be modeled. The problem is considered as 2D with the plane strain assumption.

Figure 1: The geometry of the pipe and tool.

The pipe is modeled using an elastoplastic material model and is assumed to be made of stainless steel with the following properties:

Table 1: Material data.
Property	Value
Young’s modulus	195 GPa
Poisson’s ratio	0.3
Yield stress	250 MPa
Ultimate tensile stress	616 MPa
Ultimate strain	0.52

The hardening curve is available as a text file containing pairs of data (plastic strain, stress), which can be imported as a function. The function is shown in Figure 2 below. This nonlinear hardening function σh(εpe) is defined in the node by an interpolation function. The stress is measured as true (Cauchy) stress, and the strain is measured as true (logarithmic) strain. The data can thus be used directly as input to COMSOL when large strain plasticity is used. Note that if the strain is above the ultimate strain, the curve is implicitly flat, so that deformation continues under constant stress.

Figure 2: The hardening curve.

In order to be able to assess the sensitivity to the material data, a multiplier to the hardening curve is introduced in the model, although it has the value 1 throughout this analysis.

Large plastic deformations can be taken into account in different ways. One option is to use a multiplicative decomposition of the deformation gradient

which allows a linear constitutive equation for the second Piola–Kirchhoff stress in the intermediate configuration

Here, the he elastic Green–Lagrange strain is

and the plastic deformations are thus characterized by the plastic part of the deformation gradient Fpl.

Another option is to use an additive decomposition of the logarithmic Hencky strain, for which the second Piola–Kirchhoff stress is defined as

Here, the logarithmic Hencky strain is

with F = RU = VR, and plastic deformations are characterized by the logarithmic plastic strain εH,pl. This latter formulation can be used as a less computationally demanding alternative when small elastic strains are expected, as in metal plasticity.

Results and Discussion

This example exhibits extremely large plastic strains. The deformation and stress state at the maximum compression are shown in Figure 3. The maximum stress displayed is above the ultimate tensile stress (616 MPa). This is caused by the extrapolation of the results from the integration points inside the elements, where the constitutive law is exactly fulfilled.

Figure 3: Distribution of von Mises stress at maximum compression.

The distribution of the equivalent plastic strain is shown in Figure 4 and Figure 5. As can be seen, the peak value above 1 is far above the ultimate strain (0.52). The majority of values above the ultimate strain are, however, located on the inside of the pipe, where the strain state is mainly in compression. At the outer edge of the pipe, the plastic strain approximately also exceeds the ultimate strain, but not to the same extent. There is thus a risk that cracks could start forming. The values of ultimate stress and strain are related to the specimen used for the testing (usually a cylindrical bar) and cannot directly be transferred to general multiaxial conditions.

Figure 4: Equivalent plastic strain at the end of the process.

Figure 5: Equivalent plastic strain at the end of the process, detail (multiplicative decomposition). The red contour indicates the ultimate strain (0.52).

The final state after the retraction of the compression tool is shown in Figure 6. There is some reversed yielding during the unloading process. The springback effect is very small.

Figure 6: Deformed shape and residual stresses at the end of the process (multiplicative decomposition).

The load used to compress the pipe is computed from the reaction force on the tool, and it is shown in Figure 7. When the inside of the pipe makes contact with itself, the force rises steeply.

Figure 7: Applied force as function of the compression.

The change in the cross section area during the simulation is shown in Figure 8. The original area is 96211 mm2. At the maximum compression, it is reduced to 4811 mm2, and in the final state it is 6060 mm2. Thus, the area increases by about 25% due to elastic springback. Still, the final cross section area is only 6% of the original. Due to the change in shape, the flow will be reduced even more, since the wall perimeter is essentially unchanged. The resistance to the flow is largest near the walls.

Figure 8: Development of the cross section area as function of compression.

Notes About the COMSOL Implementation

•

In order to capture when the maximum force from the indenter is reached, a state variable is used. The state variable is updated after each parameter step, and has an update expression that stores the parameter value when the maximum attainable force is reached.

•

Since symmetry is used, the possible self-contact inside the pipe is modeled as contact against the symmetry plane.

•

You can find the option in the section of the node. Note that the computation of the Hencky strain can be further simplified by choosing the option, which corresponds to the following approximation:

Nonlinear_Structural_Materials_Module/Plasticity/compressed_elastoplastic_pipe

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > .

Click .

Click

In the tree, select > .

Click

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Value	Description
Ro	200[mm]	0.2 m	Outer radius
thic	25[mm]	0.025 m	Pipe wall thickness
Sy	250[MPa]	2.5E8 Pa	Yield stress
hfact	1	1	Hardening curve multiplier
maxLoad	6[MN]	6E6 N	The maximum load from the tool
par	0	0	Solution parameter
dist	0.1[mm]	1E-4 m	Distance of tool to pipe

Geometry 1

Circle: Pipe

In the toolbar, click

In the window for , type Circle: Pipe in the text field.

Locate the section. In the text field, type Ro.

In the text field, type 90.

Click to expand the section. In the table, enter the following settings:


Layer name	Thickness (m)
Layer 1	thic

Delete Entities 1 (del1)

In the window, right-click and choose .

In the window for , locate the section.

From the list, choose .

On the object , select Domain 1 only.

Click

Rectangle: Tool

In the toolbar, click

In the window for , type Rectangle: Tool in the text field.

Locate the section. In the text field, type 0.05.

In the text field, type Ro.

Locate the section. In the text field, type Ro+dist.

Fillet 1 (fil1)

In the toolbar, click

On the object , select Point 4 only.

In the window for , locate the section.

In the text field, type Ro/10.

Polygon: Symmetry Contact Plane

In the toolbar, click

In the window for , type Polygon: Symmetry Contact Plane in the text field.

Locate the section. From the list, choose .

Locate the section. In the table, enter the following settings:


x (m)	y (m)
0	0.15
0	-0.01

Click

Form Union (fin)

In the window, under > click .

In the window for , locate the section.

From the list, choose .

In the toolbar, click

Click the

button in the toolbar.

Definitions

General Contact Pair 1 (p1)

In the toolbar, click

and choose .

Solid Mechanics (solid)

Linear Elastic Material 1

In the window, under > click .

Plasticity 1

In the toolbar, click

and choose .

Select Domain 1 only.

In the window for , locate the section.

Find the subsection. From the list, choose .

Materials

Material 1 (mat1)

In the window, under right-click and choose .

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Young’s modulus	E	195[GPa]	Pa	Young’s modulus and Poisson’s ratio
Poisson’s ratio	nu	0.3	1	Young’s modulus and Poisson’s ratio
Density	rho	8000	kg/m³	Basic
Initial yield stress	sigmags	Sy	Pa	Elastoplastic material model

Add the hardening curve for the elastoplastic material.

In the window, expand the node, then click .

In the window for , locate the section.

Click

In the dialog, type plastic strain in the text field.

In the tree, select > .

Click .

Interpolation 1 (int1)

In the toolbar, click

and choose > .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file compressed_elastoplastic_pipe_stress_strain.txt.

In the text field, type hardFcn.

Locate the section. From the list, choose .

Locate the section. In the table, enter the following settings:


Argument	Unit
t	1

In the table, enter the following settings:


Function	Unit
hardFcn	MPa

Click

Material 1 (mat1)

In the window, under > click .

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Hardening function	sigmagh	hardFcn(epe)*hfact	Pa	Elastoplastic material model

Add a state variable to store the parameter value when the maximum load is reached.

Global Definitions

Click the

button in the toolbar.

In the dialog, select > in the tree.

In the tree, select the checkbox for the node > .

Click .

State Variables 1 (state1)

In the toolbar, click

and choose > .

In the window for , locate the section.

In the table, enter the following settings:


State	Initial value (1)	Update expression (1)	Description
peakParam	0	if(peakParam>0,peakParam,if(comp1.reacF>maxLoad,par,0))	Parameter value at max compression

From the list, choose .

Add variables computing the reaction force on the tool and the current pipe cross section.

Definitions

Integration 1 (intop1)

In the toolbar, click

and choose .

In the window for , type ReacInt in the text field.

Select Domain 2 only.

Locate the section. From the list, choose .

Integration 2 (intop2)

In the toolbar, click

and choose .

In the window for , type AreaInt in the text field.

Locate the section. From the list, choose .

Select Boundary 4 only.

Variables 1

In the window, right-click and choose .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Unit	Description
reacF	-2*ReacInt(solid.RFx)	N	Applied force from tool
area	AreaInt(-xsolid.nx)4	m²	Current cross section area for flow
disp	if(peakParam,peakParam-(par-peakParam)/2,par)[mm]	m	Compression function

Solid Mechanics (solid)

Symmetry 1

In the toolbar, click

and choose .

Select Boundaries 2 and 3 only.

Prescribed Displacement: Compression

In the toolbar, click

and choose .

In the window for , type Prescribed Displacement: Compression in the text field.

Select Domain 2 only.

Locate the section. From the list, choose .

From the list, choose .

In the u 0 x text field, type -disp.

Mesh 1

In the window, under click .

In the window for , locate the section.

From the list, choose .

Free Triangular 1

In the window, under > click .

In the window for , locate the section.

From the list, choose .

Select Domain 2 only.

Distribution 1

Right-click and choose .

In the window for , locate the section.

From the list, choose .

A minimal resolution is needed for the tool, since the whole domain is under displacement control. On the curved boundary, however, a fine resolution is needed.

Locate the section. In the text field, type 1.

Distribution 2

Right-click and choose .

Select Boundary 10 only.

In the window for , locate the section.

In the text field, type 30.

Mapped 1

In the toolbar, click

Distribution 1

Right-click and choose .

Select Boundaries 4 and 5 only.

In the window for , locate the section.

In the text field, type 80.

Distribution 2

In the window, right-click and choose .

Select Boundaries 2 and 3 only.

In the window for , locate the section.

From the list, choose .

In the text field, type 8.

In the text field, type 2.

Select the checkbox.

The symmetry boundary must also be meshed because of the contact, even though there is no physics attached to it.

Edge 1

In the toolbar, click

and choose .

Select Boundary 1 only.

Distribution 1

Right-click and choose .

In the window for , locate the section.

In the text field, type 1.

Click

Study 1

Step 1: Stationary

Set up an auxiliary continuation sweep for the parameter par. We do not know exactly how many parameter steps that are needed, but the last ones will be cheap to compute, since the tool is no longer in contact with the pipe.

In the window, under click .

In the window for , click to expand the section.

Select the checkbox.

Click