COMSOL 6.4 - Torsion of a Circular Membrane

The numerical treatment of thin structures with a membrane theory is much simpler than with a shell theory due to the assumption of zero bending stiffness. However, for some load cases, this assumption is disadvantageous. For instance, when a membrane is subjected to compressive stresses, it may trigger wrinkling which is an equilibrium instability. This undesirable limitation can be overcome with the incorporation of a numerical wrinkling model that removes such instabilities.

In this example, a torque is applied to the inner edge of a circular annulus-shaped membrane while the outer edge is fixed. This causes wrinkling due to torsion. The example has been studied by different researchers (see Ref. 1, Ref. 2, Ref. 3, and Ref. 4). The results are similar albeit with some differences. The mesh pattern and discretization used in these research papers are different even though the same geometry, material, and boundary conditions are used. In this example, the model is set up in four different ways with respect to mesh pattern and discretization order, corresponding to each of the research papers.

Model Definition

The geometry consists of a thin circular membrane with an outer radius of 12.5 m and an inner radius of 5 m. it is fixed on the outer edge. A prescribed rotation of 10 degrees is applied on the inner edge.

The linear elastic material is modeled using both isotropic and orthotropic properties. The linear stress-strain relationship defined in Ref. 1 on compliance form reads

(1)

For the isotropic case, Young’s modulus is E = 100 kPa and Poisson’s ratio is ν = 0.3. Note that in Ref. 4, Poisson’s ratio is reported as ν = 0.45. This is assumed to be a typing mistake, since they refer to Ref. 1 for material properties.

In the orthotropic case, the Young’s moduli are E11 = 100 kPa and E22=1000 kPa, Poisson’s ratio is ν12 = 0.3, and the shear modulus is G12 = 38.5 kPa. The elements of the elasticity matrix in the linear stress strain relationship σ = C ε are given in Table 1.

Material Properties

Table 1: Material properties.
Property	Variable	Isotropic Model	Orthotropic Model
Elasticity matrix, 11 element	C11	109.89 kPa	100.91 kPa
Elasticity matrix, 12 element	C12	32.967 kPa	30.272 kPa
Elasticity matrix, 22 element	C22	109.89 kPa	1000.91 kPa
Elasticity matrix, 33 element	C33	76.923 kPa	77 kPa

The constitutive relation for a linear elastic material used in COMSOL Multiphysics is

where the elements of D are related to the elements in C as

It is assumed that the stiffness in the third direction is the same as in the second direction, since the out-of-plane data is not available in Ref. 1. The elements D33, D55, and D66, cannot be equal to zero because the Membrane interface includes degrees of freedom for the out-of-plane strains. However, the choice of out-of-plane data has no effect on the results.

Although the material data in principle can be entered using the or options in the subnode, the elasticity matrix computed in COMSOL Multiphysics does not match the values given in Ref. 1 even for identical Young’s modulus, Poisson’s ratio, and shear modulus. The reason for the discrepancy is the difference between three-dimensional and two-dimensional membrane theory. The authors in Ref. 1 used two-dimensional membrane theory where the thickness direction is not part of the constitutive model, while COMSOL Multiphysics employs a full three-dimensional representation where the membrane thickness enters explicitly. This gives a different elasticity matrix for the same material data. In order to match the values in Ref. 1, the elasticity matrix is first computed, and then these values are entered under the option in the node.

Finite Element Mesh

In Ref. 1 and Ref. 4, three-noded triangular mesh elements are used, but the mesh patterns differs. In Ref. 3, four-noded quadrilateral elements are used with 120 mesh elements, while Ref. 2 uses nine-noded quadrilateral elements with 180 mesh elements. This gives four different mesh pattern to investigate, as shown in Figure 1, Figure 2, Figure 3, and Figure 4.

Figure 1: The triangular mesh with 240 mesh elements from Ref. 1.

Figure 2: The triangular mesh with 240 elements from Ref. 4.

Figure 3: The quadrilateral mesh with 120 elements from Ref. 3.

Figure 4: The quadrilateral mesh with 180 elements from Ref. 2.

Results and Discussions

In all the following figures, four plots corresponding to four combinations of mesh pattern and discretization are presented. The order of the plots is the following and matches the legends position when multiple legends are used:

•

Upper left — Three-noded triangular elements as shown in Figure 1

•

Lower left — Three-noded triangular elements as shown in Figure 2

•

Upper right — Four-noded quadrilateral elements as shown in Figure 3

•

Lower right — Nine-noded quadrilateral elements as shown in Figure 4

Another point to note is that in all the following figures, and in the plot quality section are set to and respectively, in order to represent the computed values as close as possible. Using smoothing would give significantly better plots, but it would hide some interesting features in some of the solutions.

Figure 5 and Figure 6 show the wrinkled regions for isotropic and orthotropic materials, respectively. For the isotropic case, the wrinkled region is symmetric and spreads in almost all elements except the outermost ones, close to the fixed outer boundary. For the orthotropic case, the wrinkled region is not symmetric, and it is localized near the inner edge. Depending on the selected mesh pattern and discretization, the wrinkling region varies significantly.

The first principal stress and the tensile direction for the isotropic and orthotropic cases are shown in Figure 7 and Figure 8, respectively. Figure 9 and Figure 10 show the second principal stress along the wrinkling direction. The lowest value of the second principal stress is almost zero, as is expected since the membrane cannot sustain compressive stresses. The principal stress pattern and directions match those presented in Ref. 1.

The maximum values of the first principal stress obtained with COMSOL Multiphysics for the different cases are given in Table 2 and Table 3. Some papers did not provide the values of maximum first principal stress explicitly. These are marked as Not Available (NA). When only the maximum value of the first principal stress is concerned, the worst match between reference values and COMSOL values is found for the three-noded triangular elements presented in Ref. 4, while a good match is found for three-noded triangular elements presented in Ref. 1. The COMSOL values for the mesh pattern and discretizations presented in Ref. 2 and Ref. 3 are in line with target values.

It is clear that the results are highly mesh and discretization dependent. In particular, the linear elements are performing badly in this example. The four-noded quadrilateral elements show a triangular stress pattern similar to the three-noded triangular elements, see Figure 7 and Figure 8. It is clear that the nine-noded quadrilateral elements are performing better considering the obtained maximum values as well as the stress and wrinkling pattern.

Table 2: Maximum Principal Stress for Isotropic Material.
From	σp1 from Ref. 1	σp1 from Ref. 2	σp1 from Ref. 3	σp1 from Ref. 4
Reference	26.9E3	26E3	NA	25.9E3
COMSOL	25.90E3	22.88E3	26.9E3	17.48E3

Table 3: Maximum Principal Stress for Orthotropic Material.
From	σp1 from Ref. 1	σp1 from Ref. 2	σp1 from Ref. 3	σp1 from Ref. 4
Reference	270.5E3	NA	NA	230E3
COMSOL	235.7E3	247E3	296.5E3	163.4E3

The maximum values of wrinkling measure obtained with COMSOL Multiphysics are compared to the reference values in Table 4. The agreement between reference values and computed value are good.

Table 4: Maximum Values of Wrinkling Measure for Isotropic Material.
From	β in material frame from Ref. 3
Reference	0,19118
COMSOL	0.20824

The choice of mesh pattern shown in Figure 2 seems not desirable for this example. The physics is symmetric while the mesh is antisymmetric. If the counterclockwise torque (same as all references) is replaced by a clockwise torque, then the results are the same for all mesh and discretization patterns except the one used in Ref. 4. With clockwise torque the maximum first principal stress for isotropic and orthotropic material reached 26.03 kPa and 270.23 kPa, respectively. This can be compared with 17.48 kPa and 163.4 kPa obtained with counterclockwise torque. It can be noted that these values match what is reported in Ref. 4 much better.

Figure 5: Wrinkled region with isotropic material for different mesh types.

Figure 6: Wrinkled region with orthotropic material for different mesh types.

Figure 7: First principal stress with isotropic material for different mesh types.

Figure 8: First principal stress with orthotropic material for different mesh types.

Figure 9: Second principal stress with isotropic material for different mesh types.

Figure 10: Second principal stress with orthotropic material for different mesh types.

Notes About the COMSOL Implementation

A wrinkling model based on the modified deformation gradient is incorporated within the membrane theory using the feature, which solves a set of nonlinear equations using the Newton–Raphson method.

Since the unstressed membrane does not have stiffness in the normal direction, the feature is added in order to stabilize the model.

References

1. D.G. Roddeman, “Finite element analysis of wrinkling membranes,” Commun. Numer. Methods Eng., vol. 7, pp. 299–307, 1991.

2. K. Lu, M. Accorsi, and J. Leonard, “Finite element analysis of membrane wrinkling,” Int. J. Numer. Methods Eng., vol. 50, pp. 1017–1038, 2000.

3. H. Schoop, L. Taenzer, and J. Homig, “Wrinkling of nonlinear membranes,” Comput. Mech., vol. 29, pp. 68–74, 2002.

4. A. Jarasjarungkiat, R. Wuchner, and K.U. Bletzinger, “A wrinkling model based on material modification for isotropic and orthotropic membranes,” Comput. Methods Appl. Mech. Eng., vol. 197, pp. 773–788, 2008.

Structural_Mechanics_Module/Buckling_and_Wrinkling/membrane_torsion

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

Model Parameters

In the window, under click .

In the window for , type Model Parameters in the text field.

Locate the section. Click

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

Isotropic Material Properties

In the toolbar, click

and choose > .

In the window for , type Isotropic Material Properties in the text field.

Locate the section. Click

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

Orthotropic Material Properties

In the toolbar, click

and choose > .

In the window for , type Orthotropic Material Properties in the text field.

Locate the section. Click

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

For this model four different mesh patterns are used. Import the meshed geometries to save modeling time.

Triangular Mesh, Pattern 1

In the window, under click .

In the window for , type Triangular Mesh, Pattern 1 in the text field.

Import 1

In the toolbar, click

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file membrane_torsion_tria_mesh1.mphbin.

Click

Triangular Mesh, Pattern 2

In the toolbar, click and choose .

In the window for , type Triangular Mesh, Pattern 2 in the text field.

In the toolbar, click

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file membrane_torsion_tria_mesh2.mphbin.

Click

Finalize

In the window, right-click and choose .

Quadrilateral Mesh, Pattern 1

In the toolbar, click and choose .

In the window for , type Quadrilateral Mesh, Pattern 1 in the text field.

Import 1

In the toolbar, click

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file membrane_torsion_quad_mesh1.mphbin.

Click

Finalize

In the window, right-click and choose .

Quadrilateral Mesh, Pattern 2

In the toolbar, click and choose .

In the window for , type Quadrilateral Mesh, Pattern 2 in the text field.

Import 1

In the toolbar, click

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file membrane_torsion_quad_mesh2.mphbin.

Click

Finalize

In the window, right-click and choose .

The first study runs with linear shape order of the elements. Change the default discretization to .

Membrane (mbrn)

In the window, under click .

In the window for , click to expand the section.

From the list, choose .

Linear Elastic Material 1

In the window, under > click .

In the window for , locate the section.

From the list, choose .

Wrinkling 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

Thickness and Offset 1

In the window, under > click .

In the window for , locate the section.

In the d0 text field, type th.

Fixed Constraint 1

In the toolbar, click

and choose .

Select Edges 2, 3, 7, and 11 only.

Prescribed Displacement 1

In the toolbar, click

and choose .

Select Edges 4, 5, 8, and 10 only.

In the window for , locate the section.

From the list, choose .

In the u 0 x text field, type (R11-1)*X+R12*Y+R13*Z.

From the list, choose .

In the u 0 y text field, type R21*X+(R22-1)*Y+R23*Z.

From the list, choose .

In the u 0 z text field, type R31*X+R32*Y+(R33-1)*Z.

Add a feature, in order to suppress the out-of-plane singularity of the unstressed membrane.

Stabilization 1

In the toolbar, click

and choose .

The fourth study runs with quadratic shape order of the elements, so add a node. To do so, you first need to activate advanced physics options.

Click the

button in the toolbar.

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

Click .

Quadratic Discretization

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

In the text field, type Quadratic Discretization.

Add a node in order to run the same study with different materials.

Materials

Material Switch 1 (sw1)

In the window, under right-click and choose > .

Isotropic Material

In the window, right-click and choose .

In the window for , type Isotropic Material in the text field.

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


Property	Variable	Value	Unit	Property group
Elasticity matrix	{D11, D12, D22, D13, D23, D33, D14, D24, D34, D44, D15, D25, D35, D45, D55, D16, D26, D36, D46, D56, D66} ; Dij = Dji	{DD11_iso, DD12_iso, DD22_iso, 0, 0, DD33_iso, 0, 0, 0, DD44_iso, 0, 0, 0, 0, DD55_iso, 0, 0, 0, 0, 0, DD66_iso}	Pa	Anisotropic
Density	rho	0	kg/m³	Basic

Orthotropic Material

In the window, right-click and choose .

In the window for , type Orthotropic Material in the text field.

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


Property	Variable	Value	Unit	Property group
Elasticity matrix	{D11, D12, D22, D13, D23, D33, D14, D24, D34, D44, D15, D25, D35, D45, D55, D16, D26, D36, D46, D56, D66} ; Dij = Dji	{DD11_orth, DD12_orth, DD22_orth, 0, 0, DD33_orth, 0, 0, 0, DD44_orth, 0, 0, 0, 0, DD55_orth, 0, 0, 0, 0, 0, DD66_orth}	Pa	Anisotropic
Density	rho	0	kg/m³	Basic