▼Coreform IGA Verification Manual
►2 Solid Mechanics
On this page:
2.1 Linear Stress Analysis
2.1.1 Cantilever beam with end shear load
2.1.1.1 Objective
2.1.1.2 Geometry Description
2.1.1.3 External Interactions
2.1.1.4 Material Description
2.1.1.5 Method
2.1.1.6 Solution Characteristics
2.1.1.7 Results
2.1.2 Composite cantilever beam with end shear load
2.1.2.1 Objective
2.1.2.2 Geometry Description
2.1.2.3 External Interactions
2.1.2.4 Material Description
2.1.2.5 Method
2.1.2.6 Solution Characteristics
2.1.2.7 Results
2.1.3 Infinite plate with a hole
2.1.3.1 Objective
2.1.3.2 Geometry Description
2.1.3.3 External Interactions
2.1.3.4 Material Description
2.1.3.5 Method
2.1.3.6 Results
2.1.4 Stresses in long pressurized pipe with butt-welded supporting trunnion pipes
2.1.4.1 Objective
2.1.4.2 Geometry Description
2.1.4.3 External Interactions
2.1.4.4 Material Description
2.1.4.5 Method
2.1.4.6 Solution Characteristics
2.1.4.7 Results
2.1.5 Stress Concentrations:   filleted square shoulder of rectangular section
2.1.5.1 Objective
2.1.5.2 Geometry Description
2.1.5.3 External Interactions
2.1.5.4 Material Description
2.1.5.5 Method
2.1.5.6 Solution Characteristics
2.1.5.7 Results
2.2 Large Deformation
2.2.1 Large and small deformation of composite geometry with material interface connections and merged surfaces
2.2.1.1 Objective
2.2.1.2 Geometry Description
2.2.1.3 External Interactions
2.2.1.4 Material Description
2.2.1.5 Method
2.2.1.6 Solution Characteristics
2.2.1.7 Results
2.2.2 Mooney-Rivlin material verification
2.2.2.1 Objective
2.2.2.2 Geometry Description
2.2.2.3 External Interactions
2.2.2.4 Material Description
2.2.2.5 Method
2.2.2.6 Solution Characteristics
2.2.2.7 Results
2.3 Incompressible Elasticity
2.3.1 Neo-Hookean cylindrical pressure vessel
2.3.1.1 Objective
2.3.1.2 Geometry Description
2.3.1.3 External Interactions
2.3.1.4 Material Description
2.3.1.5 Method
2.3.1.6 Solution Characteristics
2.3.1.7 Results
2.4 Plasticity
2.4.1 J2 plasticity:   3D unit cube with isotropic hardening
2.4.1.1 Objective
2.4.1.2 Geometry Description
2.4.1.3 External Interactions
2.4.1.4 Material Description
2.4.1.5 Method
2.4.1.6 Results
2.5 Mechanical Contact
2.5.1 Interference fit of two pipes
2.5.1.1 Objective
2.5.1.2 Geometry Description
2.5.1.3 External Interactions
2.5.1.4 Material Description
2.5.1.5 Method
2.5.1.6 Solution Characteristics
2.5.1.7 Results
    contents  ← prev  up  next → 

2 Solid Mechanics

    2.1 Linear Stress Analysis

    2.2 Large Deformation

    2.3 Incompressible Elasticity

    2.4 Plasticity

    2.5 Mechanical Contact

2.1 Linear Stress Analysis

    2.1.1 Cantilever beam with end shear load

    2.1.2 Composite cantilever beam with end shear load

    2.1.3 Infinite plate with a hole

    2.1.4 Stresses in long pressurized pipe with butt-welded supporting trunnion pipes

    2.1.5 Stress Concentrations: filleted square shoulder of rectangular section

2.1.1 Cantilever beam with end shear load

    2.1.1.1 Objective

    2.1.1.2 Geometry Description

    2.1.1.3 External Interactions

    2.1.1.4 Material Description

    2.1.1.5 Method

    2.1.1.6 Solution Characteristics

    2.1.1.7 Results

2.1.1.1 Objective

This verification problem demonstrates that higher-order elements can capture the exact solution of a simple cantilever beam problem with coarse meshes. In this problem we analyse the normal and shear stress profiles of a cantilever beam with a square cross-section that is subjected to a tip load using quadratic and cubic basis functions. The results show that both choices of basis functions exactly capture the axial stress and that cubic elements capture both the axial and shear stresses exactly on any mesh.

Although real world engineering problems will not in general have solutions that can be captured exactly with quadratic or cubic basis functions, these higher-order functions will typically be more accurate for a given number of degrees of freedom. This allows coarser meshes to be used in simulations, which can lead to reduced computational cost in many problems.

2.1.1.2 Geometry Description

The beam geometry is an axis-aligned box with dimensions , , and in the , , and directions, respectively. The direction is considered to be the axial direction, and the tip load is applied in the direction. The overall setup is depicted in figure 1. The specific dimensions are given in table 1.

Figure 1: A sketch of the cantilever beam.

   

   

   

Table 1: Geometric parameters.

2.1.1.3 External Interactions

A transverse shear traction is applied at the free end of the beam with an average shear stress of in the direction. As the cross-sectional area is , the total applied force is . The shear traction applied to the tip is distributed according to the parabolic solution from Euler–Bernoulli beam theory given in equation , to facilitate precise quantitative accuracy testing. The cantilever support at the fixed end is modeled by fixing displacement in both transverse directions ( and ), and applying a linearly-varying pressure distribution, coming from the Euler–Bernoulli axial stress solution given in equation . The remaining rigid modes are eliminated by constraining axial displacement at only the corners of the cantilever support cross-section. Although the Euler–Bernoulli solution is formally derived through an assumption of cross-sections remaining planar (alongside other approximations), that does not hold exactly for the full three-dimensional displacement solution, and the support cross-section needs to remain free to warp in the axial direction. The remaining rigid modes are eliminated by constraining axial displacement at only the corners of the cantilever support cross-section.

Due to St. Venant’s principle, qualitatively-similar results could be obtained using a fixed boundary condition for the cantilever support and a uniformly-distributed shear traction for the tip load. However, the cubic solution would no longer be exact and significant deviations from the beam-theory solution would occur near the ends of the beam.

2.1.1.4 Material Description

The beam is modeled as an isotropic linear-elastic material. The material properties used are shown in table 2. Note that the use of a zero Poisson ratio ensures that all deformation occurs in the plane.

Young's Modulus

   

1.0

Poisson's Ratio

   

0.0

Table 2: Material properties.

2.1.1.5 Method

We approximate the beam’s displacement field using quadratic and cubic trimmed splines defined on coarse and refined background meshes. Both meshes are axis-aligned. The coarse mesh has element sizes of in the transverse directions and in the axial direction. These element sizes are halved in the refined mesh. Although the element sizes used could evenly-divide the beam geometry, the background meshes are aligned to produce a nontrivial trimmed spline, as shown in figure 2.

Figure 2: Coarse trimmed spline discretization.

2.1.1.6 Solution Characteristics

Since the cross-section of the beam is a square and is composed of a single homogenous material the neutral axis is located at the center of the beam. Thus the normal stress profile can be computed as where is the moment at the axial location being probed and is the area moment of inertia of the cross-section, which is given by for a rectangle. The shear stress profile can be calculated as In the present problem, the shear stress calculation can be simplified to which we note is a parabolic distribution with bounds

2.1.1.7 Results

The normal stress is plotted along a transverse line at a point midway along the length of the beam in figure 3. The reference solution given in equation is also plotted. All plots are coincident because both quadratic and cubic basis functions are able to represent the solution exactly.

Figure 3: Normal stress profile for immersed cantilever beam.

The computed and reference shear stresses are plotted along the same line in figure 4. Cubic basis functions are able to represent the solution exactly, so the resulting plot is coincident with the reference solution. The shear stress computed using a quadratic basis is converging at the optimal rate to the reference solution.

Figure 4: Shear stress profile for immersed cantilever beam.

2.1.2 Composite cantilever beam with end shear load

    2.1.2.1 Objective

    2.1.2.2 Geometry Description

    2.1.2.3 External Interactions

    2.1.2.4 Material Description

    2.1.2.5 Method

    2.1.2.6 Solution Characteristics

    2.1.2.7 Results

2.1.2.1 Objective

This verification problem demonstrates the accuracy of tied constraints between adjacent parts by computing the stress profile through a composite cantiliver beam composed of a stiff and flexible section. The normal and shear stress are compared to an exact reference solution.

Tied constraints play a critical role when computing results for assemblies using the flex method. In the typical body-fitted approach for assemblies, adjacent parts are often tied together by creating a conformal mesh between the two parts and specifying that coincident nodes at the interface map to the same degrees of freedom. For the flex approach, adjacent parts will each have a different background spline, which in general will not align or share coincident points at the interface. Because of this, tied constraints are used to tie adjacent parts together. In this simple problem, two adjacent cantilever beams are tied together and then subjected to an end load. Results are computed for two cases. In the first case, the background mesh for each section has the same element size. For the second case, the element sizes are different. The results for both cases show that these interfaces can be accurately modeled with the flex method’s tied constraint approach.

2.1.2.2 Geometry Description

The overall cross section of the composite beam is a square, which is then split lengthwise into an upper and lower section, as shown in figure 5. The dimensions of the stiff and flexible sections are indicated by and , respectively. The dimensions are given in table 3.

Figure 5: A sketch of the composite cantilever beam showing the stiff and flexible sections that comprise the beam.

Dimension

   

Value

   

   

   

   

   

Table 3: Geometric parameters.

2.1.2.3 External Interactions

Displacements in the transverse ( and ) directions are contrained on the fixed end of the cantilever beam. The axial displacement is only constrained at the corners of the cantilever support face of the beam’s lower section, to remove rigid modes, while the analytical distribution of axial stress from Euler–Bernoulli beam theory given in equation is applied on the fixed end instead of constraining axial displacement. This is because the analytical solution depends on allowing cross-sections to warp slightly, and we want to be able to perform precise testing of the method’s accuracy. A transverse shear traction is applied at the free end of the beam with a non-dimensional average shear stress of , as shown in figure 5. The beam has a cross-sectional area of , which results in a total applied force of . The traction is distributed according to the analtyical beam-theory solution given in equation , to facilitate accuracy testing. The interface between the two sections is modeled as a material interface, using a dimensionless penalty scale factor of . This is selected higher than the default value recommended for practical analyses, to enable precise comparisons with the analytical solution.

2.1.2.4 Material Description

Both sections of the beam are modeled as isotropic linear-elastic materials. The non-dimensional material properties for the stiff and flexible sections of the beam are given in table 4.

Property

   

Flexible

   

Stiff

Young's Modulus

   

   

Poisson's Ratio

   

   

Table 4: The non-dimensional material properties for the stiff and flexible layers of the cantilever beam.

2.1.2.5 Method

Results are shown in figure 6 for the composite beam, computed with quadratic and cubic U-splines for the two trimmed U-spline cases. Figure 6a shows a representative background spline for one of the beam sections. For the trimmed U-spline shown in figure 6b the background meshes were aligned so that surfaces elements from one section were coincident with surfaces elements on the the other section. We also computed results for the unaligned case shown in figure 6c to show robustness of tied constraints between parts with trimmed elements of differents sizes.

Figure 6a: The untrimmed background spline and immersed beam for one of the sections.

Figure 6b: The trimmed beam sections for the case where the background meshes for each section are aligned.

Figure 6c: The trimmed beam sections for the case where the background meshes are not aligned.

Figure 6: The trimming process for the composite beam.

2.1.2.6 Solution Characteristics

A reference solution for the composite beam can be calculated using the transformed section method, wherein the geometry of the beam is transformed into an imaginary section of uniform material with equivalent elastic properties. We assume here that the notional uniform beam material is equivalent to the material found in the composite beam’s flexible section, so we transform only the geometry of the stiffer material.

The normal stress profile can then be calculated as

where is the distance to the neutral axis of the section and is radius of curvature. Intermediate calculations can be computed via the following expressions:

The shear stress profile can then be calculated as

where intermediate calculations can be computed via the following expressions:

2.1.2.7 Results

The normal stress is plotted along a transverse line at a point midway along the length of the beam in figure 7. The reference solution given in equation is also plotted. All plots are coincident because both quadratic and cubic basis functions are able to represent the solution exactly.

Figure 7a: Immersed with aligned trimmed U-spline.

Figure 7b: Immersed with non-matching trimmed U-spline.

Figure 7: Normal stress profiles. The computed normal stress is coincident with the reference solution in all cases.

The computed and reference shear stress plotted along the same line is shown in figure 8. As noted, cubic basis functions are able to represent the solution exactly and the resulting plot is coincident with the reference solution. The shear stress computed using a quadratic basis is converging at the optimal rate to the reference solution.

Figure 8a: Immersed with aligned trimmed U-spline.

Figure 8b: Immersed with non-matching trimmed U-spline.

Figure 8: Shear stress profiles. The shear stress computed with the cubic basis is coincident with the reference solution in all cases.

In all cases, the interface between the stiff and flexible sections of the beam was modeled as a tied constraint. These results demonstrate that we are able to accurately represent material interfaces with tied constraints.

2.1.3 Infinite plate with a hole

    2.1.3.1 Objective

    2.1.3.2 Geometry Description

    2.1.3.3 External Interactions

    2.1.3.4 Material Description

    2.1.3.5 Method

    2.1.3.6 Results

2.1.3.1 Objective

This verification problem shows that the Flex method is able to achieve optimal theoretical mesh convergence rates on the well-known infinite plate with a hole benchmark problem.

This problem is chosen as a verification problem not only because the exact solution is known, but in particular because the geometry includes a hole. The presence of a hole in the geometry presents two problem features that are critical for verification of flex-method stress analysis. First, the hole will produce a stress concentration. Second, the hole will require elements of the background spline to be trimmed in the area of the stress concentration. Accurate representation of stresses near features that produce stress concentrations is critical for stress analysis, and such features will almost always produce trimmed elements. The results presented for this problem show that we are able to achieve theoretical mesh convergence with trimmed elements. We also show that stress gradients can be captured to a reasonable level of accuracy on relatively coarse meshes. These results show that the flex method is a poweful tool for stress analyis of complex parts.

2.1.3.2 Geometry Description

The in-plane geometry for this problem is a quarter-symmetry representation of a square plate with a hole in the center, shown in figure 9. For this problem, and . The out-of-plane thickness is selected to be the mesh element size in each computation, so that there is always only one element in the thickness direction.

Figure 9: A quarter-symmetry representation of a infinite plate with a hole.

2.1.3.3 External Interactions

The model is constrained using symmetry boundary conditions on the left and bottom surfaces, and plane strain deformations are enforced by constraining out-of-plane displacements on both the front and back surfaces. The exact stress is given by in a polar coordinate system with the origin at the center of the hole, where is the traction applied in the direction at infinity. This stress is used to apply the exact tractions on the top and right surfaces, as shown in figure 9.

2.1.3.4 Material Description

The non-dimensional material properties are given in table 5.

Dimension

   

Value

Young's Modulus

   

Poisson's Ratio

   

0.33

Table 5: Non-dimensional material properties.

2.1.3.5 Method

Results for this problem were computed in Coreform IGA for the seven levels of refinement listed in table 6. A quadratic basis with continuity was used in each case.

Mesh

   

Mesh Size

   

Edge Count

   

   

   

   

   

   

   

   

   

   

   

   

   

   

Table 6: Mesh sizes and edge counts for the convergence study.

Figure 10: The background mesh for mesh #3. There are approximately eight elements along the long edges of the part.

The trimmed U-splines for Mesh 1, Mesh 3, and Mesh 7 are shown in figure 11.

Figure 11a: Mesh 1.

Figure 11b: Mesh 3.

Figure 11c: Mesh 7.

Figure 11: The trimmed U-splines for Mesh 1, Mesh 3, and Mesh 7.

2.1.3.6 Results

The error and contour plots of the stress component are plotted in figure 12 and figure 13, respectively. Figure 12 shows that the flex method achieves optimal convergnce rates for this problem. The errors plotted in this figure are in a two-dimensional error norm computed by integrating the three-dimensional error, then normalizing by the square root of element size to account for the mesh-dependent depth of the domain. We also see in figure 13a that a very coarse mesh is able to qualitatively capture the stress for this problem quite well.

Figure 12: Convergence plot of the error in the stress component.

Figure 13a: Mesh size 1.

Figure 13b: Mesh size 0.25.

Figure 13c: Mesh size .

Figure 13: Countour plots of the stress component.

2.1.4 Stresses in long pressurized pipe with butt-welded supporting trunnion pipes

    2.1.4.1 Objective

    2.1.4.2 Geometry Description

    2.1.4.3 External Interactions

    2.1.4.4 Material Description

    2.1.4.5 Method

    2.1.4.6 Solution Characteristics

    2.1.4.7 Results

2.1.4.1 Objective

In this verification problem, the maximum von Mises stress is computed for a supported, pressurized pipe under a tension load. The pipe is supported by two trunnions that have been butt-welded to the pipe as shown in figure 14. The results are compared with the reference solution published in the NAFEMS report  [3].

The Infinite plate with a hole verification problem showed that the flex method is able to accurately represent stresses in trimmed elements for a simple geometry where an exact stress field is known. In this problem, the geometry is more representative of a real-world part wherein an exact stress is not known because of the geometric complexity. In general, as the complexity of CAD parts increases, the complexity of the trimmed elements will also increase. However, it is important to keep in mind that the trimmed elements are only used to define the region in which evaluations of the background spline will be performed. The computations themselves will always be performed on the well-behaved functions of the background spline, which provides accurate results in trimmed elements. The presented results show that the flex method is able to accurately predict the maximum von Mises stress for this problem using a relatively coarse background spline.

A comparison between a conformal body-fitted approach and the fully-immersed flex approach is also shown in the results. The body-fitted approach required additional work to decompose the geometry into hex-meshable sections. The results show that both approaches produce excellent stress approximations for relatively coarse meshes. The flex approach, however, requires much less work because there is no need to decompose the geometry.

2.1.4.2 Geometry Description

The geometry for this problem is shown in figure 14 and the dimensions are specified in table 7.

Figure 14: The geometry of the supported pressurized pipe (not to scale).

Dimension

   

Value

   

   

   

   

   

   

   

Table 7: Dimensions of the supported pressurized pipe geometry.

2.1.4.3 External Interactions

The pipe is loaded by an internal pressure and is also under a uniformly distributed tensile load at the ends of the pipe. These loads are specified in table 8. No loads are applied to the trunnions. The problem is modeled under symmetry assumptions in the three coordinate planes and no other boundary conditions or loads are applied.

   

   

Table 8: Loading conditions for the pressurized pipe.

2.1.4.4 Material Description

The material properties are given in table 9.

Property

   

Value

Mass Density

   

Young's Modulus

   

Poisson's Ratio

   

Table 9: Material properties.

2.1.4.5 Method

Results for this problem were computed in Coreform IGA for both a body-fitted and immersed mesh, each at three different mesh sizes. The body-fitted 5mm mesh is shown in figure 15a. The corresponding untrimmed 5mm background spline and trimmed U-spline for the immersed case are shown in figure 15b and figure 15c.

Figure 15a: The 5mm body-fitted mesh.

Figure 15b: The untrimmed background spline with the immersed quarter-symmetry part.

Figure 15c: The trimmed U-spline.

Figure 15: The 5mm body-fitted mesh at top, followed by the untrimmed background spline and trimmed U-spline for the 5mm immersed case.

2.1.4.6 Solution Characteristics

The reference solution for this problem was published by NAFEMS in  [3]. Convergence results for the maximum von Mises stress from this report are shown in table 10.

Mesh Convergence Iteration

   

Mesh Size (mm)

   

Maximum Von Mises Stress (MPa)

   

   

   

   

   

   

   

   

Table 10: Reference solution convergence results. The mesh size column reports the mesh edge length in the area where the maximum stress occurs.

2.1.4.7 Results

Linear elastic results were computed in Coreform IGA for both body-fitted and immersed cases. The maximum von Mises stress at the integration points was then found and is summarized in table 11. For both body-fitted and immersed cases, the results compare very well to the reference solution, with less than 5 percent error on the coarsest mesh. Displacement and von Mises stress contours are shown in figure 16.

Mesh scheme

   

Mesh size (mm)

   

Max. von Mises stress (MPa)

   

Percent difference

Body fit coarse

   

10

   

508

   

4.8

Body fit medium

   

5

   

516

   

3.4

Body fit fine

   

2.5

   

526

   

1.5

Immersed coarse

   

10

   

522

   

2.2

Immersed medium

   

5

   

533

   

0.1

Immersed fine

   

2.5

   

534

   

0.1

Table 11: Convergence results for the body-fitted and immersed simulations.

Figure 16a: Immersed von Mises stress.

Figure 16b: body-fitted von Mises stress.

Figure 16c: Immersed displacement.

Figure 16d: body-fitted displacement.

Figure 16: Comparison of immersed versus body-fitted von Mises stress and displacement for the fine mesh. The von Mises stress is shown in units of megapascals (MPa) and displacements are given in units of millimeters (mm). The minimum value for von Mises stress contours is set at 400 MPa to better show contours around the maximum value.

2.1.5 Stress Concentrations: filleted square shoulder of rectangular section

    2.1.5.1 Objective

    2.1.5.2 Geometry Description

    2.1.5.3 External Interactions

    2.1.5.4 Material Description

    2.1.5.5 Method

    2.1.5.6 Solution Characteristics

    2.1.5.7 Results

2.1.5.1 Objective

Predict the stress concentration factor in a member of rectangular section with a filleted square shoulder in axial tension.

2.1.5.2 Geometry Description

Figure 17: Roark 17.1 5A Geometry (not to scale)

Dimensions of the problem are listed below:

   

   

   

   

   

Table 12: Problem Geometry

2.1.5.3 External Interactions

One end is fixed, and the other is given an applied pressure of .

2.1.5.4 Material Description

Mass Density

   

Young's Modulus

   

Poisson's Ratio

   

Table 13: Material Properties

2.1.5.5 Method

The displacements are approximated for two cases; an immersed approach and a body fit approach. In the body fit approach the volume is meshed in Coreform Cubit with a requested mesh size of 2.0. In the immersed approach the volume is immersed in a rectilinear background mesh with a mesh size of 2.0. Each case uses a quadratic spline basis to approximate displacements.

2.1.5.6 Solution Characteristics

The stress concentration factor, , is calculated using a cubic polynomial approximation:

For the cases where:

Coefficients are tabulated for the following additional conditions:

   

   

   

   

   

   

   

   

   

   

Table 14: Approximation coefficients

For the specific problem at hand we calculate the stress concentration factor to be . The expected maximum principal stress is .

2.1.5.7 Results

Mesh Scheme

   

Maximum Principal Stress

   

Percent Difference from Analytic Solution

BodyFit

   

   

4.9 %

Immersed Rectilinear

   

   

7.8 %

Table 15: Approach Comparison

Figure 18a: Immersed rectilinear Stress

Figure 18b: BodyFit Stress

Figure 18c: Immersed rectilinear Displacement

Figure 18d: BodyFit Displacement

Figure 18: Approach Comparison

A line probe from the tip of the stress concentration out 10 mm in the x direction was used for comparison between the two methods:

Figure 19: Method Comparison - Line Probe from Tip out 10 mm in x

A line probe from the tips of each stress concentration (from tip to tip) was also used for comparison between the two methods:

Figure 20: Method Comparison - Line Probe from Tip to Tip

The table below summarizes the solution characteristics between an immersed rectilinear mesh approach and a bodyfit (traditional) approach:

Problem Time Metric

   

BodyFit

   

Immersed rectilinear

Initialization Time

   

1.097e-05

   

1.216e-05

Linear Equation Assembly Time

   

0.145

   

0.167

NonLinear Equation Assembly Time

   

0.563

   

0.458

Solver Time

   

0.833

   

0.764

Table 16: Problem Time Comparison

Problem Size Metric

   

BodyFit

   

Immersed rectilinear

Quadrature Points

   

80433

   

88263

Degrees of Freedom

   

53901

   

33525

Table 17: Problem Size Comparison

2.2 Large Deformation

    2.2.1 Large and small deformation of composite geometry with material interface connections and merged surfaces

    2.2.2 Mooney-Rivlin material verification

2.2.1 Large and small deformation of composite geometry with material interface connections and merged surfaces

    2.2.1.1 Objective

    2.2.1.2 Geometry Description

    2.2.1.3 External Interactions

    2.2.1.4 Material Description

    2.2.1.5 Method

    2.2.1.6 Solution Characteristics

    2.2.1.7 Results

2.2.1.1 Objective

The objective of this problem is to demonstrate that Coreform IGA can accurately predict small and large displacements. Additionally, it demonstrates displacement accuracy in a composite body in which the components are connected via merged geometry in one case and via material interfaces in another. Finally this problem demonstrates displacement accuracy with both linear and nonlinear material models. The problem is a simple patch test in which the displacements at the boundaries are prescribed and the displacements on the interior are easily validated against an analytical solution.

2.2.1.2 Geometry Description

The geometry for this problem consists of a single 1x1x1 cube that is made up of seven component cuboid shapes as shown in figure 21.

Figure 21a: Assembled geometry

Figure 21b: Exploded geometry

Figure 21: Geometry of the finite deformation patch test

The interfaces between the components are purposely skewed so to be non planar and non square. The corner vertex coordinates of the center cuboid are shown in table 18.

   

   

   

   

   

   

   

   

   

   

   

   

   

   

   

   

Table 18: Center cuboid corner vertex coordinates

2.2.1.3 External Interactions

The displacement on each boundary face is prescribed according the following function is the current time, and are the spatial coordinates. The output of this function is scaled by a parameter denoted by that varies between test cases to produce large and small displacements.

2.2.1.4 Material Description

The analysis is run with both a linear elastic material and a nonlinear neohookean material model. In each case the material parameters used are the same and are defined in table 19.

Young's Modulus

   

Poisson's ratio

   

Mass Density

   

Table 19: Material properties.

2.2.1.5 Method

The test is run with several different configurations enumerated in table 20

Case Number

   

Mesh Type

   

Material Model

   

Interface connection

   

Displacement Scale

   

Body-fit

   

Linear elastic

   

Merged

   

   

Body-fit

   

Neohookean

   

Merged

   

   

Body-fit

   

Neohookean

   

Merged

   

   

Body-fit

   

Neohookean

   

Merged

   

   

Body-fit

   

Neohookean

   

Material Interface

   

   

Body-fit

   

Neohookean

   

Material Interface

   

   

Body-fit

   

Neohookean

   

Material Interface

   

   

Immersed

   

Neohookean

   

Material Interface

   

Table 20: Test configurations

In each case, the displacements are approximated with a quadratic spline basis.

For the body-fit cases, the composite body is meshed in Coreform Cubit. Each component volume is assigned a requested interval size of 1, so in the resulting mesh each component cuboid is a single hex element.

For the immersed case, the composite cube is immersed in a rectilinear background grid with a mesh size of 0.5 in each coordinate direction. Note that in the immersed case we use quadrature rule of degree 5, which is higher than the default of QP1.

For the merged geometry cases, a merge operation is performed in Coreform Cubit that merges the interface surfaces. This forces the resulting body fit mesh to be contiguous and watertight.

For the material interface cases, interface surfaces are automatically detected in Coreform Flex and a tied constraint is enforced between each pair of interface surfaces.

The displacement scale is a parameter that scales the output of the displacement boundary condition function described in the interactions section. A value of corresponds to a large displacement scenario, corresponds to a small displacement and is somewhere in between.

A probe is defined to measure the displacements at the center of the cube. Another is defined to measure displacements at the location , which is the location of maximum displacement.

2.2.1.6 Solution Characteristics

The expected displacement at any point and time can be calculated with the equation used for the prescribed boundary displacements.

2.2.1.7 Results

A visualization of the displacements for each displacement scale is shown in figure 22 - figure 24. The probe results and error values for each case are tabulated in table 21 - table 28.

Figure 22: Displacement scale

Figure 23: Displacement scale

Figure 24: Displacement scale

Coordinate Direction

   

Center Expected Value

   

Center Measured Value

   

Center Relative Error

   

Max Expected Value

   

Max Measured Value

   

Max Relative Error

x

   

   

   

0.00 %

   

   

   

0.00 %

y

   

   

   

0.00 %

   

   

   

0.00 %

z

   

   

   

0.00 %

   

   

   

0.00 %

Table 21: Probe results and error values for case #1

Coordinate Direction

   

Center Expected Value

   

Center Measured Value

   

Center Relative Error

   

Max Expected Value

   

Max Measured Value

   

Max Relative Error

x

   

   

   

0.00 %

   

   

   

0.00 %

y

   

   

   

0.00 %

   

   

   

0.00 %

z

   

   

   

0.00 %

   

   

   

0.00 %

Table 22: Probe results and error values for case #2

Coordinate Direction

   

Center Expected Value

   

Center Measured Value

   

Center Relative Error

   

Max Expected Value

   

Max Measured Value

   

Max Relative Error

x

   

   

   

0.00 %

   

   

   

0.00 %

y

   

   

   

0.00 %

   

   

   

0.00 %

z

   

   

   

0.00 %

   

   

   

0.00 %

Table 23: Probe results and error values for case #3

Coordinate Direction

   

Center Expected Value

   

Center Measured Value

   

Center Relative Error

   

Max Expected Value

   

Max Measured Value

   

Max Relative Error

x

   

   

   

0.00 %

   

   

   

0.00 %

y

   

   

   

0.00 %

   

   

   

0.00 %

z

   

   

   

0.00 %

   

   

   

0.00 %

Table 24: Probe results and error values for case #4

Coordinate Direction

   

Center Expected Value

   

Center Measured Value

   

Center Relative Error

   

Max Expected Value

   

Max Measured Value

   

Max Relative Error

x

   

   

   

0.00 %

   

   

   

0.00 %

y

   

   

   

0.00 %

   

   

   

0.00 %

z

   

   

   

0.00 %

   

   

   

0.00 %

Table 25: Probe results and error values for case #5

Coordinate Direction

   

Center Expected Value

   

Center Measured Value

   

Center Relative Error

   

Max Expected Value

   

Max Measured Value

   

Max Relative Error

x

   

   

   

0.00 %

   

   

   

0.00 %

y

   

   

   

0.00 %

   

   

   

0.00 %

z

   

   

   

0.00 %

   

   

   

0.00 %

Table 26: Probe results and error values for case #6

Coordinate Direction

   

Center Expected Value

   

Center Measured Value

   

Center Relative Error

   

Max Expected Value

   

Max Measured Value

   

Max Relative Error

x

   

   

   

0.00 %

   

   

   

0.00 %

y

   

   

   

0.00 %

   

   

   

0.00 %

z

   

   

   

0.00 %

   

   

   

0.00 %

Table 27: Probe results and error values for case #7

Coordinate Direction

   

Center Expected Value

   

Center Measured Value

   

Center Relative Error

   

Max Expected Value

   

Max Measured Value

   

Max Relative Error

x

   

   

   

0.01 %

   

   

   

0.00 %

y

   

   

   

0.01 %

   

   

   

0.00 %

z

   

   

   

0.00 %

   

   

   

0.00 %

Table 28: Probe results and error values for case #8

2.2.2 Mooney-Rivlin material verification

    2.2.2.1 Objective

    2.2.2.2 Geometry Description

    2.2.2.3 External Interactions

    2.2.2.4 Material Description

    2.2.2.5 Method

    2.2.2.6 Solution Characteristics

    2.2.2.7 Results

2.2.2.1 Objective

This problem demonstrates that the results obtained when using the Mooney-Rivlin material model are consistent with analytical values derived from the underlying material formulation. The problem is setup to follow the biaxial tension test described in Cheng and Zhang [1] and uses the equations derived there to calcuated an analytical expected value for the displacment of the material.

2.2.2.2 Geometry Description

The geometry for this problem consists of a single axis-aligned cube with length 1.0 on each side as shown in figure 25

Figure 25: Patch test geometry and loading

2.2.2.3 External Interactions

The boundary and load conditions imposed on the model are described in table 29 and table 30.

Face

   

Constraint variable

   

Constraint value

-X

   

X Displacement

   

-Y

   

Y Displacement

   

-Z

   

Z Displacement

   

Table 29: Model boundary conditions

Face

   

Load type

   

Value

+X

   

Surface Pressure

   

+Y

   

Surface Pressure

   

Table 30: Model load conditions

2.2.2.4 Material Description

The analysis is run with a Mooney-Rivlin material model, with parameters shown in table 31.

Shear modulus 1

   

   

5.95522e+05

Shear modulus 2

   

   

5.03810e+04

Bulk modulus

   

   

1.0e+11

Table 31: Material properties.

2.2.2.5 Method

The cube is meshed with a single hex element in Coreform Cubit. The boundary conditions and load conditions are defined in Coreform Flex. A probe is defined at the corner formed by the intersection of the +X, +Y and +Z faces to record the displacements in each coordinate direction. The problem is solved in Coreform IGA with the continuation timestepper, resulting in a static solution.

2.2.2.6 Solution Characteristics

As outlined in Cheng and Zhang [1], the setup of this test leads to a simplfied expression for the nominal stress in the body

where is the nominal stress, is the stretch, defined as the ratio between the deformed length and the original length, and , are the 1st and 2nd shear modulii for the material. This expression is derived from an expression of the 2nd Piola Kirchoff stress tensor as a function of the partial derivatives of the Mooney-Rivlin strain energy density. The nominal stress at our probe location is a known value uniquely determined by the applied pressure, so the above expression can be used to solve for the expected stretch and by extension the expected displacements at the probe location.

2.2.2.7 Results

A visualization of the displacement field is shown in figure 26. The probe results and error values are tabulated in table 32.

Coordinate Direction

   

Expected displacement

   

Measured displacement

   

Relative error

x

   

   

   

0.02 %

y

   

   

   

0.02 %

z

   

   

   

0.00 %

Table 32: Probe results and error values

Figure 26: Displacement field

2.3 Incompressible Elasticity

    2.3.1 Neo-Hookean cylindrical pressure vessel

2.3.1 Neo-Hookean cylindrical pressure vessel

    2.3.1.1 Objective

    2.3.1.2 Geometry Description

    2.3.1.3 External Interactions

    2.3.1.4 Material Description

    2.3.1.5 Method

    2.3.1.6 Solution Characteristics

    2.3.1.7 Results

2.3.1.1 Objective

The main purpose of this test is to demonstrate the accuracy of the neo-Hookean hyperelastic material model in a finite-strain problem with a known non-trivial analytical solution. In particular, we solve the problem of a plane-strain section of a cylinder subjected to internal pressure. This problem was used to verify an academic implementation in Section 3.2 of  [1]. A secondary purpose of this test is to demonstrate the effectiveness of higher-order immersed analysis in a situation where classical low-order finite element analysis would be affected by severe "locking" due to a large bulk modulus.

2.3.1.2 Geometry Description

The geometry for this problem consists of a quarter circular annulus, as shown in figure 27 with inches and inches.

Figure 27: Geometry, loading, and boundary conditions for the cylindrical pressure vessel.

We model it as a three-dimensional object with unit depth to use the three-dimensional solver functionality in Coreform Flex IGA.

2.3.1.3 External Interactions

To model a full cylinder subject to plane-strain kinematics, we apply the symmetry boundary conditions shown in figure 27, and additional symmetry boundary conditions on the faces. The inner surface of the cylinder is subject to a pressure follower load of magnitude . This is treated as given data in the problem setup, but chosen based on an analytical solution to result in a known radial displacement of the inner surface, via equation . The outer surface of the cylinder is traction-free. The symmetry boundary conditions use a dimensionless penalty scale factor of 1000 to perform a controlled study on the accuracy of the constitutive model and potential for locking in high-order splines.

2.3.1.4 Material Description

The neo-Hookean material model is used with the shear and bulk moduli given in table 33.

Shear Modulus

   

Bulk Modulus

   

Table 33: Material properties.

The bulk modulus is selected much larger than the shear modulus to approximate the fully-incompressible case for which a closed-form analytical solution is available, and for which the neo-Hookean constitutive model is widely used in applications. While a pressure stabilization option is available as a beta feature for the nearly-incompressible regime, we focus for now on the standard formulation without any stabilization, which can still produce accurate displacement fields and reaction forces when used with higher-order splines.

2.3.1.5 Method

We immerse the geometry into an axis-aligned background spline and use static continuation in Coreform Flex IGA to solve this problem. We compute results for background splines at three resolutions, which we refer to as "coarse", "medium", and "fine". The coarse mesh has an element size of inches in the and directions, while the medium and fine meshes use elements of size inches and inches, respectively. Element size is held fixed at inches in the direction for all meshes, due to the plane-strain nature of the problem. The coarse spline mesh trimmed to the domain geometry is shown in figure 28.

Figure 28: Coarse trimmed spline discretization.

2.3.1.6 Solution Characteristics

An analytical solution is available from the literature  [1] for the fully-incompressible case, which we are approximating via . For the mesh sizes considered here, the error due to this approximate incompressibility is small relative to the approximation error inherent in using polynomial splines to approximate a transcendental exact solution. The exact solution is most conveniently parameterized in closed form by considering a target radial displacement of the inner surface, , and then calculating the corresponding pressure as where This solution is valid for arbitrarily-large deformations, and we select for calculating the results presented here, corresponding to an applied follower pressure of psi.

2.3.1.7 Results

We plot the computed radial displacement of the inner surface (viz., the component of displacement probed at the point in the reference configuration) as a function of the number of elements through the cylindrical vessel’s wall thickness (computed as ) for linear, quadratic, and cubic splines in figure 29.

Figure 29: Inner-surface radial displacement vs. refinement for different spline degrees.

These results demonstrate that using spline elements of degree can produce accurate displacement solutions, despite the large ratio of bulk modulus to shear modulus. With cubic elements, even the coarsest mesh captures the displacement solution quite accurately. The linear case shows obvious locking, and the pressure solutions for all cases contain large oscillations (not shown here). These issues can be resolved by using a stabilized formulation that is still in beta; the effect of this is illustrated for linear elements in figure 29, where we can see that it makes linear splines more accurate than non-stabilized quadratic splines on coarse meshes. However the locking relief provided by higher-order splines alone is often sufficient for problems in which only deformation and integrated reaction forces are of interest, and pointwise pressure values are not important. Linear splines should not be used in the nearly-incompressible limit without stabilization.

2.4 Plasticity

    2.4.1 J2 plasticity: 3D unit cube with isotropic hardening

2.4.1 J2 plasticity: 3D unit cube with isotropic hardening

    2.4.1.1 Objective

    2.4.1.2 Geometry Description

    2.4.1.3 External Interactions

    2.4.1.4 Material Description

    2.4.1.5 Method

    2.4.1.6 Results

2.4.1.1 Objective

This verification problem benchmarks static elasto-plastic deformation in a single 3D element with isotropic hardening against results published by NAFEMS  [2]. It demonstrates accurate stress predictions for a nonlinear material with higher order spline basis functions.

2.4.1.2 Geometry Description

The problem geometry is a 1x1x1 cube with the vertices labeled as shown in figure 30.

Figure 30: Schematic for the 3D elastoplasticity verification problem.

2.4.1.3 External Interactions

Boundary conditions prescribing displacement values are created for each face of the cube. The type of boundary condition for each face is tabulated in table 34. The prescribed displacement values vary depending on the stage of the simulation and are explained in table table 36.

Location

   

Constraint

Face AEHD

   

Translation along x-axis fixed

Face ABFE

   

Translation along y-axis fixed

Face EFGH

   

Translation along z-axis fixed

Face BCGF

   

prescribed

Face CDHG

   

prescribed

Face ABCD

   

prescribed

Table 34: Boundary conditions for the 3D elastoplasticity verification problem.

2.4.1.4 Material Description

A neo-Hookean isotropic plasticity material model is used with material properties shown in table 35. A linear hardening function is used with the modulus defined in the table.

Property

   

Value

Mass Density

   

0

Young's Modulus

   

2.5e5

Poisson's Ratio

   

0.25

Yield Stress

   

5.0

Isotropic Hardening Modulus

   

6.25e4

Table 35: Material properties.

2.4.1.5 Method

In this test, the cube is meshed with a single 3d element and displacements are calculated using a cubic spline basis. The cube is loaded in tension along a combination of its two axes, in plane strain, before being returned to zero displacement boundary conditions. The loading pattern is shown in table 36.

Step

   

   

   

1

   

R

   

0.0

   

0.0

2

   

2R

   

0.0

   

0.0

3

   

2R

   

R

   

0.0

4

   

2R

   

2R

   

0.0

5

   

2R

   

2R

   

R

6

   

2R

   

2R

   

2R

7

   

R

   

2R

   

2R

8

   

0.0

   

2R

   

2R

9

   

0.0

   

R

   

2R

10

   

0.0

   

0.0

   

2R

11

   

0.0

   

0.0

   

R

12

   

0.0

   

0.0

   

0.0

Table 36: Loading conditions for the 3D elastoplasticity verification problem.

2.4.1.6 Results

Stress results are plotted alongside benchmark values in figure 31 - figure 34. The benchmark values are from results published in the reference text  [2].

Figure 31: XX stress value plot

Figure 32: YY stress value plot

Figure 33: ZZ stress value plot

Figure 34: Von-Mises stress value plot

2.5 Mechanical Contact

    2.5.1 Interference fit of two pipes

2.5.1 Interference fit of two pipes

    2.5.1.1 Objective

    2.5.1.2 Geometry Description

    2.5.1.3 External Interactions

    2.5.1.4 Material Description

    2.5.1.5 Method

    2.5.1.6 Solution Characteristics

    2.5.1.7 Results

2.5.1.1 Objective

Interference fits are a common technique used to join parts through friction by requiring deformation in order to fit the parts together and inducing initial stress and contact forces. In analysis, the parts are initially undeformed and they overlap at the interface. A robust contact algorithm is required to resolve the overlap and produce an accurate initial displacement.

Additionally, it is very common for modeling errors to create unintentional small initial penetrations that can lead to erroneous behavior in an analysis if left unresolved. Manually resolving initial penetrations can be significantly time consuming, so it is important that a contact algorithm is able to automatically detect and resolve initial penetrations.

This current problem demonstrates the ability of Coreform IGA to accurately resolve initial penetrations. The geometry, boundary conditions and loading of this problem create a small deformation, plane strain axi-symmetric setting that allows for a relatively simple closed form analytical solution against which the analysis results can be validated.

2.5.1.2 Geometry Description

The geometry for this problem consists of two concentric quarter cylinder parts as shown in figure 35.

Figure 35: Geometry for the two pipe interference fit problem

The dimensions are specified in table 37. The thicknesses and radii are defined so that there is an interference equal to half the thickness.

Parameter

   

Value(m)

   

   

   

   

   

Table 37: Geometry parameter values

2.5.1.3 External Interactions

Symmetry boundary conditions are enforced in the and directions. Specifically, the faces of each part are constrained to have no displacement in the direction, and the faces are constrained to have no displacement in the driction. Additionally, the and faces of both parts are constrained to have no displacement in the direction to restrict rigid body translation.

The inner surface of the larger part and the outer surface of the smaller part are included in a surface to surface contact definition.

We use a dimensionless penalty scaling factor of for all displacement boundary conditions to ensure accuracy of the applied constraints. This high penalty scaling factor is used for this problem because we are comparing the computed solution to a known reference solution. In most practical applications the default penalty scaling values would produce a solution this is good enough.

2.5.1.4 Material Description

A linear elastic material model is used for both bodies, with the material properties given in table 38.

Young's Modulus

   

Poisson's ratio

   

Table 38: Material properties.

2.5.1.5 Method

We approximate the displacements in each cylinder using quadratic trimmed splines defined on a background mesh with element size of 0.03. We use a implicit statics time stepper based on numerical continuation to solve the nonlinear problem in multiple steps.

The resolution of interference fit relies on the regular contact algorithm to find the equilibrium configuration by removing the initial penetration under the quasi-static condition. See the theory manual for details. For this problem we use an initial activate distance of 0.04 and a ramp time interval of .

We define probes to measure the displacement in each coordinate direction at the points and . These point correspond to locations halfway long the cylinder width, at the outer radius of the larger cylinder and the inner radius of the smaller cylinder. The displacement at these points is compared to an analytical solution. We define an additional probe that measures the contact force that develops due to the contact defined between the two cylinder parts. The value measured is compared to an analytical solution.

2.5.1.6 Solution Characteristics

The analytical solution is derived in the plane-strain axisymmetric setting under small deformation. In the plane-strain axisymmetry, only the radial displacement component, , is nonzero. For strain, likewise, only the radial, , and circumferential strain, , are nonzero, where denotes the radial position. In the absence of Poisson’s ratio and a body force, the equilibrium in a single body is represented by the following equation:

Note that it is not necessary to have a Dirichlet boundary condition even for a static problem, as the plane-strain axisymmetric condition provides enough constraints.

The solution for the inner pipe, , and the one for the outer pipe, , that satisfy the above equation are

with the unknown constants , , and .

Next, we apply the traction boundary and interface conditions. where are defined in table 39, denotes the contact pressure, and denotes the displacement at the inner surface of the outer pipe, i.e. the interface displacement of the outer pipe.

Body

   

Inner Radius

   

Outer Radius

Inner pipe

   

   

Outer Pipe

   

   

Table 39: Pipe dimensions

By solving the above six linear equations, we obtain the six unknowns , , , , , and .

2.5.1.7 Results

Value

   

Measured value

   

Analytical solution

   

Relative error

Inner pipe displacement

   

   

   

2.8 %

Outer pipe displacement

   

   

   

1.1 %

Reaction force X

   

   

   

1.8 %

Reaction force Y

   

   

   

1.8 %

Table 40: Results comparison

Figure 36: X displacement field after interference resolution

    contents  ← prev  up  next →