Double Cantilever Beam Model to Predict Mode I Fracture for Pinned and Unpinned Laminates

This post shows the finite element (FE) model development with Abaqus CAE [1] for a double cantilever beam (DCB) to predict the mode I fracture of the z-pinned and unpinned laminates [2, 3]. This study first develops a DCB test model in Abaqus CAE with the geometric information from (ASTM D5528) [3]. The report is organized into three fundamental sections according to the project requirement.

1. Development of the DCB test model

2. Development of the Finite element model 

3. Results and discussions of the DCB model test 

3.1. Unpinned laminate simulation

3.2. Z-pinned laminate simulation

3.3. Comparison between standard test and model

        3.4. Parametric study

1. Development of the DCB test model

For the geometrical shape and size of the DCB test setup, we choose the length of the beam as 100 mm where there are two parts: 70 mm length for the adhesive material and rest 30 mm length. The adhesive layer has a depth of 1 mm, and the beam part has also a depth of 1 mm. The following schematic drawing shows the geometry for the DCB model.

The following image depicts the implementation of the schematic geometry in the Abaqus CAE.

We have chosen materials for the cohesive and beam parts separately. The following table documents the major properties that have been implemented in the DCB test. The traction type is selected for the elastic material properties and its values are shown in the following table. For the damage criterion, the ‘Quads Damage’ option is chosen. In the ‘Damage Evolution’ module, we choose the type as ‘Energy’, then softening as ‘Linear’, then degradation as ‘Maximum’, then mixed mode behaviour as ‘BK’, then mixed mode ratio as ‘Energy’, and finally power set to 2.284. For the ‘Quads Damage’ option, we select the XFEM as ‘Normal’, where the tolerance is set to 0.05 and position is centroid.

2. Development of the Finite element model 

After the geometry is created, next we move on to develop the FE model for this study. At first, we begin with meshing the geometry. The beam is meshed as ‘Quad’ and ‘Structured’ from the Abaqus meshing option, where a CPE4R beam element is selected for plain strain problem. For the cohesive part of the adhesive materials, for the mesh control, we select ‘Quad’, then we further select the ‘Sweep’, and consider the portion of highlighted orientation. We chose the cohesive element with COH2D4, where we set the viscosity to 0.00005. The following image shows the meshing for the whole structure in Abaqus.

Next, we develop to create constraint in the DCB model so that the laminate is attached in between the cantilevers. We use tie constraint between the top surfaces of the cantilever beams to the bottom surfaces of the cohesive layers. In this way, we tie the cohesive material layer to both of the cantilever beams. The following image shows the tie constrain simulation in Abaqus.

3. Results and discussions of the DCB model test

3.1. Unpinned laminate simulation

We first show the results for the unpinned laminates in the DCB model. The following simulation results show the laminate separation for the unpinned configuration of the DCB test. 

The following result show the in-plane principal stress distribution in the separated DCB beam.

The following result shows the material orientations plots on the deformed shape of the beam.

3.2. Z-pinned laminate simulation

Next, we test the condition for the Z-pins, where pins are inserted in the laminate layers of the adhesive material. The pin material is chosen copper, where the Young’s modulus is 130 GPa and Poisson’s ratio is 0.34. The following image shows the position of a pin inside the DCB laminate. The z-pin is inside the cohesive layer of the DCB setup.

The following image shows the simulation result when a Z-pin is tied to the cohesive laminate.

3.3. Comparison between standard test and model

This section of the report discusses the comparison between the standard DCB test and the model test. We compare here force vs displacement plot for both the unpinned and pinned laminate configurations. From the following plots, we see that for both the unpinned (left) and pinned (right) laminates, there is a slight discrepancy between the standard and model tests. This is due to the facts the model test has slightly different geometry, different material properties, variations in the loads and boundary conditions, choice of the solvers, and so forth. However, the profiles are similar in both cases.

3.4. Parametric study

In this part, we do the parametric study for the Z-pin location and count within the laminate through the pin pull out test. We chose two pin materials, copper, and steel, in order to see the effects of them on the displacements of the adhesive layers with respect to the applied force between the standard test and the current model results. The following two tables represent the parametric study for the two materials. 

REFERENCES

[1] ABAQUS UNIFIED FEA (https://www.3ds.com/products-services/simulia/products/abaqus/)

[2] Blacklock, M, Joosten, MW, Pingkarawat, K, and Mouritz, AP, Prediction of mode I delamination resistance of z-pinned laminates using the embedded finite element technique, Composites: Part A 91 (2016) 283–291

[3] Standard Test Method for Mode I Interlaminar Fracture Toughness of Unidirectional Fiber-Reinforced Polymer Matrix Composites, Designation: D5528 – 13

#DoubleCantileverBeam   #FiniteElement  #FractureAnalysis  #ABAQUS  #LaminateSimulation

Stress and Strain Analysis of Linear Finite Elements

In this tutorial, we will use principle of virtual work method to determine the stress and displacement of the nodes of linear one dimensional finite elements. Following figure 1 represents a cantilever beam which is stressed axially. Let's assume the parameters L, A, E, c have the same value of 1 unit for simplicity. Our goal is to compare the finite element results with the exact solutions of a typical cantilever problem.

For this problem, we need to consider two quadratic one dimensional element. So, the displacement is interpolated among the three nodes and has the following form:

Fig.1: Showing one dimensional axially load bar given in the assignment problem.

Fig.2: Showing one dimensional piecewise linear interpolation.

Mathematica Program for Element 1

Clear[x, a0, a1, a2, u1, u2, u3, L]

(*Definition of quadratic form of displacement*)

u = a0 + a1*x + a2*x^2;

eq1 = (u /. x -> 0) - u1;

eq2 = (u /. x -> L/4) - u2;

eq3 = (u /. x -> L/2) - u3;

sol = Solve[{eq1 == 0, eq2 == 0, eq3 == 0}, {a0, a1, a2}];

u = u /. sol[[1]]

N1 = Coefficient[u, u1];

N2 = Coefficient[u, u2];

N3 = Coefficient[u, u3];

(*Shape Functions*)

Nn = {N1, N2, N3}

(*Applied Axial Load*)

p = c*x;

(*Nodal Force Vectors*)

fe = Integrate[Nn*p, {x, 0, L/2}];

c = 1;

L = 1;

fe // MatrixForm

Program Outputs

Displacement Vector for Element 1:

Shape Functions for Element 1:

Nodal Force Vector for Element 1:

Mathematica Program for Element 2

Clear[x, a0, a1, a2, u1, u2, u3, L]

(*Definition of quadratic form of displacement*)

u = a0 + a1*x + a2*x^2;

eq1 = (u /. x -> 0) - u1;

eq2 = (u /. x -> L/4) - u2;

eq3 = (u /. x -> L/2) - u3;

sol = Solve[{eq1 == 0, eq2 == 0, eq3 == 0}, {a0, a1, a2}];

u = u /. sol[[1]]

N1 = Coefficient[u, u1];

N2 = Coefficient[u, u2];

N3 = Coefficient[u, u3];

(*Shape Functions*)

Nn = {N1, N2, N3}

(*Applied Axial Load*)

p = c*(x + L/2);

(*Nodal Force Vectors*)

fe = Integrate[Nn*p, {x, 0, L/2}];

c = 1;

L = 1;

fe // MatrixForm

Program Outputs

Displacement Vector for Element 2

Shape Functions for Element 2

Nodal Force Vector for Element 2

Local Stiffness Matrix for Element 1

Clear[x, a0, a1, a2, u1, u2, u3, L, EA]

(*Definition of quadratic form of displacement*)

u = a0 + a1*x + a2*x^2;

eq1 = (u /. x -> 0) - u1;

eq2 = (u /. x -> L/4) - u2;

eq3 = (u /. x -> L/2) - u3;

sol = Solve[{eq1 == 0, eq2 == 0, eq3 == 0}, {a0, a1, a2}];

u = u /. sol[[1]];

N1 = Coefficient[u, u1];

N2 = Coefficient[u, u2];

N3 = Coefficient[u, u3];

B = {{D[N1, x], D[N2, x], D[N3, x]}};

(*Local Stiffness Matrix*)

K = Integrate[EA*Transpose[B].B, {x, 0, L/2}];

L = 1;

EA = 1;

K // MatrixForm

Program Outputs

Comparison between the Exact and Solution using Four Elements

We know that the exact solution for the axially loaded bar shown in Fig. 1 may be written as,

Where, for the consecutive four elements, we may write the displacement for each element as,

From the assignment problem, the parameter values are,

From the virtual work method, we may write,

The internal virtual work (IVW) may be expressed as,

And, the external virtual work may be presented as,

MATHEMATICA PROGRAM

The following program is written to determine the nodal displacement using the finite element method considering the principle of virtual work. 

Clear[x, a0, a1, a2, u1, u2, u3, L]

(*Definition of quadratic form of displacement*)

u = a0 + a1*x + a2*x^2;

eq1 = (u /. x -> 0) - u1;

eq2 = (u /. x -> L/4) - u2;

eq3 = (u /. x -> L/2) - u3;

sol = Solve[{eq1 == 0, eq2 == 0, eq3 == 0}, {a0, a1, a2}];

u = u /. sol[[1]];

(*Shape Functions*)

N1 = Coefficient[u, u1];

N2 = Coefficient[u, u2];

N3 = Coefficient[u, u3];

Nn = {N1, N2, N3};

B = {{D[N1, x], D[N2, x], D[N3, x]}};

c = 1;

L = 1;

EA = 1;

(*Axial Force on Element 1*)

p1 = c*x;

fe1 = Integrate[Nn*p1, {x, 0, L/2}];

(*Axial Force on Element 2*)

p2 = c*(x + L/2);

fe2 = Integrate[Nn*p2, {x, 0, L/2}];

K1 = Integrate[EA*Transpose[B].B, {x, 0, L/2}];

K2 = Integrate[EA*Transpose[B].B, {x, 0, L/2}];

Ndof = 5;

Ne = 2;

ID1 = {1, 2, 3};

ID2 = {3, 4, 5};

K = ConstantArray[0, {Ndof, Ndof}];

K[[ID1, ID1]] = K1 + K[[ID1, ID1]];

K[[ID2, ID2]] = K2 + K[[ID2, ID2]];

fe = ConstantArray[0, {Ndof, 1}];

fe[[ID1]] = fe1 + fe[[ID1]];

fe[[ID2]] = fe2 + fe[[ID2]];

IDE = {1};

IDN = {2, 3, 4, 5};

KR = Delete[K, IDE];

KR = Delete[Transpose[KR], IDE];

KR = Transpose[KR];

feR = Delete[fe, IDE];

(*Calculation of Nodal Displacements*)

DisR = LinearSolve[KR, feR];

Dis = ConstantArray[0, {Ndof, 1}];

Dis[[IDN]] = DisR + Dis[[IDN]];

Dis // MatrixForm

Program Outputs

So, we see that the exact and finite element solutions are the same.

From the following two plots, we notice that although the nodal displacement is continuous from the FEA result, it is not true for the stress. The stress distribution is indeed discontinuous.

MATHEMATICA PROGRAM FOR COMPARISON

(*Exact Solution*)

uexact = c*L^2/2/EA*x - c/6/EA*x^3;

(*Finite Element Approximation*)

N1 = Piecewise[{{4 x/L, 0 <= x <= L/4}, {4/L (L/2 - x), 

     L/4 <= x <= L/2}}, 0];

N2 = Piecewise[{{4/L (x - L/4), L/4 <= x <= L/2}, {4/L (3 L/4 - x), 

     L/2 <= x <= 3 L/4}}, 0];

N3 = Piecewise[{{4/L (x - L/2), L/2 <= x <= 3 L/4}, {4/L (L - x), 

     3 L/4 <= x <= L}}, 0];

N4 = Piecewise[{{4/L (x - 3 L/4), 3 L/4 <= x <= L}}, 0];

u1 = c*L^3/EA*(47/384);

u2 = c*L^3/EA*(11/48);

u3 = c*L^3/EA*(39/128);

u4 = c*L^3/EA*(1/3);

(*Calculating Nodal Displacement*)

u = u1*N1 + u2*N2 + u3*N3 + u4*N4;

uexactp = uexact /. {c -> 1, EA -> 1, L -> 1}

sigmaexact = D[uexactp, x]

up = u /. {c -> 1, EA -> 1, L -> 1}

sigmap = D[up, x]

(*Plotting to Show the Comparison between Exact and FEA Solution*)

Plot[{uexactp, up}, {x, 0, 1}, 

 AxesLabel -> {"x", "Nodal Displacement"}, 

 PlotLegends -> {"uexact", "u_Finite Element"}]

Plot[{sigmaexact, sigmap}, {x, 0, 1}, AxesLabel -> {"x", "sigma_11"}, 

 PlotLegends -> {"sigma_exact", "sigma_Finite Element"}]

For further information on finite element analysis and Mathematica, readers are encouraged to visit this page.

#FiniteElement #StressStrain #Mathematica #Abaqus #FEA #Blog #Blogger

Simulation of Buckling Behavior of a Pipe by Abaqus

In this tutorial, I am going to show how we can simulate the buckling characteristics of a pipe in Abaqus. The length of the pipe is 4.5 m, outside diameter (OD) is 0.9 m, and thickness is 9 mm. First, let us consider the stress-strain relationship, which is defined by the following table.

Table 1: Stress-Strain Relation for a Pipe

True Strain     True Stress (N/mm^2)

0                              0

0.001                     200

0.002                     400

0.003                     600

0.004                     800

0.0055                     900

0.008                     1000

0.0125                     1100

0.026                     1200

The relationship between the true and engineering stress-strain is governed by:

σ_true = F/A = (F/A₀) (1+ε_eng) = σ_eng(1+ε_eng)

Now, to find the elastic strain, the following relation is used.

ε^e = σ/E

The following table is based on the above definitions.

Table 2: Elastic and Plastic Strain for a Pipe

True Strain	True Stress (N/mm^2)	Tr Stress	Elastic Strain (εe)	Plastic Strain (εp)
0	0	0	0	0
0.001	200	200.2	0.001001	-1e-06
0.002	400	400.8	0.002004	-4e-06
0.003	600	601.8	0.003009	-9e-06
0.004	800	803.2	0.004016	-1.6e-05
0.0055	900	904.95	0.004525	0.000975
0.008	1000	1008	0.00504	0.00296
0.0125	1100	1113.75	0.005569	0.006931
0.026	1200	1231.2	0.006156	0.019844

The internal pressure inside the pipe is given by following relation:

P = [0.8(2t)σ_y]/OD

P = 12.8 MPa

This is the ultimate internal pressure of the pipe. Now, we are ready to create the model in Abaqus.

Steps in Abaqus Modelling:

Creating the Geometry

First, a part is created in Abaqus CAE, where a 3-point arc is chosen to draw a half circle with radius 450 mm. Since, the part is selected as 3D shell with extrusion, the length is 2250 mm. We will use symmetry and that’s why, these dimensions are selected to use symmetry boundary conditions about X and Z-axes.

The S4R shell element is chosen with reduced integration for analysis.

Meshing

Next, we create the mesh. The mesh size is set to 18 mm, which is twice the size of the shell thickness.

Creating Section

After that, the sections are created. The shell thickness is set 9 mm as per the problem.

Assigning Section

After creating the section, we assign it to the part. Then, we use “Render Shell Thickness” to show the shell thickness.

Creating Partition

A partition is created with 200 mm offset from one of the pipe edges.

Setting Material Properties

For the elastoplastic behaviour, the following material properties are given:

Creating Steps with Loads and Boundary Conditions

In this step, we define the boundary condition (BC) and load. A reference point (RP) is created to fix one end of the pipe. For the RP, a rigid body constraint is defined. Three BCs are generated as follow:

BC-1: for symmetry about X-axis

BC-2: for symmetry about Z-axis

BC-3: for the reference point (RP1), which is fixed along the axes.

Then, a pressure load of magnitude 12.8 MPa is applied inside the pipe.

Results and Discussions: 

Since the model is created with necessary boundary conditions and loads, in this section, we are interested to consider two circumstances.

• Buckling when both pressure and bending are considered

• Buckling when bending is considered without pressure

We apply the static Rik’s approach as a force control. The following figure shows the Von Mises stress when only pressure is applied inside the wall of the pipe.

Buckling when both pressure and bending are considered

At first, we apply a moment of 5x10^10 Nmm at the reference point by creating a new step where Rik’s static algorithm is selected. In this scenario, the pressure is still present inside the pipe. The following figure demonstrates the three steps for the simulation.

The following figure shows two relationship when pressure is present: moment vs arc length and rotation vs arc length. The intersection point between two graphs in the figure is the convergence point for the simulation.

Showing relationships between moment and arc length, and rotation and arc length with pressure

The following plot shows the relation between moment and angular displacement when pressure is present in the pipe.

Buckling when bending is considered without pressure

In this step, we just suppress the load leaving all other conditions same. We are interested to see the effect of the bending only on the buckling rate of the pipe.

Showing 2 steps: no pressure (left), and with angular displacement (right)



The following figure shows two relationship when pressure is absent: moment vs arc length and rotation vs arc length. The intersection point between two graphs in the figure is the convergence point for the simulation.

Showing relationships between moment and arc length, and rotation and arc length without pressure

Moment without pressure (left) and rotation without pressure (right) when arc length is 0.8016 mm

The following plot shows the relation between moment and angular displacement when pressure is absent in the pipe.

From the following comparison, we see that the buckling is higher when there is no pressure than the buckling with internal pressure.

To learn in detail of FEA, you may visit here.

#VonMisesstress #Buckling #Bending #FiniteElementAnalysis #FEA #Abaqus #ModellingSimulation #Blog #Blogger