Biagini Nava-Master Thesis
Biagini Nava-Master Thesis
Biagini Nava-Master Thesis
1 Introduction 3
1.1 Cables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2 Membranes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.3 Form finding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.4 Programming language and FEM software . . . . . . . . . . . . . . . . . . . 8
1.4.1 Python and PyCharm . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.4.2 GMSH . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
1.4.3 CalculiX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.4.4 Code_Aster . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.4.5 Dlubal RFEM5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.5 Examples of tensile structures . . . . . . . . . . . . . . . . . . . . . . . . . . 11
i
5.3.1 Newton-Raphson method . . . . . . . . . . . . . . . . . . . . . . . . 44
5.3.2 Newton-Krylov method . . . . . . . . . . . . . . . . . . . . . . . . . 47
5.4 Constitutive models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
5.4.1 Isotropic liner elastic material . . . . . . . . . . . . . . . . . . . . . 49
5.4.2 Orthotropic linear elastic material . . . . . . . . . . . . . . . . . . . 51
5.4.3 Isotropic elasto-plastic material . . . . . . . . . . . . . . . . . . . . . 52
5.4.3.1 One dimensional elasto-plastic constitutive law . . . . . . . . 52
5.4.3.2 Three dimensional elasto-plastic constitutive law . . . . . . . 54
6 Cutting pattern 59
6.1 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
6.1.1 Hypar membrane . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
6.1.2 Hexagonal membrane . . . . . . . . . . . . . . . . . . . . . . . . . . 63
ii
8.3.4 Post-processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
8.3.5 Dlubal comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
8.3.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
8.4 Plastic cable under a concentrated load . . . . . . . . . . . . . . . . . . . . . 121
8.4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
8.4.2 Pre-processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
8.4.3 Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
8.4.4 Post-processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125
8.4.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125
8.5 Cable net . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
8.5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
8.5.2 Pre-processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
8.5.3 Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
8.5.4 Post-processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
8.5.5 Dlubal comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
8.5.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
8.6 Plane membrane . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142
8.6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142
8.6.2 Pre-processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142
8.6.3 Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150
8.6.4 Post-processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151
8.6.5 Dlubal comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . 153
8.6.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153
8.7 Plane orthotropic membrane . . . . . . . . . . . . . . . . . . . . . . . . . . . 168
8.7.1 Intoduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 168
8.7.2 Pre-processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 168
8.7.3 Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 176
8.7.4 Post-processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 177
8.7.5 Dlubal comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . 177
8.7.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 177
8.8 Plastic plane membrane . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 194
8.8.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 194
8.8.2 Pre-processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 194
8.8.3 Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 202
8.8.4 Post-processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 203
8.8.5 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 203
8.9 Plane membrane with collaborating cables . . . . . . . . . . . . . . . . . . . 220
8.9.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 220
8.9.2 Pre-processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 220
8.9.3 Solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 228
8.9.4 Post-processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 229
8.9.5 Dlubal comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . 229
8.9.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 246
iii
9.1.3 Input parsing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 250
9.1.4 Output parsing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 255
9.1.5 Example problem definition . . . . . . . . . . . . . . . . . . . . . . 258
9.1.6 Test functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 259
9.1.7 Example run . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 261
9.1.8 Auxiliary functions . . . . . . . . . . . . . . . . . . . . . . . . . . . 261
9.1.8.1 Lenght as norm . . . . . . . . . . . . . . . . . . . . . . . . 261
9.1.8.2 Area with Erone’s formula . . . . . . . . . . . . . . . . . . 262
9.1.8.3 Normal to the TRIA3 element . . . . . . . . . . . . . . . . . 263
9.1.8.4 Nodal forces for the membrane . . . . . . . . . . . . . . . . 265
9.1.9 Dirchlet boundary equations . . . . . . . . . . . . . . . . . . . . . . 265
9.1.10 Free equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 267
9.1.10.1 Cables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 267
9.1.10.2 Membrane . . . . . . . . . . . . . . . . . . . . . . . . . . . 268
9.1.11 Wrapper for the objective function . . . . . . . . . . . . . . . . . . . 269
9.1.12 Form-finding objective function . . . . . . . . . . . . . . . . . . . . 271
iv
10.3.6 Orchestrating scripts . . . . . . . . . . . . . . . . . . . . . . . . . . 318
10.3.7 Form-finding geometry . . . . . . . . . . . . . . . . . . . . . . . . . 321
10.3.8 Modified deck . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 324
10.3.9 Post-processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 327
10.3.10 Dlubal comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . 327
10.3.11 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 331
10.4 Plane membrane with cables on the borders . . . . . . . . . . . . . . . . . . . 336
10.4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 336
10.4.2 Test definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 336
10.4.3 Pre-processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 336
10.4.4 Problem definition . . . . . . . . . . . . . . . . . . . . . . . . . . . 341
10.4.5 Solution deck . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 342
10.4.6 Orchestrating scripts . . . . . . . . . . . . . . . . . . . . . . . . . . 343
10.4.7 Form-finding geometry . . . . . . . . . . . . . . . . . . . . . . . . . 345
10.4.8 Modified deck . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 353
10.4.9 Post-processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 356
10.4.10 Dlubal comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . 356
10.4.11 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 360
10.5 Hypar membrane with cables on the borders . . . . . . . . . . . . . . . . . . 365
10.5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 365
10.5.2 Test definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 365
10.5.3 Pre-processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 365
10.5.4 Problem definition . . . . . . . . . . . . . . . . . . . . . . . . . . . 370
10.5.5 Solution deck . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 371
10.5.6 Orchestrating scripts . . . . . . . . . . . . . . . . . . . . . . . . . . 372
10.5.7 Form-finding geometry . . . . . . . . . . . . . . . . . . . . . . . . . 375
10.5.8 Modified deck . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 378
10.5.9 Post-processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 381
10.5.10 Dlubal comparison . . . . . . . . . . . . . . . . . . . . . . . . . . . 381
10.5.11 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 385
v
12.1.8 Post-processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 438
12.1.9 Tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 438
12.2 Topological optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 441
12.2.1 Serialization class . . . . . . . . . . . . . . . . . . . . . . . . . . . . 441
12.2.2 Wrapper for the execution . . . . . . . . . . . . . . . . . . . . . . . 442
12.2.3 Optimization objective function . . . . . . . . . . . . . . . . . . . . 443
12.2.4 Problem setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 444
12.2.5 Mesh pre-processing . . . . . . . . . . . . . . . . . . . . . . . . . . 446
12.2.6 The deck in code_aster . . . . . . . . . . . . . . . . . . . . . . . . . 447
12.2.7 Example run of the program . . . . . . . . . . . . . . . . . . . . . . 452
12.2.8 Post-processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 453
12.2.9 Tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 453
Bibliography 457
Index 461
vi
List of Figures
4.1 Displacement diagram for a 1 dof system for the DRM - Edited from [Wiki411] 34
1
Form-Finding and Mechanical Behaviour of Cable and Membrane Structures
7.1 Displacement shape obtained varying the external load for a catenary cable. . 70
7.2 Displacement shape obtained varying the pretension for a catenary cable. . . . 71
7.3 Geometric quantities for the parabola cable - [Mal117] . . . . . . . . . . . . . 73
7.4 Reaction forces for the cable a under a distributed load. . . . . . . . . . . . . . 78
7.5 Tension in the cable under a distributed load. . . . . . . . . . . . . . . . . . . 79
7.6 Similarities between the parabola and catenary cables. . . . . . . . . . . . . . 80
8.1 Total displacements of the catenary cable under its own weight. . . . . . . . . 89
8.2 Displacements in the x and y-direction of the catenary cable. . . . . . . . . . . 89
8.3 Axial stress and elongation of the catenary cable. . . . . . . . . . . . . . . . . 90
8.4 Displacement of the mid-point of the catenary cable. . . . . . . . . . . . . . . 91
8.5 Axial stress and elongation of the mid-point of the catenary cable. . . . . . . . 92
8.6 Horizontal and vertical reactions of the catenary cable. . . . . . . . . . . . . . 93
8.7 Check of the displacements and tension in the cable. . . . . . . . . . . . . . . . 94
8.8 Check of the axial elongation of the cable. . . . . . . . . . . . . . . . . . . . . 95
8.9 Total displacements of the parabola cable under a distributed load. . . . . . . . 101
8.10 Displacements in the x and y-direction of the parabola cable. . . . . . . . . . . 101
8.11 Axial stress and axial elongation of the parabola cable under a distributed load. 102
8.12 Displacement of the mid-point of the parabola cable. . . . . . . . . . . . . . . 103
8.13 Axial stress and elongation of the mid-point of the parabola cable. . . . . . . . 104
8.14 Horizontal and vertical reactions of the the supports of the parabola cable. . . 105
8.15 Check of the displacements and tension in the cable. . . . . . . . . . . . . . . . 107
8.16 Check of the axial elongation of the cable. . . . . . . . . . . . . . . . . . . . . 108
2 List of Figures
Form-Finding and Mechanical Behaviour of Cable and Membrane Structures
List of Figures 3
Form-Finding and Mechanical Behaviour of Cable and Membrane Structures
4 List of Figures
Form-Finding and Mechanical Behaviour of Cable and Membrane Structures
8.107Normal stresses in the x and y-direction of the mid-point of the cabled membrane.234
8.108Normal stresses in the x and y-direction of the pin of the cabled membrane. . . 235
8.109Normal strains in the x and y-direction of the center point of the cabled membrane.236
8.110Normal strains in the x and y-direction of the mid-point of the cabled membrane. 237
8.111Normal strains in the x and y-direction of the pin of the cabled membrane. . . . 238
8.112Von-Mises stresses and strains of the center point of the membrane with cables. 239
8.113Von-Mises stresses and strains of the mid-point of the membrane with cables. . 240
8.114Von-Mises stresses and strain of the pin of the membrane with cables. . . . . . 241
8.115Check of the deformations and the 𝜎x of the membrane with cables. . . . . . . 242
8.116Check of the 𝜎y and the 𝜎mises of the membrane with cables. . . . . . . . . . . 243
8.117Check of the 𝜀x and the 𝜀y of the membrane with cables. . . . . . . . . . . . . 244
8.118Check of the 𝜀mises of the membrane and the tension in the cables. . . . . . . . 245
List of Figures 5
Form-Finding and Mechanical Behaviour of Cable and Membrane Structures
10.35Normal strains in the x and y-direction of the hypar membrane with cables. . . 383
10.36Von-Mises stresses and strains of the hypar membrane with cables. . . . . . . . 384
10.37Check of the displacement and the 𝜎x of the hypar membrane with cables. . . . 386
10.38Check of 𝜎y and the 𝜎mises of the hypar membrane with cables. . . . . . . . . . 387
10.39Check of 𝜀x and the 𝜀y of the hypar membrane with cables. . . . . . . . . . . . 388
10.40Check of 𝜀mises and of the internal action in the cable of the hypar membrane
with cables. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 389
11.1 Topology for the form-finding of a hexagonal membrane with cables. . . . . . . 393
11.2 Total displacements of the hexagonal membrane under the live loads. . . . . . . 408
11.3 Displacements in the x and y-direction of the hexagonal membrane. . . . . . . 408
11.4 Normal stresses in the x and y-direction of the hexagonal membrane. . . . . . . 409
11.5 Normal strains in the x and y-direction of the hexagonal membrane. . . . . . . 409
11.6 Von-Mises stresses and strains of the hexagonal membrane under the live loads. 409
11.7 Displacements for the front and middle membrane of the hex structure. . . . . . 410
11.8 Normal stresses in the x and y-direction of the front membrane. . . . . . . . . . 411
11.9 Normal stresses in the x and y-direction of the middle membrane. . . . . . . . . 412
11.10Normal strains in the x and y-direction of the front membrane. . . . . . . . . . 413
11.11Normal strains in the x and y-direction of the middle membrane. . . . . . . . . 414
11.12Von-Mises stresses and strains of the front membrane of the hexagonal structure. 415
11.13Von-Mises stresses and strains of the middle membrane of the hexagonal struc-
ture. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 416
11.14Check of the form-finding: shape and 𝜎x . . . . . . . . . . . . . . . . . . . . . 418
11.15Check of the form-finding: 𝜎y and 𝜎mises . . . . . . . . . . . . . . . . . . . . . 419
11.16Check of the form-finding: 𝜀x and 𝜀y . . . . . . . . . . . . . . . . . . . . . . . 420
11.17Check of the form-finding: 𝜀mises and the tension in the cables . . . . . . . . . 421
11.18Check after the application of a live load:: displacements and 𝜎x . . . . . . . . 422
11.19Check after the application of a live load: 𝜎y and 𝜎mises . . . . . . . . . . . . 423
11.20Check after the application of a live load: 𝜀x and 𝜀y . . . . . . . . . . . . . . . 424
11.21Check after the application of a live load: 𝜀mises and the tension in the cables . 425
6 List of Figures
Abstract
This master thesis is intended to inspect taut cable and membrane structures from different
perspectives.
From the structural point of view, the materials and structural elements of tensile structures are
extremely durable and the lightweight nature of membranes and cables ensures an extraordi-
nary efficiency towards ambitious projects, specifically when long spans without columns are
to be designed. The structural stiffness of this elements stems from the applied pretension,
being the flexural and shear stiffness absolutely negligible. As a consequence, out of plane
loads determine large displacements and eventually deformations which force the engineer to
employ the methods of nonlinear analysis even in the presence of linear elastic materials. For
this reason, finite element non-linear analysis is an essential ingredient in the design of this
structures.
Having analysed cables and membranes individually from the production process to the struc-
tural applications, we then concentrate on the mechanical properties of the membranes. We in-
troduce different constitutive law models such as both isotropic and orthotropic, combined with
both linear elastic and elasto-plastic behaviour. These choices have been developed due to the
composition of the fabrics which consist of two elements: the warp, which are straight threads
weaved with the weft that runs from one side to the other, determining a typical orthotropic
behaviour; all covered with an impermeable plastic material (coating) that also guarantees the
possibility to connect multiple fabric portions together.
Once the material and the structural properties have been exhaustively investigated, we focus on
the so-called form-finding, which is a pre-processing procedure intended to the define the shape
of the structure under a prescribed set of pretension loads and kinematic constrains. We also
underline how the prescription of a stress state of uniform tension in the cables and membranes
guarantees a better performance under both dead and live loads. At the same time, we stress
out how the form-finding process allows the design of beautiful architectural shapes, avoiding
the membranes to becoming slack when slightly detensioned.
As final step towards a complete analysis of taut structures, we highlight the limits of the tra-
ditional form-finding procedure which cannot perform an actual structural optimization in the
case of nonlinear constitutive laws or in presence of a significant interaction with the primary
structure. We therefore introduce a new generalized finite element approach based on the con-
strained optimization of a convex function which defines the target for one or more specific
static or geometric quantities.
Form-Finding and Mechanical Behaviour of Cable and Membrane Structures
2 List of Figures
CHAPTER 1
Introduction
Pure taut structures are made mainly with cables combined with membranes. This kind of ele-
ments allow to develop extremely various typologies of structures and architectures providing
many benefits and advantages. First, under the aesthetics point of view, the flexible character-
istics of the membrane allows to design almost unlimited forms and shapes and their unique
translucency, that offers soft diffused natural light, transform this kind of structures in authentic
artworks. On the other hand, as far as the structural aspects are concerned, the excellent dura-
bility of the materials and structural elements involved, for instance ETFE film, PVC or PTFE
fiberglass, combined with the lightweight nature of membranes and cables, make this kind of
structures an extremely efficient solution towards ambitious projects when long spans without
supports are to be designed. Moreover, a great amount of recurring elements are prefabricated,
therefore, the construction process is much easier and faster if compared to an equivalent-sized
conventional building, thanks to the light weight but stiff nature of membrane and cable el-
ements. The structural stiffness of this elements come from the pretension applied on them,
in fact they are not able to bear bending moments or compression of any kind. Finally, ten-
sile architecture determines some important environmental benefits, since the aforementioned
translucency reduces the daily costs for lighting and the high sun reflectivity and low absorption
determine a significant drop in the electrical energy required by the building.
1.1 Cables
During the sixties and early seventies, it was common thought that the architectural and eco-
nomic potential of cable roof structures would lead to an increasing demand for this type of
building. Consequently, a great deal of work was carried out world-wide to study the behaviour
of different types of structural system. The interest in cable structures also stimulated an in-
creasing interest in the use and development of numerical methods for solving large systems
of nonlinear equations, made possible by the emergence of the high-speed electronic computer
[Buc103].
3
Form-Finding and Mechanical Behaviour of Cable and Membrane Structures
Cables, as structural elements, are characterized by great flexibility and negligible bending
stiffness. However, they are able to carry transversal loads if properly pre-stressed, thus being
able to span very long distances with minimal cross-sectional area. Cables are a construction of
wire strands, laid (i.e. twisted) helically around a core, to form a tension member of symmetri-
cal cross section. Cables have a high strength-to-weight ratio as the wire strands are drawn to
high strengths and laid to carry tensile loads efficiently.
The construction and subsequent mechanical properties are determined by the wire pattern and
wire size used for the core and outer strands. Generally speaking, cables can be classified as
[Tri409]:
• Wire Rope:
For a given strength, the wire rope is larger in diameter and lower in stiffness than structural
strands or solid rods, so it may not be the most cost-effective product for static, tension load-
carrying structural members, but it is commonly used for “flexible structures”. However, in
cases where a relatively high stretch is desired, or where the cable are used with pullies wire
ropes may be the appropriate choice.
• Structural Strand:
4 Chapter 1. Introduction
Form-Finding and Mechanical Behaviour of Cable and Membrane Structures
The structural strand is the cable construction most commonly used in structural applications.
It offers an economical combination of strength and stiffness for structures. As cables become
larger in diameter, the number of strands increases.
• Full Locked Cable:
A variation for larger galvanized cables is the “full locked cable”. These cables have their
outer strands drawn in a “Z” shape so that they interlock and form a smoother outside layer
of strands. The interlocking strands also yield a denser cross section and therefore a higher
effective elastic modulus.
1.2 Membranes
The membranes taken here into consideration are fabrics which are flat, thin and with a flexible
surface created by a network of fibres perpendicular to each other according to a plain weave
pattern. Thanks to these characteristics they are suitable where large areas have to be covered.
The membrane is a composite material consisting of fabric and coating. The coating has various
roles: it protects the woven fabric against moisture, UV radiation, attack by fungi, etc. But it
also guarantees water-tightness and stabilizes the fabric configuration, giving shear stiffness
to the membrane. Furthermore, it makes possible to join two membrane portions (welded or
glued seams, etc.), with gradual load transfer between the two fabrics. The simplest possible
constitutive model that can be adopted for membrane is a plane stress elastic orthotropic model.
Such modelling is still used in the design offices even if
1.2. Membranes 5
Form-Finding and Mechanical Behaviour of Cable and Membrane Structures
Due to the low compressive and bending stiffness of their surface, membrane constructions only
withstand tensile stresses and the plane surfaces react sensitively to the deformation imposed
by external loads.
If subjected to vertical loads as snow or water, they are prone to flutter and create zones of
load concentrations, which in extreme cases can lead to the failure of the construction. Thus,
the stabilization of plane areas requires high pre-stressing and for this reason, double curved
constructions are more commonly used. In addition to the tensile stresses, tangential stresses
are developed inside the membrane surface, which contribute to support the applied load. The
combined action of traction and shear gives rise to a surface structure that can carry the loads
according to a funicular mechanism, at least as long as the applied loads do not cause the
occurrence of compressive stresses, which are caused by means of folds in the membrane.
Therefore, the fundamental principle underlying the design of membrane structures is that the
surface must be maintained in a state of traction under any conditions.
The conceptual design and realisation of taut constructions are based, at least in the first phases,
on different principles with respect to more standard building structures. In fact, the more
evident difference is that membranes and taut structures are developed in general according to
self-regulations, engineers are not guided by any normative during the design process. Another
peculiar difference is that the final shape is extremely influenced by the boundary conditions
assigned and the load to which it is subjected [Meh408].
We shall now describe the two main types of membrane structures mentioned before:
the boundary tensioned structures and the pneumatic structures.
This is a particular type of structure with a thin, flexible surface that carries loads primarily
through tensile stresses. Taut membrane structures are characterized by a shape that is deter-
mined by tension in the textile material and the geometry of the support structure (boundary
conditions). The stress demands turn to be a central point in the construction and patterning of
each taut structure element, from the membrane to the anchors.
New applications for an extremely versatile material, the Ethylene TetraFluoroEthylene
(ETFE), have been recently discovered, an acronym which stands for ETFE and is a partially
fluorinated polymer, which is a plastic material that contains fluorine. It has been used in ar-
chitecture since the ’80s due to its exceptional characteristics. In fact, ETFE is transparent like
glass but, compared to it, it is lighter, stronger, and simple and economical to install. Aside
from its lightweight nature, the ETFE is totally permeable to light and UV rays and is totally
recyclable.
6 Chapter 1. Introduction
Form-Finding and Mechanical Behaviour of Cable and Membrane Structures
ETFE cushions consist of a two layers of membrane (enclosing pressured air in between) sup-
ported by a small pressure different between the inside and outside. Since external loads are
usually quite small, the required internal pressures may be small. The air chamber also con-
tributes to the thermal insulation of the system. Along the entire closed perimeter of the area
where the cushion is located, an extruded aluminium frame is placed. The frame is connected
to the main supporting structure through plates and bolts. The cushions are fixed to valves
connected to the pumps of the pressurization system which, once the system is installed, will
start working by inflating the membranes up to the pressure required.
Seldom a field in the construction field requires such a close collaboration between all the com-
ponents involved in design, manufacture and execution, as construction with textile materials.
From the design process to the production and the erection, the collaboration between archi-
tects, engineers and the construction companies becomes unavoidable and essential to obtain
the desired result. The design of tensioned membrane and cables net follows a three-stage
procedure:
• Form-Finding;
• Cutting Pattern (procedure which determine the template from which the parts of a mem-
brane are traced onto fabric before being cut out and assembled);
• Static Analysis, generally highly nonlinear, under external (environmental) loads.
The form-finding represents the first step to be executed, being the pre-processing procedure
leading to the definition of the shape which the structure will attain according to the pretension
assigned to each element. More precisely, the form-finding calculation will provide a surface
geometry (or a length in the case of cables) which is equilibrated with the boundary condi-
tions and the desired pretension provided by the designer. The results obtained are completely
independent from the particular mechanical behaviour of the materials used, but, as already
mentioned, depend only upon the pre-stress magnitude and the boundaries conditions. This
computation is unavoidable due to the fact that, initially, membrane and cables have negligi-
ble stiffness which does not allow to have an initial shape, but it increases with the growth of
the deformations. Moreover, prescribing a uniform tension in the membranes guarantees bet-
ter performances in structural terms, in fact, they are less likely to wrinkle or fail by fatigue
[GaLe206].
This problem can be solved thanks to three different methods:
• Force Density Method (FDM): a procedure used for the design of cable structures;
• Surface Stress Density Method (SSDM): analogous method that considers only shapes
made up by mebranes;
• Dynamic Relaxation Method (DRM): step by step method that determines the equilib-
rium position starting from an initial configuration through a pseudo-dynamic analysis.
As final remark, it must be pointed out that we will refer to the term “mesh”, meaning the
connections between each node of the structure. The initial “mesh file” in the analysis code
must be interpreted as an initial topology, that should not be misunderstood with the initial
configuration, which is actually the aim of the form finding.
The mindset at the base of the development of this master’s thesis was to employ mainly open-
source software with the aim to give evidences of the potentialities and the great hard work
carried out (and freely shared) by enthusiasts programmers all around the world.
The programming code selected is Python, a general-purpose programming code, which can
be used to build almost anything, from a web site to a rocket launcher system. It is supported
by an active community and great corporate sponsors, for example, Google has developed in
Python most of its platforms and applications since 2006. This aspect is extremely relevant, in
fact, if a company like Google wants their development teams, to work with a programming
code, they will create a vast amount of guide, tutorials and tools, which contributes to enlarge
8 Chapter 1. Introduction
Form-Finding and Mechanical Behaviour of Cable and Membrane Structures
exponentially the list of documentation and support. Python has amazing libraries which allow
to save time and reduce the initial development cycle, for example in our code we take advan-
tage of NumPy, the most famous library (with SciPy) for scientific purposes. Furthermore, it
is extremely reliable and efficient, in fact Python can work in almost any environment without
any loss of performances. Last but not least, Python is incredibly easy to learn and use, it is
considered one of the most accessible programming language available, its syntax is close to
the natural language if compared with other code. For this reasons, young developers around
the world are investing their resources and time in learning this fresh and new way of coding.
Equally, we decided to follow this trend.
PyCharm is the IDE (Integrated Development Environment) which was chosen to develop the
python code for solve the Form-Finding. It is an IDE for Python full featured, it is developed
by JetBrains in two versions: The Free Community Edition, freely downloadable, and the com-
mercial one which targets enterprise developers (we used the Community Edition). PyCharm
has also become famous thanks to companies like Twitter, Groupon, Spotify and eBay, which
use it to develop their applications. Most features are available in the free version, such as in-
telligent code completion, particularly intuitive project navigation, error checking/fixing, smart
refactoring, Graphical debugger and test runner.
1.4.2 GMSH
The best description of this program comes directly from the words of its two developers,
Christophe Geuzaine and Jean-François Remacle.
“Gmsh is a free 3D finite element mesh generator with a built-in CAD engine and post-
processor. Its design goal is to provide a fast, light and user-friendly meshing tool with paramet-
ric input and advanced visualization capabilities. Gmsh is built around four modules: geometry,
mesh, solver and post-processing. The specification of any input to these modules is done either
interactively using the graphical user interface, in ASCII text files using Gmsh’s own scripting
language (.geo files), or using the C++, C, Python or Julia API” - [Gms006].
1.4.3 CalculiX
It is a powerful and versatile open-source package developed to solve field problems imple-
menting the finite element method. CalculiX takes advantage of the Abaqus input format,
which allows to use commercial pre-processors, but it is also offers a built-in interactive 3D-
tool which employs the openGL API. The solver has been designed to perform static, dynamic
and thermal analysis, solving both linear and non-linear systems. To emphasise the versa-
tility of this software, we underline that the pre-processor is able to write mesh related data
for Nastran, Abaqus, Ansys, code-aster and for the free CFD codes dolfyn, duns, ISAAC and
OpenFOAM.
The CalculiX package has been developed by a team of engineers employed at MTU Aero
Engines in Munich (Germany), who, motivated by their passion, regularly improve the code in
their spare time.
1.4.4 Code_Aster
The direct competitor of CalculiX is Code_Aster, a free software for numerical simulation
of materials and mechanical structures, developed mainly by the “Analyses Mécaniques et
10 Chapter 1. Introduction
Form-Finding and Mechanical Behaviour of Cable and Membrane Structures
The RFEM structural analysis program is a commercial finite element software which is the
foundation of a modular structural engineering software system. The main program was de-
signed to analyse structures and materials for slabs, shells, beams and rods. Additional modules
can be used to perform further analyses and projects according to different standards. This soft-
ware was chosen because of its module dedicated to the Form-Finding, in order to validate the
code written in Python. The RF-FORM-FINDING add-on module searches for the equilib-
rium shapes of taut structures that are subject to a predefined pretensioning, such as cable and
membrane structures. The shape is calculated through the balance between the surface tension
(prestressing and other loads such as own weight, pressure, etc.) and the computed constrain
reactions. For structural analysis, the new pre-stressing form is then used as in the initial state
for the application of the live loads.
The materialization of the membrane structures can also be determined with the additional RF-
CUTTING-PATTERN module. The cutting pattern module are not analysed in this master’s
thesis, since we focused on the form finding problem solution.
In recent years, taut structures have been widely employed mainly with the aim of roofing
residential and industrial areas due to all the favorable properties previously presented.
Here are listed some interesting examples of structures which implements tensile elements.
12 Chapter 1. Introduction
Form-Finding and Mechanical Behaviour of Cable and Membrane Structures
1.18: F1-YAS Marina Circuit - Abu Dhabi (United Arab Emirates) - 2008 - [Maf406]
14 Chapter 1. Introduction
Form-Finding and Mechanical Behaviour of Cable and Membrane Structures
16 Chapter 1. Introduction
CHAPTER 2
The force density method (FDM) represents an excellent approach for a direct study of the
possible shapes that a pretensioned net of cables can attain. In the initial phases of the project,
since some important data about the structure may often be lacking (for instance, the materials
and/or the cross sections), the FDM procedure could suggest the designer a possible final shape
of the structure, by imposing a certain desired pretension state. Examples of use of this method
can be found in the Mannheim Multihalle (Frei Otto, Germany, 1975) or in the Solemar Therme
(Bad Durrheim, Switzerland 1987) [Bia302].
As developed more in detail later in the thesis, the FDM concept focuses on the linearization
and decoupling of the system of equilibrium equations for the nodes through a simple escamo-
tage: rather than considering the tension assigned in each cable, a force density (for each stretch
of the ropes) is considered, defined as a relationship between the final tension and the deformed
length:
𝑇
𝑞= (2.1)
𝐿
For a system characterised by multiple nodes, as could be a cable net, the equation number
may become extremely large; therefore, it is appropriate to describe the problem with a matrix
formulation. The final coordinates of each node are influenced by the physical quantities (coor-
dinates and density of force) that comes from the immediately adjacent nodes and by the traits
that connect them. Therefore, it is necessary to construct a matrix containing the topological
information of the various nodes which describes also the way in which they are connected.
The incident matrix, C, a extension of the concept of boolean matrix, serves this purpose; it is
17
Form-Finding and Mechanical Behaviour of Cable and Membrane Structures
The resulting matrix is not necessarily symmetrical and has no prerequisites of regularity: it
strongly depends on the numbering of the nodes (that is absolutely arbitrary) and can have
a different number of elements in each column. It is important to underline the fact that the
incident matrix contains only topological and not geometric information, because according to
the FDM, the initial configuration of the system does not influence the final results, because
the initial coordinates of the nodes do not even exist; only the topology, which means how the
nodes are connected and related to each other, is crucial in defining the geometry returned by
the form-finding process.
We introduce now the following vectors whose meaning is clearly self-explanatory:
⎧
⎪
⎨x = [𝑥i ]
y = [𝑦i ] (2.3)
⎪
⎩
z = [𝑧i ]
⎧
⎪
⎨u = [𝑥k − 𝑥i ]
v = [𝑦k − 𝑦i ] (2.4)
⎪
⎩
w = [𝑧k − 𝑧i ]
With 𝑘, 𝑖 = 1, 2, ...𝑛n , the vectors x, y, z containing the nodal coordinates and u, v, w the
distances between adjacent nodes.
It is immediate to understand that the vectors u, v, w originate from the product of C and x, y,
z:
⎧
⎪
⎨u = C x
v = Cy (2.5)
⎪
⎩
w = Cz
Assuming for a certain topology n nodes with 𝑛f of them constrained, and m cables, the size of
the incident matrix will be (𝑚, 𝑛) . It is useful to partition the C-matrix (and consequently also
the coordinate vectors) into two sub-matrices:
• 𝐶N containing the the free nodes;
• 𝐶F containing the fixed nodes.
In order to reduce the computational cost of the problem, we can also write:
⎧
⎪
⎪ C = [CN CF ]
⎪
⎨x = [x x ]
N F
(2.6)
⎪
⎪ y = [yN yF ]
⎪
⎩
z = [zN zF ]
The surface stress density method is introduced to consider, in addition to cables, also taut
membranes, which are particulartly important elements of tensile structures. The purpose of
the chapter is the development of a sound structural theory capable of interpreting the strong
non-linearities induced by the mechanical properties of membranes. The SSDM follows the
same principles of the FDM but allows the combination of study of mixed cable and membrane
structures thanks to the application of a stress density for the area elements in conjunction with
a force density for the linear elements.
The stress state defined by the designer is relative to a system of local coordinates (x,y), typical
of each element and composed by two perpendicular tension: 𝜎x and 𝜎y . Without any loss of
generality, an isotropic stress state with 𝜎x = 𝜎y = 𝜎0 and zero shear is defined, in order to
guarantee that the membrane is nowhere subject to compression or torsion that may sag some
portions, compromising its structural properties against the application of live loads. This
peculiarity has also the advantage of ensuring that the resulting forces that the border elements
of the membrane exert on a node are directed orthogonally to the opposite side (this can be
verified by applying the principle of virtual works explained in section 3.2) as shown in the
following figure 3.1:
The expression of these forces can be obtained from the stress state associated with the element.
Provided that the geometric properties of the element and of the membrane are known, such as
the length of the sides and the thickness:
1
𝑇i = 𝐿i 𝑡𝜎 (3.1)
2
The corresponding matrix formulation is (considering node c):
Tc = 𝑇i n (3.2)
21
Form-Finding and Mechanical Behaviour of Cable and Membrane Structures
with n the versor perpendicular to the side opposite to the considered node. By means of this
versor, it is possible to calculate the coordinates of point H (point of intersection between the
side opposite the node and orthogonal to the side passing through the same node); they can be
calculated knowing the positions of the two extreme points of the opposite side to the consid-
ered node and making a simple proportion. Another faster alternative is to employ the definition
of scalar and vector product as explained in detail in the implementation of the algorithm.
3.2: Graphical representation of the internal force for a node in a membrane - [MaMo211]
By imposing the resultant in the three directions equal to zero for equilibrium, the coordinates
of the considered point are obtained. The iterative process consists of two cycles (one “internal”
and one “external”) rather than the classical procedure seen in many other methods. The cycle
called ”internal” works on the single node keeping unchanged the positions of the other nodes
and updating each time only the one considered until the resultant of the forces in the three
directions is null (or it is reasonable to consider it null). The “external” cycle instead operates
only when, at the end of the inner one, the i-th node is balanced so passing to the next and
repeating the internal one for each node. Once all the nodes are scrolled, the iteration restarts
checking again if the residual of the forces at the node are greater than the fixed threshold
(convergence tolerance). The convergence criterion is based on displacements, checking when
the difference between the position in an iteration and the one in the previous is lower than the
desired tolerance .
The formulation just presented is the most intuitive one, but it can be modified and written more
compactly witht he following:
In which 𝑈t+1 and 𝑈t are vectors of the nodal coordinates at the iteration t+1 and t, A is the
matrix of coefficients, a symmetric matrix that depends only on the coordinates of the other
two nodes belonging to the same element of the analysed node. Finally, B is a vector which
contains the constant terms. They are expressed as follows:
⎧ ∑︀
⎪
⎪ 𝐴 xx = 𝑞i (𝑥bi − 𝑥ai )2
⎪
⎪ ∑︀
⎪
⎪
⎪ 𝐴yy = 𝑞i (𝑦bi − 𝑥ai )2
⎪
⎨𝐴 = ∑︀
zz 𝑞i (𝑧bi − 𝑥ai )2
∑︀ (3.10)
⎪
⎪ 𝐴xy = 𝑞 i (𝑥bi − 𝑥ai ) (𝑦bi − 𝑦ai )
⎪
⎪ ∑︀
⎪
⎪ 𝐴xz = 𝑞i (𝑥bi − 𝑥ai ) (𝑧bi − 𝑧ai )
⎪
⎪ ∑︀
⎩
𝐴yz = 𝑞i (𝑦bi − 𝑦ai ) (𝑧bi − 𝑧ai )
⎧ ∑︀
⎪
⎨𝐵x = ∑︀ 𝑞i [(𝑦bi − 𝑦ai ) (𝑥bi 𝑦ai − 𝑥ai 𝑦bi ) + (𝑧ai − 𝑧bi ) (𝑥bi 𝑧ai − 𝑥ai 𝑧bi )]
𝐵y = 𝑞i [(𝑥bi − 𝑥ai ) (𝑦bi 𝑥ai − 𝑦ai 𝑥bi ) + (𝑧ai − 𝑧bi ) (𝑦bi 𝑧ai − 𝑦ai 𝑧bi )] (3.11)
⎪
⎩ ∑︀
𝐵z = 𝑞i [(𝑥bi − 𝑥ai ) (𝑧bi 𝑥ai − 𝑧ai 𝑥bi ) + (𝑦ai − 𝑦bi ) (𝑧bi 𝑦ai − 𝑥ai 𝑦bi )]
It can be shown that for infinite iterations the terms of matrix A vanish and the difference
between 𝑈t+1 and 𝑈t is reduced, so the process converges [MaMo211]. As the iterations grow,
a limit can be identified in terms of the solution described by:
𝑈t→∞ = (𝐼𝑑 − 𝐴)−1 𝐵 (3.12)
(With Id the identity matrix 3x3)
This last formulation of the problem allows to significantly streamline the process in terms of
steps to be made and consequently in terms of memory and computational cost required. The
previous “internal cycle” has been removed due to the possibility of using the aforementioned
equation to directly calculate the position to which the node tends without the need for itera-
tions. In any case, the cycles executed remain two and are one that runs through the nodes in
the same iteration and one that causes the computation to start again at the first node once the
other node is completed.
The explanation of the orthogonality of the forces can be performed in two very similar ways.
The first consists in computing the equivalent nodal loads of the stress state and checks that
the equilibrium of each element is satisfied. The other method employs the principle of virtual
works. However, in a first phase both will follow a completely identical procedure.
First of all, the domain taken in consideration must be defined. Without any loss of generality,
a generic triangle has been chosen as example of the domain, with different side lengths and
angles:
The numbering followed was to call the sides as the opposite node; a similar reasoning has
been made for the lengths of the sides and angles: the side of length c is opposite to the angle
𝛾. The first step performed is the determination of the shape functions of the nodes: a quick
rule is to divide the implicit equation of the straight line passing through the two other nodes
of the element than the one for which the shape function is being calculated, for the distance
between this straight line and the node considered in the calculation. For example, the shape
function for node 1 can be calculated with the following:
{︃ (︁ )︁
b sin γ
−𝑦 − c cos β 𝑥 + sinb γ = 0
(3.13)
ℎ = sinb γ
To highlight the forces transmitted to the nodes, it is necessary to calculate the equivalent nodal
loads associated to the distributed one. The calculation will be shown for a single node (the
others are analogous). Taking for example node 3, these are the equivalent loads derived from
side 1 and from side 2:
∫︁ c
𝐹13 = 𝑁3 𝜎0 sin 𝛽 d𝑦 (3.16)
0
where 𝐹13 is the equivalent nodal force on the 1st node due to the load applied along side 3.
Therefore:
∫︁ c
𝑦
𝐹13 = 𝜎0 sin 𝛽 d𝑦
0 𝑏 sin 𝛾
∫︁ c
𝜎0
= 𝑦 d𝑦
0 𝑐 (3.17)
[︁ 𝜎 ]︁c
0 2
= 𝑦
2𝑐 0
𝑐
= 𝜎0
2
This force has the same direction of the distributed load and with an analogous procedure it is
possible to determine also the equivalent nodal force on the 2nd node due to the load applied
along side 3:
∫︁ b
𝑦
𝐹23 = 𝜎0 sin 𝛾 d𝑦
0 𝑏 sin 𝛾
∫︁ b
𝜎0
= 𝑦 d𝑦
0 𝑐 (3.18)
[︁ 𝜎 ]︁b
0 2
= 𝑦
2𝑏 0
𝑏
= 𝜎0
2
The other forces are:
• node 1:
𝑏
𝐹21 = 𝜎0
2 (3.19)
𝑎
𝐹31 = 𝜎0
2
• node 2:
𝑏
𝐹32 = 𝜎0
2 (3.20)
𝑐
𝐹12 = 𝜎0
2
3.5: Graphical representation of the nodal forces for the SSDM - [Bia302]
3.6: Enlargement of a node for the calculation of the nodal forces - [Bia302]
Summing the parallel and perpendicular components (with respect to side 3) of each force:
⎧
⎪
⎪ 𝐹13 sin 𝛽 − 𝐹23 sin 𝛾 =
⎪
⎪
⎪
⎪ = (︁2c 𝜎0 sin 𝛽 − 2b 𝜎0 )︁
sin 𝛾 =
⎪
⎪
⎪
⎪ sin γ
⎪ c b
⎪= 2 𝜎0 − 2 𝜎0 sin β sin 𝛽 =
⎪
⎪
⎨
=0
(3.21)
⎪
⎪
⎪
⎪
⎪
⎪
⎪
⎪ 𝐹13 cos 𝛽 + 𝐹23 cos 𝛾 =
⎪
⎪
⎪= c 𝜎0 cos 𝛽 + b 𝜎0 cos 𝛾 =
⎪
⎪
⎪ 2 2
⎩= a 𝜎
2 0
This result is in accordance with what has been stated previously. For the remaining nodes the
same procedure is applied: the forces in the orthogonal and parallel directions to the opposite
side are added to check, as in this case, that the two parallel components vanish whereas the two
orthogonal ones are added together and equal to half of the resulting distributed load applied
on the opposite side.
In the setup of the problem we will refer to the first phase followed in the previous sub-section:
the studied domain is the same as reported in figure 3.3, therefore, the same shape functions
as in the equations (3.14) and (3.15) are employed. From the latter the displacement field
represented in figure 3.7 is constructed:
⎡ ⎤
𝑈1
⎢ ⎥
[︂ ]︂ [︂ ]︂ ⎢𝑈2 ⎥
𝑢(𝑥, 𝑦) 𝑁1 𝑁2 𝑁3 0 0 0 ⎢ ⎥
⎢ 𝑈3 ⎥
= (3.22)
𝑣(𝑥, 𝑦) 0 0 0 𝑁1 𝑁 2 ⎢
𝑁3 ⎢ 𝑉 1 ⎥⎥
⎣ 𝑉2 ⎦
𝑉3
According to the PVW theory, the equilubrium is guaranteed by equating the work of the
external forces with the one related to the internal ones. In purely theoretical form we have:
∫︁ ∫︁
(︀ T )︀ (︀ T )︀
𝑃 𝛿𝑢 dΓ = 𝜎 𝛿𝜀 dΩ (3.24)
Γ Ω
where P represents the field of forces applied on the edge Γ of the element and with Ω the area
of the domain. Since we want to find the forces concentrated in the nodes (and verify their
direction) the first member of the equation is developed as:
∫︁
(︀ T )︀
𝑃 𝛿𝑢 dΓ = 𝐹1T 𝛿𝑢1 + 𝐹2T 𝛿𝑢2 + 𝐹3T 𝛿𝑢3 (3.25)
Γ
in which 𝐹1 , 𝐹2 and 𝐹3 are the vectors of the forces applied respectively in nodes 1, 2 and 3;
while 𝛿𝑢1 , 𝛿𝑢2 and 𝛿𝑢3 are the virtual displacement vectors of the same nodes. At this point it
is enough to displace one degree of freedom fixing the others, in order to find the components
of the forces at the node corresponding to the activated displacement. The derivatives of the
shape functions are explicitly calculated with respect to the two variables:
⎧ A1
⎪
⎪ 𝑁1,x = − B
⎪
⎪
1
⎪
⎪ 𝑁 1,y = − 1
⎪
⎪ B1
⎨𝑁 = 1
2,x a
1
(3.26)
⎪
⎪ 𝑁 2,y = −
⎪
⎪ aA 2
⎪
⎪𝑁3,x = 0
⎪
⎪
⎩
𝑁3,y = A13
In which:
⎧
sin γ
⎪
⎪
⎪ 𝐴1 = cb cos β
⎪
⎨𝐴
2 = tan 𝛾
(3.27)
⎪𝐴3
⎪ = 𝑏 sin 𝛾
⎪
⎪
⎩𝐵 b
= sinγ
1
An example will be proposed for one node (node 3) but the procedure is similar for the others.
Following the path previously described, the horizontal displacement of node 3 is activated to
compute the horizontal component of the force applied. The correspondent deformation vector
is:
[︀ 1
]︀
𝑈3 = 1 ⇒ 𝛿𝜀 = 0 0 A3 (3.28)
As far as the vertical components are concerned, the vertical displacement is:
[︀ 1
]︀
𝑉3 = 1 ⇒ 𝛿𝜀 = 0 A3
0 (3.30)
The dynamic relaxation method is based on concepts which are fundamentally different from
the other two (FDM and SSDM). The DRM is in fact a step by step method that determines
the equilibrium position starting from an initial configuration through a pseudo-dynamic anal-
ysis. It must be underlined that the analysis performed is not for dynamic purposes, but purely
static. As a consequence, boundaries conditions and loads are applied on the structure to ob-
tain the static equilibrium. A fictitious displacement is assigned but it is completely irrelevant
for the achievement of the result. The dynamic component of the whole process is found in
the determination of the mass properties and fictitious damping assigned to all the points of
the discretization. These properties are arbitrary but a proper definition would lead to a faster
execution of of the program and a less intensive convergence towards the static solution.
Historically the DR was developed in 1965 by A.S. Day for the analysis of finite differences in
pressurized concrete tanks, although the concept was also known in earlier times by Rayleigh.
From 1970 onwards M.R. Barnes deepened the study of this method and implemented it for the
analysis and design of pure tensile structures. The method turns out to be particurarly efficient
for this structural typology due to the fact that every iteration always refers to the coordinates of
the various points of the mesh in the position calculated in the previous iteration; this allows the
method to automatically take into account large displacements (geometric linearity). Another
advantage is that it does not require the direct creation of an assembled stiffness matrix of the
structure, avoiding the formulation of a large system of nonlinear equations (non-linearity of
the material or geometry): the equations are solved point by point [Bia302].
𝐹 = 𝑀 𝑎 = 𝑀 𝑣˙ (4.1)
31
Form-Finding and Mechanical Behaviour of Cable and Membrane Structures
This equation would be written for every point of the discretization, assigning eventually dif-
ferent masses at different points if required to obtain a faster convergence. Conceptually, these
fictitious masses are associated to a fictitious damping, in particular a viscous damping, which
is proportional to the speed. The formulation becomes:
All the terms with the apex t are quantities computed at a determined time t. This equation can
be particularized in the three directions having different masses and damping coefficients in the
same point. Starting from an initial instant, the time is increased of a discretized quantity ∆𝑡,
small enough to improve the approximation of the solution. During the calculation velocities
are assumed linear on each interval ∆𝑡 and for each iteration the speed considered is the average
one, calculated starting from those in 𝑡 + ∆𝑡/2 and 𝑡 − ∆𝑡/2. Therefore, we can write:
1 (︁ t+ ∆𝑡 t− ∆𝑡
)︁
m
𝑣ix = 𝑣ix 2 + 𝑣ix 2 (4.3)
2
The accelerations are assumed constant:
t+ ∆𝑡 t− ∆𝑡
t 𝑣ix 2
− 𝑣ix 2
(4.4)
𝑣˙ ix =
∆𝑡
Therefore, the equation (4.2) becomes:
𝑀ix (︁ t+ ∆𝑡 t− ∆𝑡
)︁ 𝐶 (︁ ∆𝑡
ix t+ t− ∆𝑡
)︁
𝑅it = 𝑣ix 2 − 𝑣ix 2 + 𝑣ix 2 + 𝑣ix 2 (4.5)
∆𝑡 2
Taking advantage of the velocities computed in the time interval, the residual force on the node
has been computed and it is referred to the instant 𝑡 + ∆𝑡. From this equation it is possible to
compute also the value of the velocity for the next iteration:
(︃ )︃ (︃ )︃
M𝑖𝑥 C𝑖𝑥
t+ ∆𝑡 t− ∆𝑡
∆t
− 2 1
𝑣ix 2 = 𝑣ix 2 M𝑖𝑥 C𝑖𝑥
+ 𝑅it M𝑖𝑥 C𝑖𝑥 (4.6)
∆t
+ 2 ∆t
+ 2
The variation of the coordinate is then reconstructed and, simply adding this difference to the
previous position, the new coordinate becomes:
t+ ∆𝑡
∆𝑥it+∆t = ∆𝑡𝑣ix 2 (4.7)
t+ ∆𝑡
𝑥t+∆t
i = 𝑥ti + ∆𝑡𝑣ix 2 (4.8)
If properly assigned, an important part is played by the initial conditions (even if they are irrel-
evant for the equilibrium): fixing initial velocities different from zero is possible but not con-
venient; in fact, this would change the loads value on the membrane: loads are not the applied
ones but are modified by the introduction of the damping force. This represents undoubtedly a
useless complication in a practical sense.
Applying a null initial velocity:
∆𝑡
− ∆𝑡
𝑣ix2 = −𝑣ix 2 (4.9)
Then, introducing this equation in the (4.6) for the time ∆𝑡/2, the velocity at time ∆𝑡/2, that
is required in order to have a null initial velocity for time 𝑡 = 0 (as function of the mass and of
the external load) can be computed:
(︂ )︂ (︂ )︂
∆𝑡 𝑀ix 𝐶ix − ∆𝑡 𝑀ix 𝐶ix t
∆𝑡 ∆𝑡 𝐹ix
𝑣ix
2
+ = 𝑣ix 2
− + 𝑃ix 𝑣ix2 = (4.10)
∆𝑡 2 ∆𝑡 2 2𝑀ix
where 𝑃ix
t
represents the external load vector applied in the x-direction.
All the passages must be repeated for the other directions.
The success of the method is based on the kinetic energy of the structure remaining at the end
of the iterations: when the kinetic energy is lower than a certain fixed threshold (so that it is
possible to consider it null), the structure has reached the equilibrium condition and the iterative
process is concluded. The kinetic energy is calculated as:
{︃ ∑︀ ∑︀
𝑇k = ni m 2
j 𝑀ij 𝑣ij
(4.11)
𝑇k ≤ 𝑇0
The geometrical properties are defined a priori and are chosen based on the design requirements
and not for the purposes of the method; the characteristics that can be arbitrarily chosen during
this step are the time interval ∆𝑡, the nodal masses and the damping coefficient. The DR is gen-
erally satisfactory (it converges without encountering numerical instabilities) when the interval
∆𝑡 is not too large or when the masses are not too small; for this reason, by reducing the ∆𝑡
(i.e. by increasing the sampling points) and enlarging the masses (thus providing greater inertia
to the structure), it is possible to make the system more stable. This approach is consired “con-
ditionally stable”, whereas Barnes demonstrated that, by employing a time-dependent mass
value, the problem can become “unconditionally stable”:
∆𝑡2
𝑀ix = 𝐾ix (4.12)
2
where 𝑀ix is the x-directed stiffness in the node 𝑖th . This stiffness is calculated differently
depending on whether we are considering a network of cables or a membrane (or a combination
of the two).
• For cables or beams the elastic stiffness is found as:
𝐸𝐴i 𝑇i
𝐾i = + (4.13)
𝐿i 𝐿i
where 𝑇i is the tension in the 𝑖th element and 𝐿i its current length.
• Instead, in the case of a membrane, considering the texture discretized as triangular ele-
ments with constant thickness, it results:
(︂ )︂
3 𝐿2
𝐾iy = 𝐸𝑡𝑏 𝐿1 + (4.14)
2 𝐿 1 𝐿2
where E is Young modulus of the material, t is the thickness, b is the mean distance between
the nodes in the x-direction, 𝐿1 and 𝐿2 are the distances of the adjacent nodes in the y-direction
from the analysed node.
These process is executed for all directions (x,y,z) and the maximum of the three masses is
assumed for all the nodes, since for the form finding it is advantageous to have the largest
possible mass.
For what concerns the damping factor, it mitigates the vibration effects of the membrane, guid-
ing it towards the static equilibrium. So, based on its value, the convergence may present
different trends.
4.1: Displacement diagram for a 1 dof system for the DRM - Edited from [Wiki411]
Figure 4.1 displays the trends of the displacemenst according to different valuse of the damping
factor (𝜁 = 𝑐/𝑐cr ). The value which ensures the fastest convergence towards the solution is the
so called: “critical damping factor”, 𝜁cr = 1. For this value, the displacement pass the solution
and then tends to it asymptotically, after reversing its motion. For overdamped (𝜁 > 𝜁cr ) or
underdamped (𝜁 < 𝜁cr ) membranes the convergence is slower and therefore not optimized
from a computational point of view.
The value of the critical damping coefficient 𝑐cr is:
√︀
𝑐cr = 2 𝐾i 𝑀i (4.16)
The kinematic damping can be used as an alternative to its viscous counterpart; the kinematic
damping is in addition often priviliged by the designer because it requires the definition of less
parameters. As already shown in (4.16), the viscous damping requires to know the various
nodal masses, the stiffnesses and the integration interval; the kinematic damping requires only
the masses and the time interval ∆𝑡; in fact, in this case, the factor which guides the conver-
gence is the kinetic energy. Cundall developed this method studying mechanically unstable
rocks, finding out that the DRM is extremely stable and fast as far as the resolution of large
displacement problems are concerned.
In practical applications, no damping coefficient is considered but the kinetic energy trend in
the case of an undamped system is drawn until the first maximum: at time 𝑡 + ∆𝑡 the kinetic
energy is less than the previous instant 𝑡, so that the velocity is fixed at 0 and a new iteration will
be performed. This procedure is motivated by the behaviour of 1 degree of freedom systems,
where the peak of the kinetic energy corresponds to the static equilibrium condition (for a
multi-degree of freedom system, it is necessary to let the velocity vanish in correspondence of
more than one peak until convergence is reached).
If the position of the peak coincides with the time 𝑡+∆𝑡/2, that is the instant in which velocities
are computed, the coordinates of the actual current position will be computed in the same
instant (remember that velocities are estimated at time 𝑡 + 𝑘∆𝑡/2, whereas the coordinates at
𝑡 + 𝑘∆𝑡). It can be demonstrated that carrying on with the iterations using this new position
may possibly lead to the divergence of the solution. Therefore, the coordinates computed at
time 𝑡 − ∆𝑡/2 are set as reference coordinates (they are related to the last velocity computed
before the peak). So, it results:
t+ ∆𝑡
𝑥t+∆t
i = 𝑥ti + ∆𝑡 𝑣ix 2 (4.18)
from which:
t+ ∆𝑡
𝑥ti = 𝑥it+∆t − ∆𝑡 𝑣ix 2 (4.19)
t− ∆𝑡 ∆𝑡 t− ∆𝑡
𝑥i 2
= 𝑥ti − 𝑣 2 (4.20)
2 ix
Replacing (4.19) in (4.20):
(︁ )︁ ∆𝑡
t− ∆𝑡 t+ ∆𝑡 t− ∆𝑡
𝑥i 2
= 𝑥it+∆t − ∆𝑡 𝑣ix 2 − 𝑣ix 2 (4.21)
2
Finally, replacing now (4.21) in (4.6) (for the case with 𝑐ix = 0), it results:
(︂ t )︂
t− ∆𝑡 3 (︁ t+ ∆𝑡 )︁ 1 𝑅ix
𝑥i 2
= 𝑥it+∆t
− 𝑣 2
∆𝑡 + ∆𝑡2 (4.22)
2 ix 2 𝑀ix
Now, when the analysis is restarted, the velocities must be computed in the middle point of the
previous interval:
(︂ t )︂
t+ ∆𝑡 1 𝑅ix
𝑥i 2
= ∆𝑡 (4.23)
2 𝑀ix
4.3: Kinetic energy in the neighbourhood of the peak for the DRM - [Hen304]
The basic aim of the finite element method is to solve partial differential equations (PDEs),
or systems of coupled PDEs, numerically. Instead of finding the analytical solution of the
PDE, which is usually a function of the coordinates and therefore not trivial to develop, the
equation is solved numerically for a set of coordinates of grid points. For this purpose, the
continuum is discretized with a number of so-called finite elements which partition the body in
a mesh, which consists of both the defined grid nodes and elements. If the finite elements are
appropriately small, the solution of the PDE can be approximated with a simple function, the
so-called shape function, in each element, which acts as a contribution to the approximation
of the global solution of the PDE. Finally, a linear or non-linear system of equations must be
obtained. with the advantage that (non-)linear systems of equations can now be easily solved
by efficient computer algorithms [Hol405].
In order to develop a thorough model of the problem at hand, some fundamental steps are
strictly required; each of them involves the inclusion of numerical errors found in the final
result.
The first two passages are:
1. Modelling:
This phase is present in all engineering studies: it translates the physical system
to a mathematical model, which abstracts some aspects of interest, focusing atten-
tion on a few aggregate variables and “filtering” the remaining ones. For example,
molecular interactions are not taken into account when calculating the bending
moment of a beam. The physical system, if complex, can be subdivided into sub-
systems, for example a ship or an airplane are composed by a lot of components,
37
Form-Finding and Mechanical Behaviour of Cable and Membrane Structures
where each of them can be considered as subsystem. The subsystem will then be
divided into finite elements that “carry” a mathematical model. It is necessary that
the chosen mathematical model is adequate to the simple geometries of the finite
elements. Choosing an element type in a software program is equivalent to an im-
plicit choice of the underlying mathematical model. The error that can lead the use
of a model must be evaluated with experimental tests, an operation that is generally
expensive in terms of time and resources.
2. Discretization:
In a numerical simulation it is necessary to pass from an infinite number of degrees
of freedom (which is a prerogative condition of the “continuum”) to a finite num-
ber (nodal variables). Discretization, in space or time, has the aim of obtaining a
discrete model characterized by a finite number of degrees of freedom. The error
is given by the discrepancy with the exact solution of the mathematical model. The
most important parameter of the discretization is the characteristic length, which
prescribes the desired size of the elements in the mesh. This size field can be uni-
form, specified by values associated with points in the geometry, or defined by
general “fields”. The definition of the characteristic length depends on the element
geometry and formulation and it is a typical length of a line across an element for a
first-order element; it is half of the same typical length for a second-order element.
For beams and trusses (and cables) it is the length of the element axis. For mem-
branes and shells it is a characteristic length in the reference surface. Elements that
have aspect ratios close to one are recommended, in fact elements with large aspect
ratios will have rather different behavior depending on the direction in which loads
are applied and so the solution is said to be mesh sensitive because of this aspect
[Aba001].
We shall now give the reader a brief overlook of the basics of the finite element formulation for
cable and membrane elements, which will be later extensively employed in Calculix and Dlubal
for the solution of tensile structures, both in the pre-processing phase of the form-finding and
for the application of the additional live loads.
The essential characteristic of the cables is that they are flexible structures which can undergo
great displacements, so the mechanical analysis carried out is highly nonlinear. From the me-
chanical point of view, a cable is able to sustain axial traction and acquires stiffness with respect
to compression if it was adequately pre-stressed (otherwise compression forces would result in
failure).
5.2: Cable element: two nodes which can be subject only to axial force - [Nas007]
The first aspect to be defined is the axial force acting in this particular finite element. The stress
condition is due to two parameters:
• the pre-stress assigned as a thermal load;
𝑁prestress = (𝛼𝐸𝐴) ∆𝑇 with ∆𝑇 > 0 (5.1)
• the elongation caused by the external loads.
So, the axial force is:
𝑁 = 𝐸𝐴 [𝜂 − 𝛼 (𝑇 − 𝑇0 )] (5.2)
with:
{︃
𝐸 > 0, if (𝑇 − 𝑇0 ) < 0
(5.3)
𝐸 = 0, if (𝑇 − 𝑇0 ) ≥ 0
where 𝐸t represents the tangent apparent elastic modulus or effective elastic modulus and is
equal to:
𝐸mat
𝐸t = 1 γ 2 l2 (5.4)
1+ 12 σ 3
𝐸mat
with 𝑙 the projection of the chord length on the horizontal plane. The latter expression is valid
for cables with supports either on the same level or inclined and is demonstrated by putting
small increases in stresses in relation with small increases in strains.
From the expression of the tangent apparent elastic modulus, we understand that 𝐸t is smaller
than 𝐸mat and the following limit holds:
𝐸mat
lim 𝐸t = lim 1 γ 2 l2
= 𝐸mat (5.5)
σ→+∞ σ→+∞ 1+ 12 σ 3
𝐸mat
5.3: Graphical representation of the tangent and secant apparent modulus of elasticity. -
[Mal117]
In case of finite increments, starting from a reference configuration, we can define the secant
apparent modulus as follows:
𝐸mat
𝐸s = (5.6)
γ 2 l2 (1+µ)4
1+ 3
12σ𝑚
· 16µ2
𝐸mat
The equivalent stress 𝜎 is a function 𝜇 = 𝜎i /𝜎s (respectively initial and final stresses) and the
mean stress 𝜎m = 𝜎s (1 + 𝜇)/2:
[︂ ]︂ 13
16𝜇2
𝜎 = 𝜎m · (5.7)
(1 + 𝜇)4
In a linearized structural analysis, 𝜎 can be used to compute the elastic tangent modulus, instead
of using the secant modulus.
In order to find the displacement field, u for a static problem, the equilibrium position is sought
for a cable structure which is initially in a generic configuration and subject to a system of
forces; in order to solve this problem, it is necessary to employ an iterative non linear analysis,
as thoroughly described in 5.3.
To obtain the equations characterizing the behavior of the membranes, it is necessary to intro-
duce several hypotheses. The goal here is not to prove how these hypotheses have been defined,
but they will be given directly without going into details; a deeper insight can be gained refer-
ring to [LeV111].
The main hypothesis underlying the modellization are [DeS003b]:
• The membrane displacements coincide with the displacements of the mean surface;
• The membrane thickness 𝑡 is very small compared to the characteristic length of the
membrane;
• The stresses are constant in the thickness.
• We consider the membrane to be, locally, in plane stress situation, as reported in 5.4:
⎡ ⎤
𝜎xx 𝜎xy 0
σ = ⎣𝜎yx 𝜎yy 0⎦ (5.8)
0 0 0
An important aspect not yet analysed is how to determine the modulus of elasticity of a struc-
tural element of this kind. As widely described so far, the weaving process represents a signif-
icant source of non-linearity which should not be neglected, furthermore, the composite nature
of membranes impose to deepen the mechanical response. An interesting reference is repre-
sented by the ASCE STANDARD on Tensile Membrane Structures ([ASCE100]) in which a
procedure is illustrated to calibrate the elastic constants for an orthotropic elastic law.
The calibration is done on the basis of experimental results obtained on a cruciform membrane
specimen subjected to biaxial loading, for different load ratios in the two main direction (the
fibers direction).
During the test, the stresses are imposed in both directions keeping fixed the ratio between warp
and fill stress. In figure 5.6 and 5.7 are represented the results of the test, the curves labelled
as fill 1:2 represents stress-strain relations concerning fill stress and strain obtained by loading
according to the ratio 1:2 that is imposising a load in the fill direction which is the double of
the one in warp direction.
Both material test results display a nonlinear nature, mainly on the first loading application.
The two materials show different behaviour (mainly in the first loading): the fiberglass shows
a strain hardening, whereas the polyester a strain softening. The mechanical response tends to
more linear under successive load cycles.
Note that the stress range for which the stress-stain curve are plotted is limited with respect to
the bracking strength of the materials; in fact, for design purposes the stresses are kept lower
than one-forth of the ultimate stresses.
Due to the non linearity displayed, if the modulus of elasticity were to be found directly from
these curves, it would vary with the stress. This result is not practical to be implemented, it
increases the computational cost.
Generally, a constant value is extracted selecting a slope for each stress ratio. The Least Squares
Method is usually applied to determine the values of the moduli which embodied the best fitting
with respect to the linearized experimental curves.
The same procedure could be adopted also for the estimate of the elastic moduli representing
the material behaviour under successive load cycles.
In order to solve systems of nonlinear equations there are many methods which are mainly
based on the Newton procedure. We must also bear in mind that the the problem formulation
influences the solution algorithm to be adopted. For this reason, it is important to properly
understand and implement the most appropriate method.
With reference to a discretized structural problem, for which the nodal displacements (collected
in the vector u) are the primary unknowns, one of the most widely used algorithm to solve
nonlinear static problems is the Newton-Raphson Method [Kim114]. It is an iterative method
that allows to find u using the equilibrium condition defined as:
where Fext are the given nodal loads and Fint (u) are the internal forces developed on the nodes
by the structural elements (obviously for a linear problem 𝐹int = 𝑘𝑢 with 𝐾 being the overall
stiffness matrix of the structure).
Numerical methods for solving a system of nonlinear equations are based on an initial estimate,
u0 , with which it is possible to find a new increment, ∆u, so that the new estimate, u0 + ∆u,
is closer to the solution to the equilibrium equation (5.9). The increments are determined ap-
proximating locally the nonlinear equations by means of linearization. This process is repeated
until the original nonlinear equations are satisfied. Suppose an approximate solution at the 𝑖th
iteration is known and is equal to ui . The solution at the next iteration can be approximated
using the first-order Taylor series expansion as follows:
where
(︂ )︂
𝜕Fint
Ji (ui ) ≡ (5.11)
𝜕u i
that is the Jacobian matrix at the 𝑖th iteration, commonly known as the tangent stiffness matrix
in structural applications and ∆ui is the solution increment.
This passage is aimed to calculating ∆ui and iteratively update the solution, ∆ui+1 . Once all
the terms have been rearranged, the system of linearized equations can be obtained as:
Equation (5.12) is similar to the matrix equation of linear systems, except that:
1. The coefficient matrix, Ji (ui ), is not constant, but a function of ui ;
2. The equation is solved for the increment, ∆ui , not for the final solution, u;
3. The second member side is not the applied force, but the difference between the applied
force and the internal force (called the residual).
Once the system of equations has been solved for the displacement increment, ui , a new ap-
proximate solution is obtained as follows:
Until the correct solution is not found, this result will not exactly fulfill the system of nonlinear
equations and there will be some residual or unbalanced force which is computed as in:
The final solution can be considered attained when the residual force is smaller than a given
tolerance, the result, ui+1 , is accepted as the accurate solution, and the process exits. Otherwise,
the iterations continue until this residual becomes smaller than the tolerance. The convergence
parameter is expressed in the normalized form as follows:
∑︀n 2
j=1 (𝑅i+1,j )
𝑐𝑜𝑛𝑣 = ∑︀n 2 (5.15)
j=1 (𝐹ext,j )
The iterations are terminated when the convergence parameter, conv, becomes smaller than a
given tolerance (i.e. 0.001).
Another criteria, based on the solution increment, can be applied to determine the convergence
of the iterative procedure. When the increment of the solution is much smaller than the initial
increment, then it is assumed that the solution has converged. Thus, the convergence criterion
can also be written as:
∑︀n 2
j=1 (∆𝑢i+1,j )
𝑐𝑜𝑛𝑣 = ∑︀n 2 (5.16)
j=1 (∆𝑢1,j )
5.9: Two iterations of the Newton-Raphson procedure for a system with a single DOF - Edited
form [Kim114]
For a single DOF system, the Jacobian matrix is the slope of the nonlinear function, Fint . The
solution converges faster if the initial configuration is close to the solution, for this reason
in practice is common to apply gradually the desired load with a set of increments. When
the iteration analysed is close to the solution, the Newton-Raphson method shows a quadratic
convergence.
Assuming 𝑢exact as the exact solution, and 𝑢n and 𝑢n+1 be the two consecutive approximations
of the solution. Then, the method converges quadratically when there exists a constant 𝑐 > 0
such that:
where:
• |𝑢exact − 𝑢n+1 | is the error at the (𝑛 + 1)th iteration
• |𝑢exact − 𝑢n |2 is the square of the error at the 𝑛th iteration
In practice, the exact solution is the objective of the problem and therefore is not known ex
ante, so the the converged solution is assumed as 𝑢exact .
To demonstrate that this numerical algorithm converges quadratically, it must be proved that
the following ratio tends to a constant value 𝑐:
|𝑢exact − 𝑢n+1 |
lim =𝑐 (5.18)
n→∞ |𝑢exact − 𝑢n |2
In practice, it is often enough to prove that the convergence criterion in. (5.12) reduces quadrat-
ically at each iteration.
It is worth noticing that the Newton–Raphson method does not always guarantee convergence
to the accurate solution [Kim114]:
• This method assumes that the solution increment in (5.12) is relatively small. Therefore,
after every iterations the value of ∆𝑢 becomes smaller and may also tend to zero as the
solution accuracy improves. However, when the Jacobian matrix becomes singular, or
the determinant of matrix J is null, ∆𝑢 becomes infinite and the solution diverges. This
means, in a single DOF system, that the slope of Fint becomes zero and the residual
cannot be reduced. Numerically, a similar behavior can be observed when the matrix is
nearly singular.
• If the starting point is too far away from the exact solution, the method may diverge or
oscillate between two points. This also happens when the curvature of the Fint curve
changes its sign between two consecutive iterations. In such a case, it is possible that the
Newton–Raphson algorithm may result in an infinite loop. In order to prevent an infinite
loop, a maximum number of iterations is imposed and the algorithm exits with an error
message when the number of iterations reaches the maximum allowed number.
The Jacobian matrix of structural mechanics problems is normally positive definite. Physically,
this means that in order to increase the displacement, the applied force should increase. How-
ever, in some cases, the displacement may increase without the applied force increasing. Or,
equivalently, the displacement continuously increases while the applied force decreases. This
causes structural instability.
As well explained in [Kim114], a common examples of this problem are bifurcation and snap-
through behaviors. Figure 5.10 shows the snap-through behavior of an elastic inclined slender
beams with both ends clamped.
The plot shows the relation between the applied force and the vertical displacement at the force
application point. Initially, the applied force increases along with the tip displacement (region
AC). In this region, the beams are stable, and the Jacobian matrix is positive definite. When they
reach point C, the Jacobian matrix becomes singular, and the force cannot increase beyond 𝐹C .
Between points C and E, the tip displacement continuously increases while the force decreases.
The system is unstable in this region. Beyond E, the beams become stable again with a positive
definite Jacobian matrix. When the Newton–Raphson method is used, it can only solve for the
response in the region, AC. For example, for a given force, 𝐹B , it always yields point B, not
point D. Special techniques are required to follow the force–displacement curve.
As final remark, it is relevant to underline that starting from this method, other two strategies,
the Modified Newton-Raphson Method and the Incremental Secant Method have been devel-
oped. The first one is characterized by a less expensive procedure from a computational point
of view, but it needs more iterations to find the appropriate solution. On the other hand, the
Incremental Secant Method provides the same benefits in term of computational efficiency, but
it guarantees a greater level of convergence rate.
5.11: Modified Newton-Raphson Method (left) and Incremental Secant Method (right) - Edited
from [Kim114]
The Newton–Krylov method is a numerical method for solving non-linear problems using
Krylov subspace linear solvers.
In this extract we are going to underline the main aspects of the method, more details are
developed in [KnKe210] and [XuCo217].
As already explained in paragraph 5.3.1, the core of the Newton method is to successively
solve linear systems, (5.19), in which the n-by-n Jacobian matrix is computed at every iteration
[XuCo217].
Fint (ui ) + Ji (ui ) ∆ui = Fext (5.19)
This procedure could be computationally expensive if the problem is characterized by a huge
number of unknowns, thus, approximations of the Jacobian matrix are often used to replace
the exact Jacobian matrix in the Newton iteration (Modified Newton-Raphson Method or In-
cremental Secant Method, figure 5.11).
Although the computational cost of each iteration is thereby reduced, more iterations are usu-
ally required to converge to the root of (5.19).
However, the Jacobian matrix itself is structured in many applications, so it may not be neces-
sary to form the Jacobian matrix itself explicitly in each iteration. In 2004, Knoll and Keyes
([KnKe210]) reviewed the Jacobian-free Newton-Krylov method (JFNKM) for solving the non-
linear equations (5.19); this approach combines a matrix-free approach with nonsymmetric
Krylov subspace methods. This method is a projection, or generalized projection, method for
solving a linear system 𝐴𝑥 = 𝑏 using the Krylov subspace 𝐾j [XuCo217]:
𝐾j = 𝑠𝑝𝑎𝑛(𝑟0 , 𝐴𝑟0 , 𝐴2 𝑟0 , . . . , 𝐴j−1 𝑟0 ) (5.20)
where 𝑟0 = 𝑏 − 𝐴𝑥0 .
One of the most widely used method is the Generalized Minimal RESidual method (GMRES),
which is an Arnoldi-based method (further explanation of the Arnoldi method can be found in
[LeSo109] and [LSY113]). In the GMRES, the Arnoldi basis vectors form the trial subspace
out of which the solution is constructed. A matrix-vector product is required to create a new trial
vector per each iteration of the method. The method terminates on the basis of an estimation
of the residual. For this reason, the GMRES method ensures that is not necessary to form the
parameter matrix A explicitly. The method only requires the matrix-vector product (i.e. σx)
[KnKe210].
As properly shematized in [XuCo217], in order to construct a matrix-vector product 𝐽(𝑥)𝑢 for
(5.19), the finite difference method can be implemented. For a nonlinear function, its corre-
sponding Jacobian-vector product can be approximated as:
𝐹int (𝑥 + 𝜖𝑢) − 𝐹 (𝑥)
𝐽(𝑥)𝑢 = (5.21)
𝜖
where 𝜖 is a small positive number.
Thus, combining the GMRES method with finite differences, the JFNKM method is obtained
for solving nonlinear equations. In practice, this method works well in solving problems in
fluid dynamics, power systems, plasma physics and so on.
As a final remark, it is known that (5.21) is a first-order approximation of the Jacobian vector
product. We can increase its accuracy to second order with:
𝐹int (𝑥 + 𝜖𝑢) − 𝐹 (𝑥 − 𝜖𝑢)
𝐽(𝑥) ≈ (5.22)
𝜖
but this technique doubles the cost of forming 𝐽(𝑥)𝑢.
We consider a general deformable continuum as represented in figure 5.12, that in its deformed
configuration is boundend by a surface 𝑆 and has a volume 𝑉 . More precisely, we will dis-
cuss about the Cauchy Continuum, from the scholar to whom are traditionally associated the
definitions of the measurement of stress and strain [CoTa105].
The external actions acting on the body, supposed known, are:
• Volume force, F: force per unit of volume as, self-weight, inertia forces, etc.;
• Surface force, f: force per unit of surface, of known value on the unbounded surface, 𝑆f ;
• Imposed displacement, s0 , of the constrained surface, 𝑆u .
All the actions acting on the body are generally dependent on the point in which they are
applied. We will indicate with x the vector containing the generic point position. Therefore,
the generic notation for functions of position will be F(x). The components of the vectorial
quantity just introduced will be referred to a system of orthogonal Cartesian axes, and all the
components will be grouped in column vector:
⎡ ⎤ ⎡ ⎤ ⎡ ⎤
𝐹x 𝑓x 𝑠0x
F = 𝐹y f = 𝑓y s0 = 𝑠0y ⎦
⎣ ⎦ ⎣ ⎦ ⎣ (5.23)
𝐹z 𝑓z 𝑠0z
We assume the external actions vary in time in a quasi-static manner, so that the dynamic effect
can be neglected. Moreover, all the displacement which the body undergoes are of small entity,
so that it is possible to assume the small strains and displacements hypotesis, which allows to
consider the deformed configuration equal to the undeformed one.
The unknown fields characterizing the stress-strain behaviour of the body are:
• Displacement, s(x):
⎡ ⎤
𝑠x
s = ⎣ 𝑠y ⎦ (5.24)
𝑠z
• Deformation, ε(x):
⎡⎤
𝜀xx
⎢ 𝜀yy ⎥
⎢ ⎥
⎢ 𝜀zz ⎥
ε=⎢
⎢𝛾xy ⎥
⎥ (5.25)
⎢ ⎥
⎣𝛾yz ⎦
𝛾zx
• Stress, σ(x):
⎡⎤
𝜎xx
⎢𝜎yy ⎥
⎢ ⎥
⎢ 𝜎zz ⎥
σ=⎢
⎢𝜎xy ⎥
⎥ (5.26)
⎢ ⎥
⎣ 𝜎yz ⎦
𝜎zx
In order to find these fields we employ different equations:
• Indefinite Equilibrium Equation:
These are 3 first order, linear, differential equations which in tensorial notation can be written
as:
𝜎ij,i + 𝐹j = 0 (5.27)
𝜎ij 𝑛i = 𝑓j (5.28)
• Kinematic Compatibility:
These are 6 first order, linear, differential equations, which in tensorial notation can be written
as:
1
𝜀ij = (𝑠i,j + 𝑠j,i ) (5.29)
2
with the associated boundary condition, which represents the displacement imposed or allowed
by the constrains (= 0 in case they are fixed):
𝑠j = 𝑠0j (5.30)
• Constitutive Law:
These are 6 linear algebric equations, that link together stress and strain by means of the elastic
tensor:
{︃
𝜎ij = 𝑑ijhk 𝜀hk Direct Form
(5.31)
𝜀ij = 𝑐ijhk 𝜎hk Inverse Form
Where 𝑑 represent the stiffness of the system and 𝑐 the flexibility. This tensor is symmetric and
its independent components (elastic constants) are at most 21:
σ = dε (5.34)
Orthotropic materials are a subset of anisotropic materials which, in contrast to the most gen-
eral case, have two preferred directions of strength which are mutually perpendicular. The
properties along these directions (also known as principal directions) are the extreme values of
the elastic coefficients [SoWo008].
The formulation follows the same passages explained in paragraph 5.4.1; the constitutive law
is:
⎡ ⎤ ⎡ E𝑥 ν𝑥𝑦 E𝑦
⎤⎡ ⎤
𝜎xx 1−ν𝑥𝑦 ν𝑦𝑥 1−ν𝑥𝑦 ν𝑦𝑥
0 𝜀xx
⎣𝜎yy ⎦ = ⎢ ν𝑦𝑥 E𝑥 E𝑦 ⎥⎣ ⎦
⎣ 1−ν𝑥𝑦 ν𝑦𝑥 1−ν𝑥𝑦 ν𝑦𝑥 0 ⎦ 𝜀yy (5.35)
𝜏xy 0 0 𝐺 𝛾 xy
Note that in order to satisfy the symmetry in the matrix, the following relationship linking the
Poisson’s ration and the elastic modulus in the two principal directions must be imposed:
The elastoplastic behaviour of material occurs when a structural element is subject to both elas-
tic and plastic deformation. The hypothesis of small strains and displacements is maintained
such that geometric nonlinear effects can be neglected and the focus is centred on the material
nonlinearity.
In order to consider the elastoplastic behaviour of a given structure, strains are assumed com-
posed by two quantities:
These governing equations hold ∀𝑡, but in this contest the time is a pseudo-time, with often the
meaning of load multiplier. In fact, plasticity is a time independent phenomenon, the response
does not depend on the velocity with which we load the material.
Before proceeding with the discussion, it is important to define some quantities:
1. Plastic Multiplier, 𝜆(𝑡). It is a parameter which by definition can only grow and depends
on time. So its characteristics can be resumed as follow:
⎧
∫︁ t
⎪
⎨𝜆(𝑡 = 0) = 0
𝜆= ˙
𝜆𝑑𝑡 ⇒ 𝜆(𝐸) = 0 (5.39)
0 ⎪
⎩
d𝜆 ≥ 0
2. Plastic Modulus, 𝐻 which is the slope of the strain-hardening portion of the stress–strain
curve after removal of elastic strain components.
⎧
∆𝜎 ⎨𝐻 > 0 Hardening
⎪
𝐻= → 𝐻 = 0 Perfect Plasticity (5.40)
∆𝜀p ⎪
⎩
𝐻 < 0 Softening
3. Linear Hardening, ℎ(𝜆):
ℎ(𝜆) = 𝐻𝜆 (5.41)
4. Yield Function, 𝜙 that is a function which describes the state of the structure at a given
time:
{︃
𝜙 < 0 material in elastic domain
𝜙 = |𝜎| − (𝜎0 + 𝐻𝜆) ≤ 0 ⇒ (5.42)
𝜙 = 0 material on the border
Thanks to this definitions it is possile to compute the yield stress:
and the consistency condition, a useful relation to describe the loading/unloading paths:
{︃
𝜙≤0
⇒ 𝜙𝜆˙ = 0 ∀𝑡 (5.44)
˙𝜆 ≥ 0
d𝜙 d|𝜎| d𝜎
=
𝜙˙ = − 𝐻 𝜆˙
d𝑡 d𝜎 d𝑡
= sign(𝜎) [𝐸 (𝜀˙ − 𝜀˙p )] − 𝐻 𝜆˙ (5.46)
The plastic strain rate (𝜀˙p ) can be different from zero only for 𝜙 = 0 and 𝜙˙ = 0, due to the
consistency condition (5.44)
Therefore, it results:
(︂ )︂
𝐸
𝜆˙ = sign(𝜎)𝜀˙ (5.48)
𝐸+𝐻
𝜎˙ = 𝐸 𝜀˙e
= 𝐸 (𝜀˙ − 𝜀˙p )
(︂ )︂
𝐸
= 𝐸 𝜀˙ − 𝐸 𝜀˙
𝐸+𝐻
(︂ )︂
𝐸2 (5.50)
= 𝐸− 𝜀˙
𝐸+𝐻
𝐸𝐻
⇒ 𝐸T =
𝐸+𝐻
Starting from the differential equation (5.38), substituting the value of 𝐸T , eq. (5.50), we obtain
the elastoplastic constitutive law for one-dimensional systems:
(︂ )︂ (︂ )︂
𝐸𝐻 𝐸𝐻
d𝜎 = d𝜀 ⇒ 𝜎 = 𝜀 + const
𝐸+𝐻 𝐸+𝐻 (5.51)
⇒ 𝜎 = 𝐸T (𝜀 − 𝜀0 ) + 𝜎0
The basic concepts presented in the previous section for one-dimensional systems can be gen-
eralized for multidimensional systems. Obtaining the the required stress–strain relation for
one-dimensional system is relatively straightforward, in fact, experiments are most commonly
uniaxial. At the same time, one-dimensional strain-hardening models are also relatively easy to
develop because either the magnitude of yields stress increases (isotropic hardening) or main-
tains the same range between tension and compression (kinematic hardening). However, one-
dimensional elastoplasticity has a restricted range of practical applications, i.e. bars and trusses.
In the case of two or three dimensional structural systems, the elastoplastic theory is way more
complex to be implemented: the stress becomes a tensor with up to six components instead of
a scalar quantity. Therefore, to test the material with a procedure similar to the one used for the
previous case, different combinations of stress components must be applied. Of course, this is
is practically impossible because the possible combinations are almost infinite.
To better understand the difficulties at the base of this problem, let us consider biaxial loading
of a plane stress structure. The only nonzero stress components are 𝜎xx and 𝜎yy . A possible
procedure to set different stress combination is to fix 𝜎xx at 500 MPa (for example) and gradu-
ally increase 𝜎yy beyond the yielding limit of the material. At this point, the procedure must be
repeated fixing 𝜎xx at different values. To obtain the stress–strain relationship for biaxial load-
ings, it is necessary to test all these possible combinations. Therefore, it is barely impossible to
determine the constitutive law for multidimensional systems except for very limited cases.
To overcome the problem and define the multidimensional elastoplasticity stress-strain relation-
ship a physic-based model able to represents all possible cases is implemented. Furthermore, it
is assumed a scalar measure of stress and strain as representative of the multidimensional stress
status.
Given that the multidimensional model must be consistent also for the one-dimensional case,
the stress–strain relationship can be obtained from uniaxial tests. The main aspect is to properly
define the form of yield criterion with a hardening rule and the elastoplastic constitutive law
relating incremental stress to strain in the plastic region [Kim114].
In 3D, the elastoplastic model is based on:
⎧
⎪
⎪ 𝜀ij = 𝜀eij + 𝜀pij Elastic + Plastic II Order Symmetric Tensor
⎪
⎪
⎪ e
Constitutive Law(𝜈, 𝐸)
⎨𝜎ij = 𝑑ijhk 𝜀hk
⎪
𝜙(𝜎ij , 𝜆) = ‖𝜎ij ‖ − (𝜎0 + ℎ(𝜆)) Scalar Yield Function (5.53)
⎪
⎪
⎪
⎪ 𝜙 ≤ 0; d𝜆 ≥ 0; 𝜙 d𝜆 = 0 Complementary Equations
⎪
⎪
⎩d𝜀p = 𝑛 d𝜆 = ∂ϕ d𝜆 Assumed Associated Flow Rule
ij ij ∂σ𝑖𝑗
Since it is possible to choose different shapes for 𝜙, we assume the Von Mises cylinder, which
is a proper choice for isotropic materials, since the invariant stresses are used fo the domain:
√︂
3 (5.54)
𝜙(𝜎ij , 𝜆) = 𝑆ij 𝑆ij − (𝜎0 + ℎ(𝜆))
2
√︁
where 32 𝑆ij 𝑆ij is the Von Mises equivalent stress (or norm), and 𝑆ij is the deviatoric tensor:
⎧ (︀ )︀
⎪
⎪ 𝑆11 = 32 𝜎11 − 31 𝜎22 − 31 𝜎33
⎪
⎪ (︀ 1 )︀
⎪
⎪ 𝑆 22 = − 𝜎 11 + 2
𝜎 22 − 1
𝜎 33
⎪
⎪ (︀ 3 3 3 )︀
1 ⎨ 𝑆33 = − 13 𝜎11 − 31 𝜎22 + 23 𝜎33
𝑆ij = 𝜎ij − 𝜎kk 𝛿ij ⇒ (5.55)
3 ⎪
⎪ 𝑆 12 = 𝜎12
⎪
⎪
⎪
⎪ 𝑆13 = 𝜎13
⎪
⎪
⎩
𝑆23 = 𝜎23
Having fixed ℎ(𝜆) = 𝐻 𝜆, our model represents linear isotropic hardening, with the Von Mises
yield criterion orientation and with an associated flow rule.
Since the flow rule depends on the choice of 𝜙, we need to show the consequence of this chioce.
We underline the diagonal and extradiagonal term:
⎛ ⎞
𝜕𝜙 3 𝑆12 ⎠
d𝜀p12 = d𝜆 = ⎝ √︁ d𝜆
𝜕𝜎12 2 3
𝑆 𝑆
2 ij ij
⎛ ⎞ (5.56)
𝜕𝜙 3 𝑆11 ⎠
d𝜀p11 = d𝜆 = ⎝ √︁ d𝜆
𝜕𝜎11 2 3
𝑆 𝑆 2 ij ij
Thus:
⎛ ⎞ √︂
3 𝑆ij 2 p p
d𝜀pij = ⎝ √︁ ⎠ d𝜆 = d𝜀ij d𝜀ij = d𝜀pequiv (5.57)
2 3
𝑆 𝑆 3
2 ij ij
given that d𝜀p11 + d𝜀p22 + d𝜀p33 = 0, because of the cylinder is isotropic and plasticity can
not provoke any volume variation. As a consequence, it is possible to define the Von Mises
plasticity 3D model, by means of the following quantities:
• Stress and deviatoric stress tensor;
• Strain and deviatoric strain tensor;
• Elastic strain energy density (the material is considered isotropic linear elastic)
• Principal stress space.
First of all, we can determine:
• Hydrostatic Stress (mean stress):
𝜎11 + 𝜎22 + 𝜎33 𝜎kk
𝜎m = = (5.58)
3 3
• Hydrostatic Stress Tensor:
𝜎kk
𝑝ij = 𝛿ij (5.59)
3
• Deviatoric Stress Tensor:
𝜎kk
𝑆ij = 𝜎ij − 𝛿ij (5.60)
3
So, the stress tensor can be written as:
𝜎kk
𝜎ij = 𝑆ij + 𝛿ij (5.61)
3
The same decomposition can be applied to the strain tensor obtaining:
𝜀kk
𝜀ij = 𝑒ij + 𝛿ij (5.62)
3
This result is the sum of the deviatoric strain tensor and the volumetric/spherica/mean strain
tensor. So, the dilatation is:
∆𝑉
𝜀kk = (5.63)
𝑉0
Due to the previous definitions, under the hypothesis of small strains and displacements and
isotropic linear elastic material, the elastic strain energy density (strain energy per unit of vol-
ume) results:
1 1 (︁ 𝜎kk )︁ (︁ 𝜀kk )︁
Ω = 𝜎ij 𝜀ij = 𝑆ij + 𝛿ij 𝑒ij + 𝛿ij
2 2 3 (︂ 3 )︂
1 1
= 𝑆ij 𝑒ij + 𝜎kk 𝜀kk
2 3
[︂ (︂ )︂ ]︂ (5.64)
1 𝑆ij 1 1 (︁ 𝜎kk )︁
= 𝑆ij + 𝜎kk
2 2𝐺 3 𝐾 3
1 11 1
= 𝑆ij 𝑆ij + 𝜎hh 𝜎kk
4𝐺 29𝐾
where 𝐾 is the Bulk Modulus, defined ad : 𝐾 = 𝐸/3(1 − 2𝜈)
The distortion energy theory is also called the Von Mises Yield Criterion, which states that
yielding occurs when the equivalent stress reaches the yield stress of the material in uniaxial
tension [Kim114].
Before proceeding with the discussion the invariants of the tensors must be defined:
⎧
⎪
⎨𝐼1 = 𝜎kk
𝐼2 = 21 𝜎ij 𝜎ji (5.65)
⎪
⎩
𝐼3 = 13 𝜎ij 𝜎jh 𝜎ki
⎧
⎪
⎨𝐽1 = 𝑆11 + 𝑆22 + 𝑆33
𝐽2 = 12 𝑆ij 𝑆ji (5.66)
⎪
⎩
𝐽3 = 13 𝑆ij 𝑆jh 𝑆ki
It is possible to write the Von Mises equivalent stress, the norm shown in equation (5.54), and
the definitions provided in (5.66), as:
√︀
‖𝜎‖ = 3𝐽2 (5.67)
Using the equivalent stress, the yield function and the corresponding yield criterion can be
defined as:
Since any isotropic scalar function of a symmetric tensor can be described as a function of the
principal values of its argument, it follows that any iso-surface (i.e. any subset of the function
domain with fixed function value) of such functions can be graphically represented as a surface
in the space of principal values of the argument. This allows, in particular, the yield surface of
any isotropic yield criterion to be represented in a particularly simple and useful format as a
three-dimensional surface in the space of principal stresses [DePeOw106].
5.15: The Tresca and Von Mises yield surfaces in principal stress space - [DePeOw106]
5.16: The 𝜋-plane in principal stress space (a), and, the 𝜋-plane representation of the Tresca
and von Mises yield surfaces (b) - [DePeOw106]
Cutting pattern
As described in chapter 1.3, the last step before the static analysis is the patterning of the mem-
brane resulting from the form finding. Structural fabrics are manufactured as flat panels with a
typical width of 2–3 m, and a maximum width of 5 m. Tensile structure made of membranes
are generally implemented in buildings characterized by larger dimension. In order to fulfill
the production requirements it is necessary to “split” the membrane in patterns and these can-
not be flattened in a plane without some distortion being in general the surface undevelopable
(Gaussian curvature not null). The extracted pieces from the cutting pattern procedure will be
sewed together and consequently tensile fabric structures cannot be produced without inducing
stresses in the surface which appear even before the application of the live loads. If poorly
designed, the final shape of the panels may affect the final form and stress distribution of the
whole structure.
In terms of structural properties, the most efficient shape for the membrane is a combination
of double curvature. This form allows to reduce the displacements due to the vertical loads but
requires a special design procedure. Patterning seeks to determine the arrangement of the planar
59
Form-Finding and Mechanical Behaviour of Cable and Membrane Structures
fabric panels such that, when the panels are assembled, the desired 3D shape is achieved, and
the stress distribution is as close as possible to the one imposed during form finding [GaLe206].
Historically patterning was conducted using physical models. S. Gale and W.J. Lewis in 2015
describe how computational methods of patterning could be implemented for this aim. They
underlined five steps in which the general process of patterning could be divided:
1. Seam definition, concerns the division of the form-found membrane surface into panels.
The divisions are defined by seams which generally take the form of sewn and/or welded
overlapping panels of fabric. It is considered good practice to establish seam lines along
geodesics. Geodesics are lines of minimum distance over a surface, and are the path
adopted by a constant stress cable stretched over the surface. Consequently, geodesic
seams do not introduce undesirable stresses into the membrane. Seams physically dictate
the panel size, and affect the subsequent flattening, stress reduction and compensation of
the individual panels. For practical reasons, seams should run through the regions of low
curvature to avoid possible wrinkling during sewing and welding. At the same time, they
should be spaced reasonably, to ensure the curvatures across the panels remain low.
2. Flattening concerns the development of a portion of the 3D form-found membrane sur-
face into the 2D plane.
3. Stress reduction concerns the application of iterative methods to the flattened panel
geometry to reduce stresses.
4. Compensation concerns the shrinking of the pattern to account for tensioning of the
membrane during construction.
5. Pattern assembly can be included in patterning schemes to calculate the final geometry
and stresses in the constructed membrane. The cutting pattern is assembled according to
the physical boundary conditions, and relaxed into its equilibrium shape, giving the final
geometry and stresses.
6.1 Examples
Dlubal RFEM5 provides the user with a module that performs the cutting pattern, where the
flattening process is performed according to the minimum energy theory.
The main advantage of Dlubal RFEM5 is that, for each pattern, the compensation can be ap-
plied in both the warp and weft direction, giving the user also the possibility to set a special
compensation value for each boundary line as well as overlaps for manufacturing processes.
6.1. Examples 61
Form-Finding and Mechanical Behaviour of Cable and Membrane Structures
6.1. Examples 63
Form-Finding and Mechanical Behaviour of Cable and Membrane Structures
Before diving into the numerical formulation of the form-finding procedure of mixed cable and
membrane structures, the natural behaviour of cables subjected to own-weight and distributed
load are inspected analytically.
An analytical approach based on the solution of ordinary differential equations in closed form
is presented to determine the shape developed by a cable hinged at its edges and subjected to its
own weight. On the one hand, the method is limited to very simple cases where both the bound-
ary conditions and the statical loads allow to find a closed form solution of the general integral
of the ODE, but on the other hand the analytical formulation allows to capture straightforwardly
the interaction between the various geometric and static parameters of the problem.
This approach could also be extended to cases for which an analytical solution is generally not
trivial, but the effort for the identification of the internal boundary conditions in the presence
of heterogeneity or discontinuities is not justified in comparison to the numerical approach of
the force density method which results in a much easier development of the solving algorithm.
7.1.1 Introduction
Let us imagine a taut cable with an assigned horizontal reaction 𝐻 due to the pretensioning
of the cable itself. Our objective is the identification of the displaced shape of the cable when
subject to its own weight and the pretensioning. In order to investigate the relationship between
the pretensioning 𝐻 and the dead load 𝑔, which is related to both the mechanical (density 𝜌)
and geometric (the area 𝐴) properties of the cable, the solution of the equation of the catenary
cable is repeated multiple times varying first the applied prestress having fixed the external
load, and subsequently viceversa.
65
Form-Finding and Mechanical Behaviour of Cable and Membrane Structures
7.1.3 Solution
The solution of the equation of the catenary is performed multiple times to generate the plots
shown after the discussion of the alogorithms. Every solution is delegated to the following
function which automatically performs the spatial discretization of the cable and then employs
the analytical formulas to retrieve the axial force and the horizontal and vertical support reac-
tions to be passed to the main function for the plotting.
solve_catenary.solve_catenary(ax, ay, bx, by, g, tension, space=1000)
Solves the form finding for the catenary equation based on the given applied horizontal
force.
The tension is represented by the horizontal force applied to both ends of the taut cable;
if no additional horizontal loads are present in the form-finding process, the latter is kept
constant.
Parameters
• ax (float) – x-coordinate of the left pin of the cable.
• ay (float) – y-coordinate of the left pin of the cable.
• bx (float) – x-coordinate of the right pin of the cable.
• by (float) – y-coordinate of the right pin of the cable.
• g (float) – distributed load equivalent to the specific weight per
unit length of the cable.
• tension (float) – horizontal force to be prescribed in the form-
finding process.
• space (float) – discretization of the cable to compute the dis-
placements and internal actions.
Returns dictionary grouping the retvals of the problem.
Return type
dictionary containing the return values with the following standard and
convention:
• x (ndarray) - numpy array containing the location of the nodes span-
ning along the length of the cable representing the numerical dis-
cretization of the geometrical entity.
• y (ndarray) - numpy array containing the displacements of the nodes
under the applied gravitational load.
• t (ndarray) - numpy array containing the tension of the cable func-
tion of the x-coordinate spanning along the length of the cable.
respected.
• ta (ndarray) - tension of the cable at the left pin.
• tb (ndarray) - tension of the cable at the right pin.
• va (ndarray) - vertical reaction of the support at the left pin.
• vb (ndarray) - vertical reaction of the support at the right pin.
• g (ndarray) - load for which the problem has been solved with the
current run.
• tension (ndarray) - horizontal external force for which the form-
finding has been performed.
The script is invoked with the main function, which is responsible for the the generation of the
plots for the comparison between two significant static conditions:
1. A prescribed horizontal force is assigned and the dead load insiting on the cable, stem-
ming from the section of the considered element, is varied.
2. Keeping the external load constant, the horizontal component of the pretension is varied
incrementally.
catenary.main()
Main entry point for the form-finding of the catenary cable. No special interface is pro-
vided due to simplicity of the study; the user must supply the values directly in the
command file.
The target of the form-finding procedure is the horizontal force in the cable. The script
investigates the problem in two different situations:
1. The horizontal force is kept constant and the gravitational load is varied.
2. The gravitational load is kept constant and a different horizontal force for the
form-finding is prescribed.
The exhibited behaviour is highly nonlinear as pointed out by the figures produced by
matplotlib.
Returns None
Return type None
Code block 7.2: Main entry for the program solving the
form-finding for a catenary cable.
1 from solve_catenary import solve_catenary
2
3 import matplotlib.pyplot as plt
4
5
6 def main():
7 # Coordinates of the left support.
8 ax = 0.0
9 ay = 0.0
10 # Coordinates of the right support
11 bx = 100.0
12 by = 2.0
13 # List containing the applied gravitational loads.
14 g = [20.0, 25.0, 30.0, 35.0, 40.0]
15 # Vector containing the horizontal forces to be applied
16 # symmetrically at the supports.
17 tension = [25.0e3, 30.0e3, 35.0e3, 40.0e3, 45.0e3]
18 # Discretization of the cable.
19 spacing = 10000
20
21 plt.figure(
22 num="Displacements depending on the applied load for a catenary")
23 for grav in g:
24 retvals = solve_catenary(ax, ay, bx, by, grav, tension[-1],
25 spacing)
26 plt.plot(retvals["x"], -retvals["y"],
27 label="g = %.0f N/m" % retvals["g"])
28 plt.title("Deformed shape for a catenary")
29 plt.xlabel("Abscissa [m]")
30 plt.ylabel("Displacement [m] for H = %.0f N" % tension[-1])
31 plt.legend()
32 plt.grid()
33
34 plt.figure(
35 num="Displacements depending on the horizontal applied force for "
36 "a catenary")
37 for load in tension:
38 retvals = solve_catenary(ax, ay, bx, by, g[0], load, spacing)
39 plt.plot(retvals["x"], -retvals["y"],
(continues on next page)
7.1.4 Post-processing
The post-processing of the results is delegated to the python package matplotlib, which provides
the user with not only a programmimg interface similar to Matlab, but also a more refined
interaction based on object classes.
7.1: Displacement shape obtained varying the external load for a catenary cable.
7.1.5 Tests
The program is tested with pytest, a python module which allows to write better programs by
systematically testing the code against significant situations to avoid, in particular for numerical
algorithms, annoying regressions.
test_catenary.test_catenary()
Tests the form-finding algorithm for the catenary cable in the case of a taut cable, whose
supports are located at (0,0)m and (100, 2)m. The cable is subjected to a gravitational
load of 40 N/m computed as:
𝑔 = 𝐴 · 𝛾 = 𝐴 · 𝑑𝑔 (7.8)
7.2: Displacement shape obtained varying the pretension for a catenary cable.
Code block 7.3: Test case for the program solving the form-
finding for a catenary cable.
1 from solve_catenary import solve_catenary
2 import pytest
3
4
5 def test_catenary():
6 # Coordinates of the left pin support.
7 ax = 0.0
8 ay = 0.0
9 # Coordinates of the right pin support.
10 bx = 100.0
11 by = -2.0
12 # Applied gravitational load.
13 g = 40.0
14 # Horizontal applied force.
15 tension = 17157.2875
16
17 # Solve the form-finding problem.
18 retvals = solve_catenary(ax, ay, bx, by, g, tension)
19 # Extract the tensions and vertical reactions at the ends for the
20 # comparison in the test.
21 tensions_test = {k: retvals[k] for k in ('ta', 'tb', 'va', 'vb')}
(continues on next page)
7.1.6 Conclusions
The first conclusion that can be drawn from the previous examples is that the cable does not
have a configuration of its own, since, without pretensioning, the theoretical displacement at-
tained by the element is unlimited:
𝐻 (︁ 𝑔 · 𝑥 )︁
lim 𝑦(𝑥) = lim − · cosh − − 𝐶 + 𝐷 −→ ∞ (7.9)
H→0 H→0 𝑔 𝐻
Even if with the case of a distributed load, one can understand that the cables assumes the
shape of the funicular curve of the insisting forces; when the pretension 𝐻 increases, the rays
of the polygon flatten, causing the increase of the axial action in the cable. On the other hand,
the cables becomes extremely slack for low values of pretensioning. As discussed later in the
thesis, this behaviour is endowed with serious drawbacks from a numerical point of view, since
the extremely large rotations of a slack cable may lead to the divergence of the nonlinear solver
if the latter is employed with the average parameters for the number of iterations and maximum
number of line searches.
Finally, we can observe that the behaviour of the cable is linear for small loads or high values
of the pretensioning, whereas the load-displacement relationship becomes severly nonlinear
for large loads or especially for moderate values of the pretensioning that lead to the aforemen-
tioned slack deformation of the element.
7.2.1 Introduction
The example of the parabola cable under a distributed load illustrates the concept of “inverse
form-finding”: having assigned a certain shape to the cable, we shall find the actions that
need to be introduced in the system in order to guarantee the static equilibrium, being the
kinematic compatibility already assured since it has been imposed as the starting hyphotesis of
the problem.
The final shape of the cable is assigned by means of the three quantities 𝑓a , 𝑓b and ℎ with the
meaning as reported in 7.3.
The rotational equilibria with respect to the left and right pins are:
𝑝𝑎2 𝑝𝑏2
𝐻𝑓a − =0 𝐻𝑓b − =0 (7.10)
2 2
Subtracting the first equation from the second one:
𝑝 (𝑏2 − 𝑎2 )
𝐻 (𝑓b − 𝑓a ) − =0 (7.11)
2
and considering that 𝑎 + 𝑏 = 𝑙, then:
𝑙 𝐻ℎ 𝑙 𝐻ℎ
𝑎= − 𝑏= + (7.12)
2 𝑝𝑙 2 𝑝𝑙
Putting the new definition of 𝑎 in the previous equation and solving the second degree equation
we obtain:
(︂ )︂
𝑝𝑙2 𝑓a + 𝑓b √︀
𝐻= 2 ± 𝑓a 𝑓b (7.13)
ℎ 2
then the vertex is comprised between A and B; otherwise, the solution is found with the plus
and being the cable very taut, the vertex is outside the segment AB:
(︂ )︂
𝑝𝑙2 𝑓a + 𝑓b √︀
𝐻= 2 + 𝑓a 𝑓b (7.15)
ℎ 2
𝑝𝑙2 (︀ )︀
𝑦(𝑥) = · 𝜉 − 𝜉 2 + ℎ𝜉 (7.16)
2𝐻
where 𝜉 = xl .
The derivative of the displacement function can be computed as:
(︂ )︂
𝑝𝑙2 1 1 ℎ
′
𝑦 (𝑥) = · − 𝜉 + (7.17)
2𝐻 𝑙 𝑙 𝑙
to finally get:
𝐻 𝐻
𝑇A = 𝑇b = 𝑉A = 𝐻 tan 𝜃A 𝑉B = 𝐻 tan 𝜃B (7.19)
cos 𝜃A cos 𝜃B
7.2.3 Solution
The solution of the parabola cable under a distributed load is conducted following the analytical
formulation explained in the previous paragraph by operating a discretization of the cable to
map the displacement continuosly.
solve_parabola.solve_parabola(l, fa, fb, h, p, spacing=1000)
Solve the inverse form-finding problem for the taut cable under a distributed load, whose
supports are located at different levels and the vertex is at fixed level defined by the
quantities 𝑓a , 𝑓b and ℎ.
Parameters
• l (float) – length of the projection of the cable.
• fa (float) – distance of the vertex from the left support in the y-
direction.
• fb (float) – distance of the vertex from the right support in the
y-direction.
• h (float) – difference in height between the supports.
• p (float) – distributed load acting on the taut cable.
• spacing (int) – discretization of the cable.
Returns dictionary grouping the retvals of the problem.
Return type
dictionary containing the return values with the following standard and
convention:
• x (ndarray) - numpy array containing the location of the nodes span-
ning along the length of the cable representing the numerical dis-
cretization of the geometrical entity.
• y (ndarray) - numpy array containing the displacements of the nodes
under the applied distributed load.
• H (float) - horizontal force at the supports of the cable.
• ta (ndarray) - tension of the cable at the left pin.
• tb (ndarray) - tension of the cable at the right pin.
• va (ndarray) - vertical reaction of the support at the left pin.
• vb (ndarray) - vertical reaction of the support at the right pin.
• p (ndarray) - load for which the problem has been solved with the
current run.
• t (ndarray) - numpy array containing the tension of the cable func-
tion of the x-coordinate spanning along the length of the cable.
Main entry point for the inverse form-finding of the cable under distributed load. No
special interface is provided due to simplicity of the study; the user must supply the
values directly in the command file.
The target of the form-finding procedure is the final configuration of the cable. The script
investigates the problem for different levels of the load, showing that the deformed shape
is kept constants, whereas the tension varies to equilibrate the applied external load.
Returns None
Return type None
Code block 7.5: Main entry for the program solving the ca-
ble under a distributed load.
1 import matplotlib.pyplot as plt
2 from solve_parabola import solve_parabola
3
4
5 def main():
6 # Projected length of the cable.
7 length = 100.0
8 # Difference in height between the two supports.
9 h = 2.0
10 # Vertical distance between the left pin and the vertex of the
11 # parabola.
12 fa = 0.35
13 # Vertical distance between the right pin and the vertex of the
14 # parabola.
15 fb = fa + h
16 # Vector containing the loads to be tested.
17 p = [20.0, 25.0, 30.0, 35.0, 40.0]
18 # Discretazione of the cable.
19 spacing = 10000
20
21 plt.figure(
22 num="Fixed displacement for a cable under distributed load for "
23 "the inverse form-finding problem")
24 for grav in p:
25 retvals = solve_parabola(length, fa, fb, h, grav, spacing)
26 plt.plot(retvals["x"], -retvals["y"],
27 label="p = %.0f N/m - H = %.0f N/m" % (
28 retvals["p"], retvals["H"]))
29
30 plt.title("Deformed shape for a cable under distributed load")
31 plt.xlabel("Abscissa [m]")
32 plt.ylabel("Displacement [m]")
33 plt.legend()
34 plt.grid()
35
36 plt.figure(
37 num="Tension in a cable under distributed load for the inverse "
38 "form-finding problem")
39 for grav in p:
40 retvals = solve_parabola(length, fa, fb, h, grav, spacing)
41 plt.plot(retvals["x"], retvals["t"],
(continues on next page)
7.2.4 Post-processing
The program produces two graphical outputs by using the simple interface provided by mat-
plotlib:
1. Displacement shape of the cable, to show the fact that during the inverse form-finding
process the deformed shape of the cable is actually kept constant.
2. Tension in the cable along the whole span: the tension varies according to the imposed
external distributed load but is rather constant along the span being the vertical offset
between the supports is very limited.
7.2.5 Tests
The tests are invoked using python -m pytest and are executed against the example reported in
the literature in [Mal117].
test_parabola.test_parabola()
Tests the inverse form-finding algorithm for a taut cable under distributed load, whose
supports are located at (0,0)m and (100, 2)m. The cable is subjected to a distributed load
of 40 N/m and its vertex is at fixed level defined by the quantities 𝑓a = 2𝑚 and 𝑓b = 4𝑚.
The supports are distanced 2m in height.
The tests are validated against a case found in the literature by Prof. Malerba from
Politecnico di Milano.
Returns asserts the positiveness of the test.
Return type None
Code block 7.6: Test case for the program solving the cable
under a distributed load.
1 import pytest
2 from solve_parabola import solve_parabola
3
4
5 def test_parabola():
6 p = 40.0
7
8 # Solve the inverse form-finding problem.
9 tensions = solve_parabola(length, fa, fb, h, p)
10
11 # Extract the tensions and vertical reactions at the ends for the
12 # comparison in the test.
13 tensions_test = {k: tensions[k] for k in ('ta', 'tb', 'va', 'vb')}
14 # Asserts the positiveness of the test.
15 assert tensions_test == pytest.approx(
16 {"tb": 17316.548,
17 "ta": 17237.102,
(continues on next page)
7.2.6 Conclusions
The data used for the example of the cable under a distributed load is studied to produce the
same results, in terms of tension, as for the catenary cable. At the same time, the pretension as-
signed to the catenary cable produces a displacement which is almost equal to the one imposed
for the parabola cable.
Therefore, the deformed shape and the statical reactions for the case of the catenary or parabola
cable are very similar, being however the procedure to retrieve the results slightly different:
the catenary cable expects the solution of an ordinary differential equation with hyperbolic
functions, whereas for the parabola cable we can proceed by pure equilibrium.
As logical subsequent step, the non-linear behaviour of cable and membrane structures are
individually analysed with the package CalculiX. Firstly, the same examples analytically de-
veloped in the previous chapter has been performed and checked, later some different load cases
and configurations have been inspected. With the aim to completely understand the possible
non-linear behaviour of taut elements, elasto-plastic constitutive law has been assigned to the
structural elements, therefore the analysis performed presented both material and geometrical
non-linearity. Finally, has been produced an example of a plane membrane bounded with a
collaborating cable which provides the pre-stress.
8.1.1 Introduction
This example has been performed to show the natural behaviour developed by a taut cable
hinged at its ends under its own weight. This effect is more evident when the cable covers
long span; in fact, in order to avoid the sag effect for this particular situation, a high level of
pre-stress should be imposed.
The cable analysed presents a Young modulus equal to 160 · 103 𝑀 𝑃 𝑎 and weight density of
78.5𝑘𝑔/𝑚3 , with a sectional area of 519𝑚𝑚2 and it is pre-stressed with an equivalent thermal
load of 17.14𝑘𝑁 .
8.1.2 Pre-processing
The generation of the mesh is delegated to cgx due to the simplicity of the topology. Also the
solution and post-processing phase are automated directly in Calculix.
83
Form-Finding and Mechanical Behaviour of Cable and Membrane Structures
79 ds -3 e 1
80 plot fv all
81 sys sleep 1
82 hcpy png
83 sys mv hcpy_5.png cable_catenary_stress_xx.png
84
85 ds -2 e 1
86 plot fv all
87 sys sleep 1
88 hcpy png
89 sys mv hcpy_6.png cable_catenary_strain_xx.png
90
8.1.3 Solution
The solution deck is written in ccx for Calculix. We would like to underline again the fact
that the initial prestress is prescribed as an inelastic thermal load computed starting from the
coefficient of thermal expansion.
8.1.4 Post-processing
The characteristic catenary shape is put in evidence by the plots generated automatically by
cgx.
The same exercise is solved in the commercial software Dlubal showing a perfect match with
respect to the results obtained in the open-source software Calculix.
8.1.6 Conclusions
The shape of the catenary cable has been chosen to match the results of the analytical form-
finding performed in the next chapter of the thesis; therefore, the numerical calculations show
the effectiveness of the analytical procedure.
8.1: Total displacements of the catenary cable under its own weight.
1
0.9
Node=51
0.8
0.7
Values at Nodes (cable-catenary.fbd)
0.6
Time
0.5
0.4
0.3
0.2
0.1
0
-200
-300
-400
-500
-600
-700
-800
-900
-1000
-1100
-1200
D2
0.0004
0.00035
0.0003
0.00025
EXX
0.0002
0.00015
0.0001
5x10-5
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
80
70
60
50
SXX
40
30
20
10
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
8.5: Axial stress and elongation of the mid-point of the catenary cable.
2500
2000
1500
F2
1000
500
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
-5000
-10000
-15000
-20000
F1
-25000
-30000
-35000
-40000
-45000
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
8.2.1 Introduction
We now draw a comparison between the catenary and parabola shape for a cable. As discussed
later in the next chapter, the two structural cases show strong similarities both in terms of
displacements and internal actions, but the analytical formulation to obtain the results will
widely diverge: in the case of the catenary cable, an ordinary differential equation must be
solved, whereas for the cable under a distributed load, we can proceed by pure equilibrium.
The cable analysed presents a Young modulus equal to 160 · 103 𝑀 𝑃 𝑎 and weight density of
78.5𝑘𝑔/𝑚3 , with a sectional area of 519𝑚𝑚2 and it is pre-stressed with an equivalent thermal
load of 17.14𝑘𝑁 .
8.2.2 Pre-processing
Being the initial configuration of the cable a straight line, the mesh is directly generated in cgx
without referring to external meshers. At the same time, the solution and the post-processing
phase are launched from Calculix.
4 # point geometry to define the two pins for the cable and also
5 # a node in the middle to test for maximum displacement.
6 pnt p1 0.0 0.0 0.0
7 pnt p2 100e3 -2000.0 0.0
8 pnt p3 50e3 -1000.0 0.0
9
72
73 ds -4 e 3
74 plot fv all
75 sys sleep 1
76 hcpy png
77 sys mv hcpy_4.png cable_parabola_displacement_z.png
78
79 ds -3 e 1
80 plot fv all
81 sys sleep 1
82 hcpy png
83 sys mv hcpy_5.png cable_parabola_stress_xx.png
(continues on next page)
92
93 graph loadv time DISP D2
94 sys mv graph_0.ps cable_parabola_displacement_load.eps
95 graph loadv time STRESS SXX
96 sys mv graph_1.ps cable_parabola_stress_load.eps
97 graph loadv time TOSTRAIN EXX
98 sys mv graph_2.ps cable_parabola_strain_load.eps
99
100 sys ps2pdf cable_parabola_stress_load.eps
101 sys ps2pdf cable_parabola_strain_load.eps
102 sys pdfnup --no-landscape --a4paper --nup 1x2 --outfile cable_parabola_ss_
˓→load.pdf cable_parabola_strain_load.pdf cable_parabola_stress_load.pdf
103
104 graph psx time FORC F2
105 sys mv graph_3.ps cable_parabola_reaction_pin.eps
106 graph psx time FORC F1
107 sys mv graph_4.ps cable_parabola_horizontalf_pin.eps
108
109 sys ps2pdf cable_parabola_reaction_pin.eps
110 sys ps2pdf cable_parabola_horizontalf_pin.eps
111 sys pdfnup --no-landscape --a4paper --nup 1x2 --outfile cable_parabola_
˓→reactions.pdf cable_parabola_reaction_pin.pdf cable_parabola_horizontalf_
˓→pin.pdf
112
113 # let's go back to the deformed shape.
114 rot -z
115 view disp
116 ds -4 e 4
117 plot fv all
118
119 # exit the program
120 # quit
8.2.3 Solution
The distributed load is transformed in an equivalent nodal load either in the deck or in the
pre-processing phase.
8.11: Axial stress and axial elongation of the parabola cable under a distributed load.
1
0.9
Node=51
0.8
0.7
Values at Nodes (cable-parabola.fbd)
0.6
Time
0.5
0.4
0.3
0.2
0.1
0
-200
-300
-400
-500
-600
-700
-800
-900
-1000
-1100
-1200
D2
0.0004
0.00035
0.0003
0.00025
EXX
0.0002
0.00015
0.0001
5x10-5
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
80
70
60
50
SXX
40
30
20
10
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
8.13: Axial stress and elongation of the mid-point of the parabola cable.
2500
2000
1500
F2
1000
500
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
-5000
-10000
-15000
-20000
F1
-25000
-30000
-35000
-40000
-45000
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
8.14: Horizontal and vertical reactions of the the supports of the parabola cable.
8.2.4 Post-processing
The post-processing phase as usual reports the countour plots for the displacements, stresses
and strains.
The same case study is also analysed in the commercial software Dlubal to verify the obtained
results. There is no discrepancy between the two results.
8.2.6 Conclusions
In particular, we focus our attention on the time graphs reporting the maximum dis-
placement for the mid span point: as in the previous examples, independently from
the applied pretension, the behaviour of the cable appears to be highly nonlinear for
the first increments, after which the slope of the graph plotting the displacements
versus the applied load tends to decrease.
We shall now analyse the response of cable structures under the application of live loads. Since
the beginning, we would like to point out that the parameter time reported in the time graphs
represents the load multiplier which amplifies the reference applied load. The amplitude func-
tion assigned by Calculix ranges from 0 to 1 and automatically imposes a load multiplier equal
to the currrent time incrementation.
8.3.1 Introduction
The studied cable structure consists of a single taut cable hinged at both ends, which are sepa-
rated by a span of 8m. The cable is subject to a vertical load applied in the middle of the span
equal to 10𝑘𝑁 . The steel section of the cable has an area of 60𝑚𝑚2 with effective elastic mod-
ulus of 140 · 103 𝑀 𝑃 𝑎. The pretension simulated by means of an equivalent inelastic thermal
load is 20𝑘𝑁 .
8.3.2 Pre-processing
The pre-post and solution steps are run automatically by means of a bash-like script written in
cgx for Calculix.
Code block 8.11: Pre and post-processing file for the cable
under a concentrate load.
1 # each part of the cable is represented by just one finite element.
2 valu ds 1
3
4 # point geometry defining the pins and middle node for the cable.
5 pnt p1 0.0 0.0 0.0
6 pnt p2 8000.0 0.0 0.0
7 pnt p3 4000.0 0.0 0.0
8
74 ds -3 e 1
75 plot fv all
76 sys sleep 1
77 hcpy png
78 sys mv hcpy_5.png cable_concentrated_stress_xx.png
79
80 ds -2 e 1
(continues on next page)
Code block 8.12: Mesh and element sets for the cable under
a concentrated load.
1 *NODE, NSET=Nall
2 1,0.000000000000e+00,0.000000000000e+00,0.000000000000e+00
3 2,4.000000000000e+03,0.000000000000e+00,0.000000000000e+00
4 3,8.000000000000e+03,0.000000000000e+00,0.000000000000e+00
5 *ELEMENT, TYPE=T3D2, ELSET=Eall
6 1, 1, 2
7 2, 2, 3
8.3.3 Solution
The solution deck is written in Calculix, which has a similar syntax to the one adopted by the
commercial solver Abaqus.
8.3.4 Post-processing
The post-processing phase, as already state, is fully automated in cgx, which directly produces
the countour plots for the displacements (magnitude, x and y-direction), the axial stress and the
axial elongation. We also output some interesting time graphs for the same quantities referred
to the mid span versus the load multiplier represented on the abscissa by time.
8.18: Displacements in the x and y-direction of the cable under a concentrated load.
8.19: Stresses and strains for the cable under a concentrated load.
1
0.9
Node=2
0.8
0.7
Values at Nodes (cable-concentrated.fbd)
0.6
Time
0.5
0.4
0.3
0.2
0.1
0
-50
-100
-150
-200
-250
-300
-350
-400
D2
0.004
0.0035
0.003
0.0025
EXX
0.002
0.0015
0.001
0.0005
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
900
800
700
600
500
SXX
400
300
200
100
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
8.21: Axial stress and elongation of the loaded point of the cable.
4500
4000
3500
3000
2500
F2
2000
1500
1000
500
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
-10000
-20000
-30000
F1
-40000
-50000
-60000
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
Performing the same nonlinear geometric analysis in the commercial software Dlubal, the same
results are obtained, even though Dlubal does not activate by default large deformations; how-
ever, being the end rotations of the cable less than 10%, this effect can be considered negligible.
8.3.6 Conclusions
This simple example introduces the reader to the concept of nonlinear elastic analysis, show-
ing the need for an iterative approach to the solution. Being the increase of the load strictly
monotonic, the convergence of the solution is attained in a few steps, considering also the fact,
as outlined in the paragraphs describing the theory behind the method, that no second order
effects, which could lead to instability, are activated.
8.4.1 Introduction
The same example of the cable under a concentrated load has been thoroughly analysed chang-
ing the mechanical properties of the elements; an isotropic hardening law for the cable is as-
signed by means of a table linking the level of stress and the attained Mises equivalente plastic
strain (in our case, the yielding is activated at 800𝑀 𝑃 𝑎 and at 1000𝑀 𝑃 𝑎 the element under-
goes an equivalent plastic deformation of 15%).
The structure is tested against a cycle of loading and unloading to show the permanent defor-
mations in the cable after the discharge; the characteristic sag effect is put in evidence.
This case will not be analysed in Dlubal due to the impossibility to model the cycle without the
necessary additional modules.
8.4.2 Pre-processing
The pre-post processing phase is automated in Calculix with a bash script; the solution by ccx
is also directly invoked by the .fbd script.
44 rot -z
45 view disp
46 scal d 5
47 ds -5 e 4
48 plot fv all
49 sys sleep 1
50 hcpy png
51 sys mv hcpy_1.png cable_plastic_displacement_magnitude_3d.png
52
53 view disp off
54 ds -5 e 1
55 plot fv all
56 sys sleep 1
57 hcpy png
58 sys mv hcpy_2.png cable_plastic_displacement_x.png
59
60 ds -5 e 2
61 plot fv all
62 sys sleep 1
63 hcpy png
64 sys mv hcpy_3.png cable_plastic_displacement_y.png
65
66 sys montage cable_plastic_displacement_x.png cable_plastic_displacement_y.
˓→png -tile 2x1 -geometry +0+0 cable_plastic_displacement_xy.png
67
68 ds -5 e 3
69 plot fv all
70 sys sleep 1
71 hcpy png
72 sys mv hcpy_4.png cable_plastic_displacement_z.png
73
74 ds -4 e 1
75 plot fv all
76 sys sleep 1
77 hcpy png
78 sys mv hcpy_5.png cable_plastic_stress_xx.png
79
80 ds -3 e 1
81 plot fv all
82 sys sleep 1
83 hcpy png
(continues on next page)
113
114 graph vload time PE PE
115 sys mv graph_5.ps cable_plastic_pe_vload.eps
116
117 # reshow the deformed shape
118 rot -z
119 view disp
120 ds -5 e 4
121 plot fv all
8.4.3 Solution
The consitutive law in terms of equivalent plastic deformations is assigned by means of a table
in the deck file; no dependency with regards to the temperature is activated.
8.4.4 Post-processing
The same graphs as in the previous paragraph are here reported; the contour plots are drawn at
the final time after the discharge.
8.4.5 Conclusions
The example clearly shows that the cable should never, under service conditions, yield. In fact,
even though the final sag could be avoided by increasing again the pretension, the permanent
plastic strains have deformed the cable so significantly that the original shape imposed by the
prestressing has been completely lost. The application of additional live loads would only
worsen the situation.
8.27: Axial stresses and strains of the plastic cable under a concentrated load.
8.28: Von-Mises equivalent plastic strain of the cable under a concentrated load.
2.5
Node=2
2
Values at Nodes (cable-plastic.fbd)
1.5
Time
1
0.5
0
0
-50
-100
-150
-200
-250
-300
-350
-400
-450
D2
0.005
0.004
0.003
EXX
0.002
0.001
0
0 0.5 1 1.5 2 2.5
Time
700
600
500
400
SXX
300
200
100
-100
0 0.5 1 1.5 2 2.5
Time
8.30: Axial stress and elongation of the loaded point of a plastic cable.
4500
4000
3500
3000
2500
F2
2000
1500
1000
500
0
0 0.5 1 1.5 2 2.5
Time
-5000
-10000
-15000
-20000
F1
-25000
-30000
-35000
-40000
-45000
-50000
0 0.5 1 1.5 2 2.5
Time
2.5
Node=2
2
Values at Nodes (cable-plastic.fbd)
1.5
Time
1
0.5
0
0.0025
0.002
0.0015
0.001
0.0005
PE
8.5.1 Introduction
We now introduce the example of a simple plane cable net, hinged at every edge.
The cables, with a section of 2𝑚𝑚2 , used for the construction, have an elastic modulus of
110 · 103 𝑀 𝑃 𝑎 and are pretensioned with a load of 500𝑁 . The forces applied to the free nodes
are 50𝑁 .
8.5.2 Pre-processing
The grid points defining the start and end nodes of the cable elements are defined directly in
cgx without the need to employ any external mesher. ccx for the finite element solution and cgx
for the post-processing are automatically invoked by the script.
72 ds -4 e 2
73 plot fv all
74 sys sleep 1
75 hcpy png
76 sys mv hcpy_3.png cable_net_displacement_y.png
77
8.5.3 Solution
We report here the Calculix deck defining the materials, boundary conditions and loads for the
definition of the problem at hand.
8.5.4 Post-processing
The post-processing phase outputs the contour plots for both displacements, stresses and
strains, but the concentration should be focused on the time graphs (i.e. load multiplier) es-
pecially for the displacement.
The same structural case is analysed in the commercial software dlubal to verify the results
obtained with the open-source finite element software Calculix.
8.5.6 Conclusions
The analysis of a structure composed by multiple cables has highlighted the fundamental role
played by the pretensioning: in fact, even though the cables are all prestressed before the ap-
plication of the live loads, a strongly nonlinear relationship between the external load and the
mid displacement is evident. The cables, bearing no flexural stiffness on their own, in order to
collaborate among them efficiently, must be properly pretensioned to allow the structure to re-
spond slightly super-linearly to the external actions. Otherwise, the initial nonlinear behaviour
could worsen as the load increases, causing the typical sag effect of slack cables.
-40
-50
-60
D3
-70
-80
-90
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
-0.1
-0.2
-0.3
D2
-0.4
-0.5
-0.6
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
8.36: Vertical and lateral displacement of the loaded point of the cable net.
8.6.1 Introduction
The analysis of a plane membrane, in addition to showing the general behaviour of this type of
structure, allows the engineer to understand that the membrane can be assimilated to a contin-
uous curtain of cables in two directions.
In this following examples, no pretension is assigned to the membrane for two main reasons:
1. The pretension should be imposed by the action of the cables as later discussed for the
form-finding procedure.
2. We want to highlight the nonlinear behaviour of membranes which are not subject to
prestressing, to show that after an initial severe nonlinear phase, the membrane recovers
substantial stiffness to carry additional loading.
The membrane is characterized by an elastic modulus of 8 · 103 𝑀 𝑃 𝑎 and Poisson’s ratio of 0.4.
It is subject to a distributed load of 1𝑘𝑃 𝑎.
8.6.2 Pre-processing
The mesh geometry is created with the open-source code gmsh, drawing a square with side
equal to 5𝑚 and thickness 2𝑚𝑚.
117 ds -2 e 2
118 plot fv all
119 max 2e-3
120 min -7e-5
121 sys sleep 1
122 hcpy png
123 sys mv hcpy_11.png membrane_plane_strain_yy.png
124
125 sys montage membrane_plane_strain_xx.png membrane_plane_strain_yy.png -
˓→tile 2x1 -geometry +0+0 membrane_plane_strain_xxyy.png
126
127 ds -2 e 7
128 plot fv all
129 max 1.75e-3
130 sys sleep 1
131 hcpy png
132 sys mv hcpy_12.png membrane_plane_strain_m.png
133
(continues on next page)
136
137 # time graphs
138 seta frex n 842
139 graph frex time DISP D3
140 sys mv graph_0.ps membrane_plane_displacement_frex.eps
141 graph frex time STRESS SXX
142 sys mv graph_1.ps membrane_plane_stress_xx_frex.eps
143 graph frex time STRESS SYY
144 sys mv graph_2.ps membrane_plane_stress_yy_frex.eps
145
146 sys ps2pdf membrane_plane_stress_xx_frex.eps
147 sys ps2pdf membrane_plane_stress_yy_frex.eps
148 sys pdfnup --no-landscape --a4paper --nup 1x2 --outfile membrane_plane_
˓→stress_xxyy_frex.pdf membrane_plane_stress_xx_frex.pdf membrane_plane_
˓→stress_yy_frex.pdf
149
150 graph frex time STRESS Mises
151 sys mv graph_3.ps membrane_plane_stress_mises_frex.eps
152 graph frex time TOSTRAIN EXX
153 sys mv graph_4.ps membrane_plane_strain_xx_frex.eps
154 graph frex time TOSTRAIN EYY
155 sys mv graph_5.ps membrane_plane_strain_yy_frex.eps
156
188
189 graph pin time TOSTRAIN Mises
190 sys mv graph_12.ps membrane_plane_strain_mises_pin.eps
191
192 sys ps2pdf membrane_plane_stress_mises_mid.eps
193 sys ps2pdf membrane_plane_strain_mises_mid.eps
194 sys pdfnup --no-landscape --a4paper --nup 1x2 --outfile membrane_plane_se_
˓→mid.pdf membrane_plane_stress_mises_mid.pdf membrane_plane_strain_mises_
˓→mid.pdf
195
196 seta mid n 611
197 graph mid time STRESS SXX
198 sys mv graph_13.ps membrane_plane_stress_xx_mid.eps
199 graph mid time STRESS SYY
200 sys mv graph_14.ps membrane_plane_stress_yy_mid.eps
201
202 sys ps2pdf membrane_plane_stress_xx_mid.eps
203 sys ps2pdf membrane_plane_stress_yy_mid.eps
204 sys pdfnup --no-landscape --a4paper --nup 1x2 --outfile membrane_plane_
˓→stress_xxyy_mid.pdf membrane_plane_stress_xx_mid.pdf membrane_plane_
˓→stress_yy_mid.pdf
205
206 graph mid time STRESS Mises
207 sys mv graph_15.ps membrane_plane_stress_mises_mid.eps
208 graph mid time TOSTRAIN EXX
209 sys mv graph_16.ps membrane_plane_strain_xx_mid.eps
210 graph mid time TOSTRAIN EYY
211 sys mv graph_17.ps membrane_plane_strain_yy_mid.eps
212
8.6.3 Solution
The solution deck written in Calculix imposes, after having defined the material properties, the
boundary conditions of the pinned membrane on four edges.
13 ** Define the thickness for the membrane and assign the material.
14 *SHELL SECTION, ELSET=Ememb, MATERIAL=membr
15 2
16
17 ** Boundary conditions to be applied to the infinitely rigid edges.
18 *BOUNDARY
19 Ncables, 1, 3
20
21 ** Static nonlinear step for the analysis: request a fixed time increment,
22 ** apply the distributed load and retrieve nodal displacements and
23 ** reaction forces, and element stresses and strains.
24 *STEP, NLGEOM
25 *STATIC, SOLVER=SPOOLES
26 1.0e-6, 1.0, , 0.05
27 * DLOAD
28 Ememb, P, -0.001
29 * NODE FILE
30 U, RF
31 *EL FILE
32 S, E
33 *END STEP
8.6.4 Post-processing
The post-processing is automated using cgx to extract from the solution file the contour plots
of the displacements (magnitude, x and y-direction), stresses and strains (x and y-direction and
Mises). We focus our attenion also on the time graphs (i.e. load multiplier) which are drawn
for three significant points:
1. mid point of the membrane representing the point of maximum displacement.
2. pin point for the analysis of a support.
3. mid point of an edge which is usually subject to the maximum stress.
The same example has been developed in the commercial software Dlubal to test the three quan-
tities of interest: displacements, stresses and strains. The results are reported in the attached
figures.
8.6.6 Conclusions
As far as the vertical displacements are concerned, the behaviour displayed by the membrane is
strongly nonlinear in the initial phase, where the convergence of each step is difficult due to the
lack of pretensioning. After the load multiplier reaches the value of 0.3, the membrane regains
sufficient stiffening to substain additional loading without excessive further deformation and
the behaviour becomes almost linear.
As for the stresses, the stress in the coordinate directions for the center and mid points are
slightly nonlinear for all the analysis, differently from the trend of displacements, where the
initial contribution deviates strongly from linearity.
An important remark for the stresses in the support is the initial tensioning of the membrane
in the local zone, which however inverts its trend after a load multiplier equal to 0.10. The
subsequent minor compression of the membrane can be explained considering the same case
for pinned plate on four edges. This behaviour disappears, as shown later, when the structure
is pre-processed by the form-finding procedure.
Finally, the Von-Mises stresses present a similar graph with respect to the coordinate stresses,
but the nonlinearity is accentuated.
1
0.9
Node=842
0.8
0.7
Values at Nodes (membrane-plane.fbd)
0.6
Time
0.5
0.4
0.3
0.2
0.1
0
0
-10
-20
-30
-40
-50
-60
-70
-80
-90
-100
D3
10
6
SXX
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
10
6
SYY
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
8.46: Normal stresses in the x and y-direction of the center point of the plane membrane.
12
10
8
SXX
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
4.5
3.5
2.5
SYY
1.5
0.5
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
8.47: Normal stresses in the x and y-direction of the mid-point of the plane membrane.
0.1
0.05
SXX
-0.05
-0.1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
0.1
0.08
0.06
0.04
0.02
SYY
-0.02
-0.04
-0.06
-0.08
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
8.48: Normal stresses in the x and y-direction of the pin of the plane membrane.
0.0008
0.0007
0.0006
0.0005
EXX
0.0004
0.0003
0.0002
0.0001
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
0.0008
0.0007
0.0006
0.0005
EYY
0.0004
0.0003
0.0002
0.0001
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
8.49: Normal strains in the x and y-direction of the center point of the plane membrane.
0.0012
0.001
0.0008
EXX
0.0006
0.0004
0.0002
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
-1x10-6
-2x10-6
-3x10-6
EYY
-4x10-6
-5x10-6
-6x10-6
-7x10-6
-8x10-6
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
8.50: Normal strains in the x and y-direction of the mid-point of the plane membrane.
-5x10-6
-1x10-5
EXX
-1.5x10-5
-5
-2x10
-2.5x10-5
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
-5x10-6
EYY
-1x10-5
-1.5x10-5
-2x10-5
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
8.51: Normal strains in the x and y-direction of the pin of the plane membrane.
10
8
Mises
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
0.0012
0.001
0.0008
Mises
0.0006
0.0004
0.0002
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
8.52: Von-Mises stresses and strains of the center point of the plane membrane.
10
8
Mises
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
0.0012
0.001
0.0008
Mises
0.0006
0.0004
0.0002
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
8.53: Von-Mises stresses and strains of the mid-point of the plane membrane.
1.4
1.2
1
Mises
0.8
0.6
0.4
0.2
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
0.00018
0.00016
0.00014
0.00012
Mises
0.0001
8x10-5
6x10-5
4x10-5
2x10-5
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
8.54: Von-Mises stresses and strains of the pin of the plane membrane.
8.7.1 Intoduction
The logical conseguence of the previous analysis is to improve the model in order to better
represent a real case. In fact, as described in chapter 1.2, the memebrane presents a typical
orthotropic mechanical behaviour.
The analysis follows the same principles adopted for the plane membrane, to which an extensive
comparison is drawn in the conclusion. The only modification is the elastic modulus in the y-
direction is reduced to 4 · 103 𝑀 𝑃 𝑎 and the Poisson’s ratio accordingly to 0.2.
8.7.2 Pre-processing
The actual mesh is created with open-source software gmsh and imported in the Calculix ap-
plication cgx for the pre-processing; at the same time, the solution and the post-processing is
conducted with the reported bash-like script with .fbd extension.
˓→orthotropy_displacement_xy.png
74
75 ds -4 e 3
76 plot fv all
77 max 0.0
78 min -110
79 sys sleep 1
80 hcpy png
81 sys mv hcpy_6.png membrane_plane_orthotropy_displacement_z.png
82
83 ds -3 e 1
84 plot fv all
(continues on next page)
125
126 ds -2 e 7
127 plot fv all
128 max 1.75e-3
129 sys sleep 1
130 hcpy png
131 sys mv hcpy_12.png membrane_plane_orthotropy_strain_m.png
132
133 sys montage membrane_plane_orthotropy_strain_m.png membrane_plane_
˓→orthotropy_stress_m.png -tile 2x1 -geometry +0+0 membrane_plane_
˓→orthotropy_ss_m.png
134
(continues on next page)
˓→pdf membrane_plane_orthotropy_strain_yy_frex.pdf
158
159 graph frex time TOSTRAIN Mises
160 sys mv graph_6.ps membrane_plane_orthotropy_strain_mises_frex.eps
161
162 sys ps2pdf membrane_plane_orthotropy_stress_mises_frex.eps
163 sys ps2pdf membrane_plane_orthotropy_strain_mises_frex.eps
164 sys pdfnup --no-landscape --a4paper --nup 1x2 --outfile membrane_plane_
˓→orthotropy_se_frex.pdf membrane_plane_orthotropy_stress_mises_frex.pdf
˓→membrane_plane_orthotropy_strain_mises_frex.pdf
165
166 seta pin n 536
167 graph pin time STRESS SXX
168 sys mv graph_7.ps membrane_plane_orthotropy_stress_xx_pin.eps
169 graph pin time STRESS SYY
170 sys mv graph_8.ps membrane_plane_orthotropy_stress_yy_pin.eps
171
172 sys ps2pdf membrane_plane_orthotropy_stress_xx_pin.eps
173 sys ps2pdf membrane_plane_orthotropy_stress_yy_pin.eps
174 sys pdfnup --no-landscape --a4paper --nup 1x2 --outfile membrane_plane_
˓→orthotropy_stress_xxyy_pin.pdf membrane_plane_orthotropy_stress_xx_pin.
˓→pdf membrane_plane_orthotropy_stress_yy_pin.pdf
175
176 graph pin time STRESS Mises
177 sys mv graph_9.ps membrane_plane_orthotropy_stress_mises_pin.eps
178 graph pin time TOSTRAIN EXX
179 sys mv graph_10.ps membrane_plane_orthotropy_strain_xx_pin.eps
180 graph pin time TOSTRAIN EYY
181 sys mv graph_11.ps membrane_plane_orthotropy_strain_yy_pin.eps
182
(continues on next page)
˓→pdf membrane_plane_orthotropy_strain_yy_mid.pdf
214
215 graph mid time TOSTRAIN Mises
216 sys mv graph_18.ps membrane_plane_orthotropy_strain_mises_mid.eps
217
218 sys ps2pdf membrane_plane_orthotropy_stress_mises_pin.eps
219 sys ps2pdf membrane_plane_orthotropy_strain_mises_pin.eps
220 sys pdfnup --no-landscape --a4paper --nup 1x2 --outfile membrane_plane_
˓→orthotropy_se_pin.pdf membrane_plane_orthotropy_stress_mises_pin.pdf
˓→membrane_plane_orthotropy_strain_mises_pin.pdf
221
222
223 # switch back to the deformed shape
224 rot -z
225 rot u -60
226 rot c -45
227 view disp
228 ds -4 e 4
(continues on next page)
8.7.3 Solution
The orthotropic paramters are specified in the deck file in the section concerning the elastic
properties of the membrane. The deck then follows as usual instructing the nonlinear Newton
solver for the iterative process.
Code block 8.33: Calculix deck for the solution of the or-
thotropic membrane.
1 ** Include the mesh produced by gmsh.
2 *INCLUDE, INPUT=all.msh
3
4 ** Include the physical entities produced by the .fbd file.
5 *INCLUDE, INPUT=memb.nam
6 *INCLUDE, INPUT=cables.nam
7
8 ** Orthotropic elastic material for the membrane.
9 *MATERIAL, NAME=membrorto
10 *ELASTIC, TYPE=ENGINEERING CONSTANTS
11 8e3, 4e3, 8e3, 0.4, 0.2, 0.4, 2.9e3, 2.9e3
12 2.9e3
13
14 ** Define the thickness for the membrane and assign the material.
15 *SHELL SECTION, ELSET=Ememb, MATERIAL=membrorto
16 2
17
Code block 8.34: Bash script to start Calculix from the com-
mand line for the orthotropic membrane.
1 #!/usr/bin/env bash
2 cgx_2.14 -b membrane-plane-orthotropy.fbd
8.7.4 Post-processing
The post-processing is modified to force the contour plots to report the same maximum and
minimum scale of the previous example in order to draw a thouroughly comparison with the
linear isotropic case.
As benchmark, the same model improvement has been introduced in the commercial software
Dlubal, to verify accuracy of the results in terms of displacements (magnitude), stresses and
strains (x and y-direction, Mises).
8.7.6 Conclusions
The most important conclusion coming from the analysis at hand is that the orthotropy does not
change the pattern of the displacements, stresses and strains, if the membrane is loaded out of
8.62: Normal stresses in the x and y-direction of plane orthotropic plane membrane.
8.63: Normal strains in the x and y-direction of the orthotropic plane membrane.
its plane. This behaviour is essentially different from the structural case of a plane membrane
loaded in its plane, where the orthotropy causes the stress pattern to differentiate in the principal
directions.
The orthotropy is mainly reflected in the value of the maximum and minimum attained vari-
ables for the three quantities of interest reported. The most significant deviation is for the
z-displacement, as shown by the time graphs in the previous paragraph.
1
0.9
Node=842
0.8
Values at Nodes (membrane-plane-orthotropy.fbd)
0.7
0.6
Time
0.5
0.4
0.3
0.2
0.1
0
0
-20
-40
-60
-80
-100
-120
D3
10
6
SXX
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
4
SYY
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
8.66: Normal stresses in the x and y-direction of the center point of the orthotropic membrane.
12
10
8
SXX
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
2.5
1.5
SYY
0.5
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
8.67: Normal stresses in the x and y-direction of the mid-point of the orthotropic membrane.
-0.1
-0.2
-0.3
-0.4
SXX
-0.5
-0.6
-0.7
-0.8
-0.9
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
-0.05
-0.1
-0.15
-0.2
SYY
-0.25
-0.3
-0.35
-0.4
-0.45
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
8.68: Normal stresses in the x and y-direction of the pin of the orthotropic membrane.
0.001
0.0008
0.0006
EXX
0.0004
0.0002
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
0.001
0.0008
0.0006
EYY
0.0004
0.0002
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
8.69: Normal strains in the x and y-direction of the center-point of the orthotropic membrane.
0.0014
0.0012
0.001
0.0008
EXX
0.0006
0.0004
0.0002
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
-1x10-6
-2x10-6
-3x10-6
EYY
-4x10-6
-5x10-6
-6
-6x10
-7x10-6
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
8.70: Normal strains in the x and y-direction of the mid-point of the orthotropic membrane.
-1x10-5
-2x10-5
-3x10-5
-4x10-5
EXX
-5x10-5
-6x10-5
-7x10-5
-8x10-5
-9x10-5
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
-1x10-5
-2x10-5
-3x10-5
-4x10-5
EYY
-5x10-5
-6x10-5
-7x10-5
-8x10-5
-9x10-5
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
8.71: Normal strains in the x and y-direction of the pin of the orthotropic membrane.
6
Mises
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
0.0012
0.001
0.0008
Mises
0.0006
0.0004
0.0002
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
8.72: Von-Mises stresses and strains of the center point of the orthotropic membrane.
12
10
8
Mises
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
0.0012
0.001
0.0008
Mises
0.0006
0.0004
0.0002
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
8.73: Von-Mises stresses and strains of the mid-point of the orthotropic membrane.
1.5
Mises
0.5
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
0.00025
0.0002
Mises
0.00015
0.0001
5x10-5
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
8.74: Von-Mises stresses and strains of the pin of the orthotropic membrane.
8.8.1 Introduction
We present here the same case of the plane membrane with distributed pin constraints on the
borders, changing the mechanical behaviour of the membrane to elastoplastic. The equivalent
plastic strains are assigned directly in the Calculix deck by defining a table linking the level of
stress to the Von Mises equivalent plastic strain.
The membrane is subject to a considerable load which extensively plasticizes it; the load is then
removed to observe the permanent plastic strains.
8.8.2 Pre-processing
The simple square membrane is generated with gmsh to produce both the geometrical entitiy
identifying the plastic membrane and the physical entities defining the borders representing the
infinitely rigid cables.
Code block 8.35: Gmsh script to produce the mesh for the
plastic plane membrane.
1 // Characteristic length of the triangular mesh to be created on the
2 // membrane.
3 lc = 500.0;
4
5 // Length of the edge of the square membrane.
6 l = 5000.0;
7
8 // Geometrical points defining the pins where the membrane is anchored.
9 Point(1) = {0, 0, 0, lc};
10 Point(2) = {l, 0, 0, lc};
11 Point(3) = {l, l, 0, lc};
12 Point(4) = {0, l, 0, lc};
13
14 // Geometrical lines defining the edges or boundaries of the membrane.
15 Line(1) = {1, 2};
16 Line(2) = {2, 3};
17 Line(3) = {3, 4};
18 Line(4) = {4, 1};
19
20 // Line loop defining the edges which enclose the membrane.
21 Line Loop(1) = {1, 2, 3, 4};
22
23 // Geometrical surface identifying the membrane.
24 Surface(1) = {1};
25
26 // Physical line representing the edges, as though they were infinitely
27 // rigid cables
28 Physical Line("cables", 200) = {1, 2, 3, 4};
29
30 // Physical surface representing the membrane which will be meshed.
31 Physical Surface("memb", 300) = {1};
(continues on next page)
63 ds -5 e 2
64 plot fv all
65 sys sleep 1
66 hcpy png
67 sys mv hcpy_5.png membrane_plane_plasticity_displacement_y.png
68
˓→plasticity_ss_m.png
118
119 ds -1 e 1
120 plot fv all
121 sys sleep 1
122 hcpy png
123 sys mv hcpy_13.png membrane_plane_plasticity_pe.png
124
125 # time graphs
126 seta frex n 842
127 graph frex time DISP D3
128 sys mv graph_0.ps membrane_plane_plasticity_displacement_frex.eps
129 graph frex time STRESS SXX
130 sys mv graph_1.ps membrane_plane_plasticity_stress_xx_frex.eps
131 graph frex time STRESS SYY
132 sys mv graph_2.ps membrane_plane_plasticity_stress_yy_frex.eps
133
˓→pdf membrane_plane_plasticity_stress_yy_frex.pdf
137
138 graph frex time STRESS Mises
139 sys mv graph_3.ps membrane_plane_plasticity_stress_mises_frex.eps
140 graph frex time TOSTRAIN EXX
141 sys mv graph_4.ps membrane_plane_plasticity_strain_xx_frex.eps
142 graph frex time TOSTRAIN EYY
143 sys mv graph_5.ps membrane_plane_plasticity_strain_yy_frex.eps
144
145 sys ps2pdf membrane_plane_plasticity_strain_xx_frex.eps
146 sys ps2pdf membrane_plane_plasticity_strain_yy_frex.eps
147 sys pdfnup --no-landscape --a4paper --nup 1x2 --outfile membrane_plane_
˓→plasticity_strain_xxyy_frex.pdf membrane_plane_plasticity_strain_xx_frex.
˓→pdf membrane_plane_plasticity_strain_yy_frex.pdf
148
149 graph frex time TOSTRAIN Mises
150 sys mv graph_6.ps membrane_plane_plasticity_strain_mises_frex.eps
151
˓→pdf membrane_plane_plasticity_stress_yy_mid.pdf
199
200 graph mid time STRESS Mises
201 sys mv graph_17.ps membrane_plane_plasticity_stress_mises_mid.eps
202 graph mid time TOSTRAIN EXX
203 sys mv graph_18.ps membrane_plane_plasticity_strain_xx_mid.eps
204 graph mid time TOSTRAIN EYY
205 sys mv graph_19.ps membrane_plane_plasticity_strain_yy_mid.eps
206
207 sys ps2pdf membrane_plane_plasticity_strain_xx_mid.eps
208 sys ps2pdf membrane_plane_plasticity_strain_yy_mid.eps
209 sys pdfnup --no-landscape --a4paper --nup 1x2 --outfile membrane_plane_
˓→plasticity_strain_xxyy_mid.pdf membrane_plane_plasticity_strain_xx_mid.
˓→pdf membrane_plane_plasticity_strain_yy_mid.pdf
210
211 graph mid time TOSTRAIN Mises
212 sys mv graph_20.ps membrane_plane_plasticity_strain_mises_mid.eps
213
214 sys ps2pdf membrane_plane_plasticity_stress_mises_pin.eps
215 sys ps2pdf membrane_plane_plasticity_strain_mises_pin.eps
216 sys pdfnup --no-landscape --a4paper --nup 1x2 --outfile membrane_plane_
˓→plasticity_se_pin.pdf membrane_plane_plasticity_stress_mises_pin.pdf
˓→membrane_plane_plasticity_strain_mises_pin.pdf
217
8.8.3 Solution
The plastic behaviour, as previously done for the cables, is defined via the plastic card in the
material properties of the membrane: only one table is defined, so that the properties of the
material are independent of the temperature.
8.8.4 Post-processing
The post-processing is fully automated by the previous bash-like script written for cgx in Cal-
culix.
8.8.5 Conclusions
The first important conclusion to be drawn in the structural case of the plastic membrane is that
the membrane should never reach the plastic domain; the plastic permanent strain deforms the
2.5
Node=842
2
Values at Nodes (membrane-plane-plasticity.fbd)
1.5
Time
1
0.5
0
0
-20
-40
-60
-80
-100
-120
-140
-160
-180
-200
D3
30
25
20
SXX
15
10
0
0 0.5 1 1.5 2 2.5
Time
30
25
20
SYY
15
10
0
0 0.5 1 1.5 2 2.5
Time
8.87: Normal stresses in the x and y-direction of the center point of the plastic membrane.
30
25
20
SXX
15
10
0
0 0.5 1 1.5 2 2.5
Time
16
14
12
10
SYY
0
0 0.5 1 1.5 2 2.5
Time
8.88: Normal stresses in the x and y-direction of the mid-point of the plastic membrane.
14
12
10
8
SXX
-2
0 0.5 1 1.5 2 2.5
Time
30
25
20
15
SYY
10
-5
0 0.5 1 1.5 2 2.5
Time
8.89: Normal stresses in the x and y-direction of the pin of the plastic membrane.
0.0035
0.003
0.0025
0.002
EXX
0.0015
0.001
0.0005
0
0 0.5 1 1.5 2 2.5
Time
0.0035
0.003
0.0025
0.002
EYY
0.0015
0.001
0.0005
0
0 0.5 1 1.5 2 2.5
Time
8.90: Normal strains in the x and y-direction of the center point of the plastic membrane.
0.007
0.006
0.005
0.004
EXX
0.003
0.002
0.001
0
0 0.5 1 1.5 2 2.5
Time
-1x10-5
-2x10-5
EYY
-3x10-5
-4x10-5
-5x10-5
-6x10-5
0 0.5 1 1.5 2 2.5
Time
8.91: Normal strains in the x and y-direction of the mid-point of the plastic membrane.
-1x10-5
-5
-2x10
-3x10-5
EXX
-4x10-5
-5
-5x10
-6x10-5
-7x10-5
0 0.5 1 1.5 2 2.5
Time
0.006
0.005
0.004
EYY
0.003
0.002
0.001
0
0 0.5 1 1.5 2 2.5
Time
8.92: Normal strains in the x and y-direction of the pin of the plastic membrane.
30
25
20
Mises
15
10
0
0 0.5 1 1.5 2 2.5
Time
0.006
0.005
0.004
Mises
0.003
0.002
0.001
0
0 0.5 1 1.5 2 2.5
Time
8.93: Von-Mises stresses and strains of the center-point of the plastic membrane.
30
25
20
Mises
15
10
0
0 0.5 1 1.5 2 2.5
Time
0.007
0.006
0.005
Mises
0.004
0.003
0.002
0.001
0
0 0.5 1 1.5 2 2.5
Time
8.94: Von-Mises stresses and strains of the mid-point of the plastic membrane.
30
25
20
Mises
15
10
0
0 0.5 1 1.5 2 2.5
Time
0.006
0.005
0.004
Mises
0.003
0.002
0.001
0
0 0.5 1 1.5 2 2.5
Time
8.95: Von-Mises stresses and strains of the pin of the plastic membrane.
2.5
Node=842
2
Values at Nodes (membrane-plane-plasticity.fbd)
1.5
Time
1
0.5
0
0.003
0.0025
0.002
0.0015
0.001
0.0005
PE
8.96: Von-Mises equivalent plastic strain of the center point of the membrane.
2.5
Node=611
2
Values at Nodes (membrane-plane-plasticity.fbd)
1.5
Time
1
0.5
0
0.004
0.0035
0.003
0.0025
0.002
0.0015
0.001
0.0005
PE
2.5
Node=553
2
Values at Nodes (membrane-plane-plasticity.fbd)
1.5
Time
1
0.5
0
0.003
0.0025
0.002
0.0015
0.001
0.0005
PE
membrane so extensively that not only the structural properties are compromised, but also the
aesthetic characteristic are significantly deteriorated.
In addition, the membrane is significantly compressed in the zones adjacent to the pins: the
plastic deformations have made the fabric slack; we should also question whether the material
can bear this compressive stresses or we should modify the constitutive law to cut them off.
Finally, observing the graphs of the Von Mises equivalent plastic strain, all the membrane finite
elements reach the plastic domain and present permanent deformations after the discharge.
8.9.1 Introduction
In a practical contest, membranes are enhanced with sourrounding cables which are responsible
for the application of the total initial prestress to the entire mixed cable and membrane structure.
The characteristics of the membrane are the same ones as reported in the case of the plane
membrane pinned along its borders (chapter 8.6). As far as the cables are concerned, the
Young modulus of the steel is equal to 160 · 103 𝑀 𝑃 𝑎 with a sectional area of 314𝑚𝑚2 , with a
pretension (force, and not force density) of 60𝑘𝑁 .
The pretension is applied at the initial instant, whereas the load is increased monotonically until
time 1.0, when the pretension is further increased to show the originated stiffening behaviour.
8.9.2 Pre-processing
gmsh provides the user with all the necessary cabilities to create meshes with elements of
different geometric order; for instance, we can mesh at the same time both 1D cable elems and
2D membranes which may also share common nodes.
The initial geometry stems from an approximate form-finding calculated by hand and also
verified with Dlubal. Additional details for the form-finding can be found in the next chapters.
8.99: Mesh produced by Calculix of the plane membrane with collaborating cables.
106
107 ds -4 e 3
108 plot fv all
109 sys sleep 1
110 hcpy png
111 sys mv hcpy_6.png membrane_plane_cables_displacement_z.png
112
113 ds -3 e 1
114 plot fv all
115 min 20
116 max 50
117 sys sleep 1
118 hcpy png
119 sys mv hcpy_7.png membrane_plane_cables_stress_xx.png
120
121 ds -3 e 2
122 plot fv all
123 min 20
124 max 50
125 sys sleep 1
126 hcpy png
127 sys mv hcpy_8.png membrane_plane_cables_stress_yy.png
128
153 ds -2 e 7
154 plot fv all
155 sys sleep 1
156 hcpy png
157 sys mv hcpy_12.png membrane_plane_cables_strain_m.png
158
159 sys montage membrane_plane_cables_strain_m.png membrane_plane_cables_
˓→stress_m.png -tile 2x1 -geometry +0+0 membrane_plane_cables_ss_m.png
160
161
162 # time graphs
163 seta frex n 2966
164 graph frex time DISP D3
165 sys mv graph_0.ps membrane_plane_cables_displacement_frex.eps
166 graph frex time STRESS SXX
167 sys mv graph_1.ps membrane_plane_cables_stress_xx_frex.eps
168 graph frex time STRESS SYY
169 sys mv graph_2.ps membrane_plane_cables_stress_yy_frex.eps
170
171 sys ps2pdf membrane_plane_cables_stress_xx_frex.eps
172 sys ps2pdf membrane_plane_cables_stress_yy_frex.eps
173 sys pdfnup --no-landscape --a4paper --nup 1x2 --outfile membrane_plane_
˓→cables_stress_xxyy_frex.pdf membrane_plane_cables_stress_xx_frex.pdf
˓→membrane_plane_cables_stress_yy_frex.pdf
174
175 graph frex time STRESS Mises
176 sys mv graph_3.ps membrane_plane_cables_stress_mises_frex.eps
177 graph frex time TOSTRAIN EXX
178 sys mv graph_4.ps membrane_plane_cables_strain_xx_frex.eps
179 graph frex time TOSTRAIN EYY
180 sys mv graph_5.ps membrane_plane_cables_strain_yy_frex.eps
181
182 sys ps2pdf membrane_plane_cables_strain_xx_frex.eps
183 sys ps2pdf membrane_plane_cables_strain_yy_frex.eps
184 sys pdfnup --no-landscape --a4paper --nup 1x2 --outfile membrane_plane_
(continues on next page)
˓→cables_strain_xxyy_frex.pdf membrane_plane_cables_strain_xx_frex.pdf
˓→membrane_plane_cables_strain_yy_frex.pdf
192
193 seta pin n 2426
194 graph pin time STRESS SXX
195 sys mv graph_7.ps membrane_plane_cables_stress_xx_pin.eps
196 graph pin time STRESS SYY
197 sys mv graph_8.ps membrane_plane_cables_stress_yy_pin.eps
198
199 sys ps2pdf membrane_plane_cables_stress_xx_pin.eps
200 sys ps2pdf membrane_plane_cables_stress_yy_pin.eps
201 sys pdfnup --no-landscape --a4paper --nup 1x2 --outfile membrane_plane_
˓→cables_stress_xxyy_pin.pdf membrane_plane_cables_stress_xx_pin.pdf
˓→membrane_plane_cables_stress_yy_pin.pdf
202
203 graph pin time STRESS Mises
204 sys mv graph_9.ps membrane_plane_cables_stress_mises_pin.eps
205 graph pin time TOSTRAIN EXX
206 sys mv graph_10.ps membrane_plane_cables_strain_xx_pin.eps
207 graph pin time TOSTRAIN EYY
208 sys mv graph_11.ps membrane_plane_cables_strain_yy_pin.eps
209
210 sys ps2pdf membrane_plane_cables_strain_xx_pin.eps
211 sys ps2pdf membrane_plane_cables_strain_yy_pin.eps
212 sys pdfnup --no-landscape --a4paper --nup 1x2 --outfile membrane_plane_
˓→cables_strain_xxyy_pin.pdf membrane_plane_cables_strain_xx_pin.pdf
˓→membrane_plane_cables_strain_yy_pin.pdf
213
214 graph pin time TOSTRAIN Mises
215 sys mv graph_12.ps membrane_plane_cables_strain_mises_pin.eps
216
217 sys ps2pdf membrane_plane_cables_stress_mises_mid.eps
218 sys ps2pdf membrane_plane_cables_strain_mises_mid.eps
219 sys pdfnup --no-landscape --a4paper --nup 1x2 --outfile membrane_plane_
˓→cables_se_mid.pdf membrane_plane_cables_stress_mises_mid.pdf membrane_
˓→plane_cables_strain_mises_mid.pdf
220
221 seta mid n 2939
222 graph mid time STRESS SXX
223 sys mv graph_13.ps membrane_plane_cables_stress_xx_mid.eps
224 graph mid time STRESS SYY
225 sys mv graph_14.ps membrane_plane_cables_stress_yy_mid.eps
226
248
249
250 # switch back to the deformed shape
251 rot -z
252 rot u -60
253 rot c -45
254 view disp
255 ds -4 e 4
256 plot fv all
257
258
8.9.3 Solution
We would like to underline the fact that the numerical solution of the system of equations for the
nonlinear static analysis is delegated to the open-source solver SPOOLES, which is extremely
handy at dealing with the strong nonlinearities stemming from the collaboration of the two
distinct elements.
34 ** The membrane and cables are anchored at the four side pins.
35 *BOUNDARY
36 Npins, 1, 3
37
Code block 8.46: Bash script to start Calculix for the analy-
sis of the plane membrane with collaborating cables.
1 cgx_2.14 -b membrane-plane-cables.fbd
8.9.4 Post-processing
In order to plot only the stresses of the membrane, being the latter of a different order of
magnitude compared to the cables. There are two choices to tackle the problem:
1. Suppress the output for the cables and show only the “S6” elements which represent the
finite element of the membrane.
2. Cap the maximum value of the quantity of interest thus excluding the results originating
from the collaborating cables.
With the commercial software Dlubal we first perform the form finding of the structure, then
load it and finally increase the pretensioning. The results are comparable with the ones obtained
in Calculix, where the form-finding has been delegated, for the time being, to calculations by
hand.
8.102: Normal stresses in the x and y-direction of the membrane with cables.
8.103: Normal strains in the x and y-direction of the membrane with cables.
2
1.8
Node=2966
1.6
Values at Nodes (membrane-plane-cables.fbd)
1.4
1.2
Time
1
0.8
0.6
0.4
0.2
0
20
-20
-40
-60
-80
-100
-120
D3
30
25
20
SXX
15
10
0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time
30
25
20
SYY
15
10
0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time
8.106: Normal stresses in the x and y-direction of the center point of the cabled membrane.
35
30
25
20
SXX
15
10
0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time
35
30
25
20
SYY
15
10
0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time
8.107: Normal stresses in the x and y-direction of the mid-point of the cabled membrane.
35
30
25
20
SXX
15
10
0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time
35
30
25
20
SYY
15
10
0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time
8.108: Normal stresses in the x and y-direction of the pin of the cabled membrane.
0.002
0.0015
EXX
0.001
0.0005
0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time
0.002
0.0015
EYY
0.001
0.0005
0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time
8.109: Normal strains in the x and y-direction of the center point of the cabled membrane.
6x10-5
5x10-5
4x10-5
3x10-5
EXX
2x10-5
1x10-5
-1x10-5
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time
0.0002
0.00015
0.0001
EYY
5x10-5
-5x10-5
-0.0001
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time
8.110: Normal strains in the x and y-direction of the mid-point of the cabled membrane.
-0.0005
-0.001
-0.0015
EXX
-0.002
-0.0025
-0.003
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time
-0.0005
-0.001
-0.0015
EYY
-0.002
-0.0025
-0.003
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time
8.111: Normal strains in the x and y-direction of the pin of the cabled membrane.
30
25
20
Mises
15
10
0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time
0.0035
0.003
0.0025
Mises
0.002
0.0015
0.001
0.0005
0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time
8.112: Von-Mises stresses and strains of the center point of the membrane with cables.
35
30
25
Mises
20
15
10
0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time
0.004
0.0035
0.003
0.0025
Mises
0.002
0.0015
0.001
0.0005
0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time
8.113: Von-Mises stresses and strains of the mid-point of the membrane with cables.
40
35
30
25
Mises
20
15
10
0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time
0.0045
0.004
0.0035
0.003
Mises
0.0025
0.002
0.0015
0.001
0.0005
0
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
Time
8.114: Von-Mises stresses and strain of the pin of the membrane with cables.
8.115: Check of the deformations and the 𝜎x of the membrane with cables.
8.116: Check of the 𝜎y and the 𝜎mises of the membrane with cables.
8.118: Check of the 𝜀mises of the membrane and the tension in the cables.
8.9.6 Conclusions
The example of a plane membrane with collaborating cables shows the importance of the form-
finding procedure to identify a suitable shape of the membrane to allow the correct pretension-
ing of the membrane itself. The cables, which adhere to the border of the membrane, must
be positioned such that the membrane is properly tensioned; the additional prerequisite of the
form-finding to impose an almost uniform stress state guarantees also that the membrane will
be devoid of slack portions and most importantly will react practically super-linearly to the
additional live loads.
This procedure has been developed by Bernard Maurin and René Motro in 1997 in order to pro-
vide a valid method which implements the form-finding theory to tensile membranes combined
with reinforcing cables at the edges or above the fabric [MaMo211].
Both the Force Density Method (Chapter 2) and the Surface Stress Density Method (Chapter 3)
are based on the same solving approach: nodal equilibrium. Therefore, the following procedure
will compute the nodal load due to the membrane and the cables and will set the sum of the
two equal to zero.
n
∑︁ m
∑︁
Rres = Ti + Fj (9.1)
i=1 j=1
The code is structured in different files whose classes and functions are responsible for:
• main.py: defines the interface of the program and starts the execution of the form-finding
procedure, after having parsed the command line arguments.
• parser.py: parses the input and output files, and writes the tests in hdf5 binary format for
later execution by the test suite in pytest.
• ff.py: is the actual form-finding algorithm on top of the Newton-Krylov algorithm, which
solves iteratively the nonlinear system of equations which compose the problem.
247
Form-Finding and Mechanical Behaviour of Cable and Membrane Structures
The form finding algorithm is invoked using the command line program written in Python based
on the following interface:
The main function accepts only the arguments passed by the parser and starts the execution of
the program:
ff_clc.main(args)
Main entry point for the form-finding algorithm for membrane and cable structures.
Takes care of passing the control to the functions delegated to the pre-processing of the
mesh and problem files, to the actual form-finding algorithm, and finally to the post-
processor to write the output for the study with additional loadings and nonlinear consti-
tutive laws. Optionally, the user can request a test file to be created for future reference
or inspection.
Parameters args (list) – arguments parsed by the command line parser.
Returns system exit code. 0 if successful.
The parser is written with argparse, a module that makes it easy to write user-friendly
command-line interfaces. We define which arguments we require, and argparse will figure
out how to parse those out of sys.argv.
Code block 9.3: Parser for the geometry and problem files.
1 import numpy as np
2 import h5py
3 import os
4
5
6 def read_topology(file_path, prob_path):
7 # List container for the nodes in matrix form, where the index
8 # represents the node number minus one, and the columns the x,
9 # y and z-coordinates of the node
10 nodes = []
11 # List container for the cable elements in matrix form, where the
12 # first column represents the element number, whereas the second
13 # and third columns are the nodes connected by the element itself.
14 cable = []
15 # List container for the membrane elements in matrix form, where the
16 # first column represents the element number, whereas the second,
17 # third and fourth are the nodes connected by a three-nodes membrane
18 # plane first order membrane element.
19 membr = []
20 # Set containing the boundary conditions, instead of list to speed
21 # up the lookup during the definition of the boundary equations.
22 bcc = set()
23
24 # Open the file containing the topology and mesh.
25 with open(file_path, 'r') as fin:
26 # Enter the file and retrieve the first line.
27 line = next(fin)
28 # Skip the heading lines.
29 for ii in range(0, 3):
30 line = next(fin)
31 # The first part of the file contains the node definitions to be
32 # passed to the node container. Read until we meet the header
33 # of the element section.
34 while not line.startswith("******* E L E M E N T S *************"):
35 # Split the line according to the formatting of Abaqus and
36 # retrieve the node coordinates.
37 checked_line = line.split(', ')
38 nodes.append([float(checked_line[1]), float(checked_line[2]),
39 float(checked_line[3])])
40 line = next(fin)
41 # Advance a couple of lines to enter the element section.
42 line = next(fin)
43 line = next(fin)
44 # Go on until we find the physical entities which we are not
45 # interested in for the time being.
46 while not line.startswith("*ELSET"):
47 # If a card is placed at the beginning of the line, that is a
48 # header for the set.
49 if not line.startswith("*ELEMENT"):
50 # Split the line according to the Abaqus format and
51 # retrieve the needed information.
52 checked_line = line.split(', ')
53 checked_line_len = len(checked_line)
54 # If the line consists of three entries (element number,
(continues on next page)
The form-finding topology is written with the following function of the parser:
Returns None
Raises
• FileNotFoundError – if any of the input files is not found.
• PermissionError – insufficient permissions to write to the cor-
responding path.
Code block 9.4: Parser for the output of the updated geome-
try and new deck in Calculix
1 def write_topology(input_path, solve_path, nodes, cable, prec,
2 # Based on the name of the file supplied by the user, output the
3 # resulting topology from the form-finding process to the same
4 # directory appending '_resu' to the original name of the file.
5 working_dir_input = os.path.dirname(input_path)
6 output_file = os.path.join(working_dir_input, "form_finding_resu.inp")
7 # At the same time, open the input file with the original topology
8 # and a blank file to write the new geometry after the form-finding
9 # process.
10 with open(input_path, "r") as fp, open(output_file, 'w') as fo:
11 # Initialize the node counter to keep track of the written nodes
12 # in the output file.
13 node_counter = 0
14 line = next(fp) # get the first line of the file
15 fo.write(line)
16 # Copy the header of the input file to the output file.
17 for ii in range(0, 3):
18 line = next(fp)
19 fo.write(line)
20 # Until we reach the elements cards, copy the node to the new
21 # destination with the updated final coordinates.
22 while not line.startswith("******* E L E M E N T S *************"):
23 fo.write("%d, %f, %f, %f\n" % (
24 node_counter + 1, nodes[node_counter, 0],
25 nodes[node_counter, 1], nodes[node_counter, 2]))
26 node_counter += 1 # advance the counter for the node numbering
27 line = next(fp)
28 fo.write(line)
29 # Write everything else to the file in its original form,
30 # in particular the element connections and the physical groups.
31 for line in fp:
32 fo.write(line)
33
34 # Retrieve the number of nodes and cable elements.
35 nodes_len = nodes.shape[0]
36 cable_len = cable.shape[0]
37 # Pre-allocate an empty vector with the nodal temperatures
38 # representing the inelastic contribution of the prestressing.
39 delta_t = np.zeros(nodes_len)
40 for ii in range(0, cable_len):
41 # Extract the starting and end node for the cable elements.
42 start_cable = cable[ii, 1] - 1
43 end_cable = cable[ii, 2] - 1
44 # For the current cable, compute the length in the deformed
45 # configuration.
(continues on next page)
Code block 9.6: Test writing for the output of the verification
with the hdf5 file.
1 def write_test(test_path, nodes, cable, membr, bcc, pretension_cable,
2 # Create a blank hdf5 file to write the initial topology and static
3 # conditions. The final coordinates after the form-finding process
4 # are added to allow pytest to test without relying on the custom
5 # parser written for Calculix; hdf5 acts as transparent layer to
6 # abstract the input and output operations and also to speed them up
7 # for extremely large files.
8 # 'libver' is set to 'latest' to gain maximum speed to skip legacy
9 # checks.
10 with h5py.File(test_path, 'w', libver='latest') as test_file:
11 # hdf5 files act as regular numpy arrays after having created
12 # the groups with a simple dictionary syntax.
13 test_file["mesh/nodes"] = nodes
14 test_file["mesh/cable"] = cable
15 test_file["mesh/membr"] = membr
16 test_file["static/bcc"] = np.array(list(bcc))
17 test_file["static/pretension_cable"] = pretension_cable
(continues on next page)
𝐿 = ||𝑥2 − 𝑥1 || (9.3)
4
5 def _side_length(x_a, y_a, z_a, x_b, y_b, z_b):
6 length = np.sqrt(
7 (x_a - x_b) ** 2.0 + (y_a - y_b) ** 2.0 + (z_a - z_b) ** 2.0)
8 return length
Parameters
• len1 (float) – length of the first side of the triangle.
• len2 (float) – length of the second side of the triangle.
• len3 (float) – length of the third side of the triangle.
𝑇i = 𝐿2i 𝑡i 𝑄s 𝑛i (9.14)
Parameters
• length (float) – length of the side opposite to the node whose
nodal force we are computing.
• thick (float) – thickness of the considered membrane element.
• area (float) – area of the considered membrane element resup-
plied after having been computed with the auxiliary function.
• pre (float) – stress density for the investigated element retrieved
in the definition of the problem.
• ni (float) – component of the force in the considered direction;
the user must supply the projection of the normal vector in the x, y or
z-direction.
Returns the tension to be employed in the stress density method.
Return type float
Parameters
9.1.10.1 Cables
Parameters
• x_a (float) – x-coordinate of the first node of the cable.
• y_a (float) – y-coordinate of the first node of the cable.
• z_a (float) – z-coordinate of the first node of the cable.
• x_b (float) – x-coordinate of the second node of the cable.
• y_b (float) – y-coordinate of the second node of the cable.
• z_b (float) – z-coordinate of the second node of the cable.
• prec (ndarray) – numpy array containing the pretension of the
cables as a force in vector form; the index represents the position of
the cable element in the cable array, whereas the vector component
the pretension itself.
9.1.10.2 Membrane
where the index ‘i’ spans all the sides opposite to the considered node.
Parameters
• node_x (int) – x-position of the node in the stacked vector of both
unknowns and equations.
This section is devoted to the verification of the results of the algo deeply described in the
previous Chapter 9. The logical reasoning followed is to initially compute the shape of a plane
cable net and to check that the FDM (Chapter 2) is properly implemented in the the code. After
that a three-dimensional cable net has been tested, the attention was moved to the membrane
structures. In order to verify the correct coding of the SSDM (Chapter 3) a 3D hypar membrane,
constrained at its edges, has been analysed with our algorithm and compared with the results
coming from Dlubal. Once the previous examples have been proved to match the benchmarks,
two other cases of membrane with collaborating cables have been developed and tested.
Tests are run with the following wrapper function that retrieves the desired test case in binary
format in hdf5 and asserts the positiveness of the test.
test_ff_clc._retrieve_and_test_ff(test_path)
Wrapper function to be used in the actual test functions. It executes the form finding
procedure and checks if the prescribed tolerance is met with respect to the test case.
Parameters test_path (str) – path to the case test in hdf5 format.
Returns asserts the positiveness of the test if the latter passed.
Return type None
277
Form-Finding and Mechanical Behaviour of Cable and Membrane Structures
10.1.1 Introduction
The first example is a plane square cable net which is pretensioned with the same force density
for all the cables; clearly, the cables on the borders will have a higher axial force since the
form-finding will privilege their elongation with respect to the central ones.
The test suite includes a test case for a plane square cable net which is also compared graphi-
cally to the same result in Dlubal.
test_ff_clc.test_net()
The form-finding procedure is tested against a simple square cable net with dimension
8m x 8m with 5 cables in each direction. The cable net is hinged at its corners. As for
the material, the cables share a cross-section of 50𝑚𝑚2 and Young modulus of 160GPa.
The coefficient of thermal expansion is set to the real value 1.2 · 10−5 even though it is
just a fictitious value to simulate the inelastic deformation due to the prestressing of the
cables.
The test stresses only the part of algorithm concerning the force density method.
Returns asserts the positiveness of the test.
Return type None
10.1.3 Pre-processing
The cable net is drawn in gmsh and then exported the Abaqus compatible format used by
Calculix.
Code block 10.3: Mesh in gmsh of the cable net before the
form-finding process.
1 // Mesh for a cable net with dimension 8m x 8m with 5 cables in each
2 // direction.
3
4 // Characteristic length that sets the target element size of each point.
5 lc = 100;
6
7 // Nodes geometry defining the connections between the cables.
8 Point(1) = {0, 0, 0, lc};
9 Point(2) = {2, 0, 0, lc};
10 Point(3) = {4, 0, 0, lc};
11 Point(4) = {6, 0, 0, lc};
12 Point(5) = {8, 0, 0, lc};
13 Point(6) = {0, 2, 0, lc};
14 Point(7) = {0, 4, 0, lc};
15 Point(8) = {0, 6, 0, lc};
16 Point(9) = {0, 8, 0, lc};
17 Point(10) = {2, 2, 0, lc};
18 Point(11) = {2, 4, 0, lc};
19 Point(12) = {2, 6, 0, lc};
20 Point(13) = {2, 8, 0, lc};
21 Point(14) = {4, 2, 0, lc};
22 Point(15) = {4, 4, 0, lc};
23 Point(16) = {4, 6, 0, lc};
24 Point(17) = {4, 8, 0, lc};
25 Point(18) = {6, 2, 0, lc};
26 Point(19) = {6, 4, 0, lc};
27 Point(20) = {6, 6, 0, lc};
28 Point(21) = {6, 8, 0, lc};
29 Point(22) = {8, 2, 0, lc};
30 Point(23) = {8, 4, 0, lc};
31 Point(24) = {8, 6, 0, lc};
32 Point(25) = {8, 8, 0, lc};
33
34 // Geometrical lines defining each part of the cables.
35 Line(1) = {1, 6};
36 Line(2) = {6, 7};
37 Line(3) = {7, 8};
38 Line(4) = {8, 9};
39 Line(5) = {9, 13};
40 Line(6) = {13, 12};
41 Line(7) = {12, 8};
42 Line(8) = {12, 11};
43 Line(9) = {11, 7};
44 Line(10) = {6, 10};
45 Line(11) = {1, 2};
46 Line(12) = {2, 3};
47 Line(13) = {3, 14};
48 Line(14) = {14, 10};
(continues on next page)
The mechanical and geometric properties of the cable net are here assigned. No specific over-
rides are employed.
The solution deck is a Calculix file which will act as a placeholder for processing phase of
the form-finding in Python. The thermal loads to simulate the prestressing will be added here
automatically.
Both the pre and post-processing and the finite element solution are automated by a script in
cgx in Calculix.
4 view fill
5 rot z
6 hcpy png
7 sys sleep 1
8 sys mv hcpy_1.png form_finding_cable_net_resu_mesh_fill.png
9
10 view line
11 rot z
12 hcpy png
13 sys sleep 1
14 sys mv hcpy_2.png form_finding_cable_net_resu_mesh_line.png
15
22 # solve the finite element model with ccx after the edits of the
23 # form-finding algo.
24 sys ccx_2.14 form_finding_solve
25
26 # read the results of the analysis and load them in memory.
27 read form_finding_solve.frd
(continues on next page)
53
54 ds -4 e 3
55 plot fv all
56 sys sleep 1
57 hcpy png
58 sys mv hcpy_6.png form_finding_cable_net_resu_uz.png
59
60 ds -3 e 3
61 plot fv all
62 sys sleep 1
63 hcpy png
64 sys mv hcpy_7.png form_finding_cable_net_resu_sxx.png
65
66 ds -2 e 3
67 plot fv all
68 sys sleep 1
69 hcpy png
70 sys mv hcpy_8.png form_finding_cable_net_resu_exx.png
71
72 sys montage form_finding_cable_net_resu_sxx.png form_finding_cable_net_
˓→resu_exx.png -tile 2x1 -geometry +0+0 form_finding_cable_net_resu_se.png
73
74 # switch back to deformed shape.
75 rot z
76 ds -4 e 4
77 plot fv all
Code block 10.9: Modified Calculix mesh with the new co-
ordinates of the cable net after the form-finding.
1 *Heading
2 /home/anava/Documents/Projects/tension-structures/src/python/fsd_calculix/
˓→tests/cable_net/cable_net_gmsh.inp
3 *NODE
4 1, 0, 0, 0
5 1, 0.000000, 0.000000, 0.000000
6 2, 2.212766, 1.531917, 0.000000
7 3, 4.000001, 1.957451, 0.000000
8 4, 5.787236, 1.531917, 0.000000
9 5, 8.000000, 0.000000, 0.000000
10 6, 1.531918, 2.212766, 0.000000
11 7, 1.957453, 4.000001, 0.000000
12 8, 1.531918, 5.787234, 0.000000
13 9, 0.000000, 8.000000, 0.000000
14 10, 2.638303, 2.638304, 0.000000
15 11, 2.808515, 4.000003, 0.000000
16 12, 2.638304, 5.361698, 0.000000
17 13, 2.212765, 6.468072, 0.000000
18 14, 4.000003, 2.808515, 0.000000
19 15, 4.000005, 4.000005, 0.000000
20 16, 4.000003, 5.191478, 0.000000
21 17, 4.000001, 6.042540, 0.000000
22 18, 5.361698, 2.638303, 0.000000
23 19, 5.191481, 4.000003, 0.000000
24 20, 5.361699, 5.361698, 0.000000
25 21, 5.787234, 6.468072, 0.000000
26 22, 6.468074, 2.212766, 0.000000
27 23, 6.042534, 4.000001, 0.000000
28 24, 6.468074, 5.787236, 0.000000
29 25, 8.000000, 8.000000, 0.000000
30 ******* E L E M E N T S *************
31 *ELEMENT, type=T3D2, ELSET=Line1
32 26, 1, 6
33 *ELEMENT, type=T3D2, ELSET=Line2
34 27, 6, 7
35 *ELEMENT, type=T3D2, ELSET=Line3
36 28, 7, 8
37 *ELEMENT, type=T3D2, ELSET=Line4
38 29, 8, 9
39 *ELEMENT, type=T3D2, ELSET=Line5
(continues on next page)
The modified deck reports the thermal loads simulating the prestressing of the cable.
10.1.9 Post-processing
The post processing is also fully automated in Calculix; the contour plots of displacements,
stresses and strains are extracted, merged and printed on paper.
Choosing the force density method in the graphical interface of the software, we show the per-
fect adherence between the algorithm developed in the thesis and the commercial algo proposed
by Dlubal.
10.4: Axial force and axial elongation of the plane cable net.
10.1.11 Conclusions
The force density method proves itself to be very effective for two main reasons:
1. The shape developed through the form-finding process are more visually pleasant with
respect to the pure force method. Inner cables are automatically less tensioned than
cables on the borders without having to assign each pretension in the cable.
2. The solution from a computational point of view is extremely fast, since the iterative
solver needs to solve a preconditioned linear system.
10.2.1 Introduction
This example shows the cabilities of the developed algorithm in the case of the analsyis of a 3D
cable structure of an anticlastic net, which has opposite curvatures at the center points, in partic-
ular curved convexly along a longitudinal plane section and concavely along the perpendicular
section.
The test is started with the wrapper function of pytest and the initial guess of the solution,
prescribed while laying down the mesh in gmsh, is a linear interpolation of the nodes between
the position of the lower and higher supports.
test_ff_clc.test_net_anti()
The form-finding procedure is tested against a mesh for an anticlastic cable net with
dimension 5m x 5m with four cables for each direction, hinged at its edges. As for the
material, the cables share a cross-section of 50𝑚𝑚2 and Young modulus of 160GPa. The
coefficient of thermal expansion is set to the real value 1.2 · 10−5 even though it is just a
fictitious value to simulate the inelastic deformation due to the prestressing of the cables.
The test stresses only the part of algorithm concerning the force density method.
Returns asserts the positiveness of the test.
Return type None
10.2.3 Pre-processing
As usual, the mesh is prepared with a script in gmsh and then exported to Calculix in Abaqus
format. The same file is also exported to MED format for an eventual verification in code_aster.
In the file for the problem definition, the section and coefficient of thermal expansion of the
cables are assigned to simulate the initial pretensioning of the cables. In addition, the boundary
conditions for the pins are prescribed. No particulr overrides, for both the static boundary
conditions and the geometric and mechanical properties of the material, are assigned.
We report here the initial solution deck that act as a placeholder for the edits that will be
operated by the form-finding algorithm.
The solution process, both in terms of preprocessing in gmsh, solution in ccx and post-
processing in cgx is fully automated.
60
61 ds -4 e 3
62 plot fv all
63 sys sleep 1
64 hcpy png
65 sys mv hcpy_6.png form_finding_cable_anticlastica_resu_uz.png
66
67 ds -3 e 3
68 plot fv all
69 sys sleep 1
70 hcpy png
71 sys mv hcpy_7.png form_finding_cable_anticlastica_resu_sxx.png
(continues on next page)
˓→cable_anticlastica_resu_es.png
80
81
82 # switch back to the deformed shape
83 rot -z
84 rot u -45
85 rot c -30
86 view disp
87 ds -4 e 4
88 plot fv all
The initial mesh is modified according to the form-finding algorithm to prepare the input file
for subsequent analysis for possible live loads.
The modified deck in Calculix reports the fictitious temperature loads to be assigned to the
cables to simulate the force density which is at the core of the method.
10.2.9 Post-processing
The post-processing in cgx plots the contour maps for the displacements, stresses and strains,
in both coordinate directions; in this case, no graphs for the same quantities at the pins and in
the mid-point of the span for different values of the external load are plot since the prestressing
is applied instantaneously.
As in the previous example of the plane cable net, we employ again the Dlubal software to
verify the developed algorithm; the same force density is assigned to all cables and the linear
system of the force density method is solved returning the exact results retrieved with the
algorithm discussed in this thesis.
10.9: Axial force and axial elongation of the anticlastic cable net.
10.2.11 Conclusions
The example of the 3D anticlastic cable net shows off the cabalities of the algorithm to tackle
3D problems too. The proposed method, in the case of a structured composed of cables only,
needs to solve only a linear system, which, independently from the initial guess, the Newton-
Krylov algorithm solves iteratively in no more than three function evaluations.
10.3.1 Introduction
The example of the hypar membrane with fixed edges tests the part of the developed algorithm
related to the application of the stress density method. In this case, the actual value of the
assigned stress density does not influence the final result, because, without cables, only the
ratio between the stress densities of the elements modifies the final geometry (different values
of the coefficients for the elements on the borders and center for instance).
The test suite includes a case for a hypar membrane with fixed edges with uniform stress den-
sity. The test function is as usual suitable for testing in pytest.
test_ff_clc.test_hypar()
The form-finding procedure is tested against a simple square membrane, hinged at its
corners, but with a couple of hinges located 2m higher than the other two.
The size of the membrane is 5x5m, and its thickness is 2mm. The cables are laid at the
edges of the membrane. The target force in the cables is set to 30 kN (which will be
adjusted according to the force density method) and the isotropic stress in the membrane
to 3 MPa (it will also be corrected according to the stress density method).
In this example, the combined method of the stress density and force density method is
tested.
Returns asserts the positiveness of the test.
Return type None
10.3.3 Pre-processing
The preprocessing for the generation of the mesh is operated in gmsh, where an approximate
shape of the finaly hypar is drawn, thus allowing the Netwon-Krylov algorithm to converge
faster.
The problem definition file contains only the information about the boundary conditions; the
mechanical and geometric properties, being uniform on the all domain in this particular case,
do not influence the final shape.
The solution deck acts as a placeholder for the following verification in Calculix, which will
show that the shape obtained after the form-finding is in a stress-free state since no pretension-
ing is induce being the cables absent.
18 ** Dirichlet boundary conditions for the cables that are fixed at the
19 ** edges.
20 *BOUNDARY
21 Ncables, 1, 3
22
23 ** Impose a null temperature field on the whole model at T0.
24 *INITIAL CONDITIONS, TYPE=TEMPERATURE
25 Nall, 0.0
26
27 ** Static nonlinear analysis: request fixed time incrementation.
28 ** Ask for nodal displacements and reaction forces, and element
29 ** stresses and strains after having imposed the prestress.
30 *STEP, NLGEOM
31 *STATIC, SOLVER=SPOOLES, DIRECT
32 0.01, 1.0, , 0.05
33 *NODE FILE
34 U, RF
35 *EL FILE
36 S, E
37 *END STEP
38 ** No load is applied in this analysis; it is just a skeleton for
39 ** the imposition of additional loads.
The pre-processing in gmsh and cgx, the solution in ccx and the final post-processing in cgx is
fully automated by a script written in Calculix and a bash script to start the analysis program-
matically.
65 ds -4 e 3
66 plot fv all
67 sys sleep 1
68 hcpy png
69 sys mv hcpy_6.png form_finding_membrane_hypar_fixed_resu_uz.png
70
71 ds -3 e 1
(continues on next page)
The geometry obtained after the form-finding procedure is very close to the one employed as
starting guess from gmsh; this ensures that the convergence of the Newton-Krylov algorithm is
fast.
3 *NODE
4 1, 0, 0, 2000
5 1, 0.000000, 0.000000, 2000.000000
6 2, 5000.000000, 0.000000, 0.000000
7 3, 5000.000000, 5000.000000, 2000.000000
8 4, 0.000000, 5000.000000, 0.000000
9 5, 454.545455, 0.000000, 1818.181818
10 6, 909.090909, 0.000000, 1636.363636
11 7, 1363.636364, 0.000000, 1454.545455
12 8, 1818.181818, 0.000000, 1272.727273
13 9, 2272.727273, 0.000000, 1090.909091
14 10, 2727.272727, 0.000000, 909.090909
15 11, 3181.818182, 0.000000, 727.272727
16 12, 3636.363636, 0.000000, 545.454545
17 13, 4090.909091, 0.000000, 363.636364
18 14, 4545.454545, 0.000000, 181.818182
19 15, 5000.000000, 454.545455, 181.818182
20 16, 5000.000000, 909.090909, 363.636364
21 17, 5000.000000, 1363.636364, 545.454545
22 18, 5000.000000, 1818.181818, 727.272727
23 19, 5000.000000, 2272.727273, 909.090909
24 20, 5000.000000, 2727.272727, 1090.909091
25 21, 5000.000000, 3181.818182, 1272.727273
26 22, 5000.000000, 3636.363636, 1454.545455
27 23, 5000.000000, 4090.909091, 1636.363636
28 24, 5000.000000, 4545.454545, 1818.181818
29 25, 4545.454545, 5000.000000, 1818.181818
(continues on next page)
The solution is modified but only null temperatures are introduced since there are no thermal
loads to be applied because pretensioned cables are totally absent. For this reason, the mem-
brane, independently from the assigned value of the stress density, will be in a stress-free state.
10.3.9 Post-processing
Post processing is fully automated in cgx to produce the contour plots of the displacements in
the three directions, the stresses and strains for the coordinate plane directions and the Mises
stresses and strains.
Being the distorsion of the new mesh found by means of the stress density method, the result
obtained with Dlubal using a target stress value (and not a density) is very similar to the one
produced with the algorithm discussed in this thesis. However, Dlubal shows for the the stresses
10.13: Displacements in the x and y-direction of the hypar membrane with fixed edges.
10.14: Normal stresses in the x and y-direction of the hypar membrane with fixed edges.
10.15: Normal strains in the x and y-direction of the hypar membrane with fixed edges.
10.16: Von-Mises stresses and strains of the hypar membrane with fixed edges.
the target value, even though no loads are applied. On the other hand, the membrane has null
strains everywhere.
10.3.11 Conclusions
The verification of the algorithm for the part concerning the stress density method highlights
that, in particular for cable or membrane only structures, the value of the assigned stress density
is not important, as far as its absolute value is concerned, but only the ratio of the stress densities
of the elements in the structure influences the final shape. In fact, different shapes, even though
belonging to the same class, can be obtained by varying the stress density coefficient for the
center and border elements.
10.17: Check of the displacements and the 𝜎x of the hypar membrane with fixed edges.
10.18: Check of 𝜎y and the 𝜎mises of the hypar membrane with fixed edges.
10.19: Check of 𝜀x and the 𝜀y of the hypar membrane with fixed edges.
10.4.1 Introduction
The example concerning the plane membrane with cables on the borders is the first verification
test for mixed cable and membrane tensile structures. The results show all the potentialities
of the studied method, which combines the force and stress density method to solve hybrid
structures. Not only the ratio between the force or stress density method can be respectively
changed, but also the ratio between the stress and force, giving rise to different shapes due to
the different interaction between the taut cables and membranes.
The test suite includes a case for a plane membrane with cables on the borders which are re-
sponsible for the application of the initial prestressing by means of, computationally speaking,
a simulated thermal load.
test_ff_clc.test_plane()
The form-finding procedure is tested against a simple square plane membrane, hinged at
its corners.
The size of the membrane is 5x5m, and its thickness is 2mm. The cables are laid at the
edges of the membrane. The target force in the cables is set to 60 kN (which will be
adjusted according to the force density method) and the isotropic stress in the membrane
to 5 MPa (it will also be corrected according to the stress density method).
In this example, the mixed method (force density plus stress density method) is tested.
Returns asserts the positiveness of the test.
Return type None
10.4.3 Pre-processing
We start by producing the mesh with a gmsh script which first draws the geometrical entities
(points, lines, and surfaces) and then assigns significant physical entities to be retrieved in
Calculix to distinguish between the inner surface and the border lines representing respectively
the membrane and the cables.
48 // Save the file also to the gmsh ASCII format for inspection.
49 Save "membrane_plane.msh";
50 Mesh.Format = 39; // abaqus format
51 Save "membrane_plane_gmsh.inp";
52 Mesh.Format = 33; // med-hdf5 format
53
(continues on next page)
3 *NODE
4 1, 0, 0, 0
5 2, 5000, 0, 0
6 3, 5000, 5000, 0
7 4, 0, 5000, 0
8 5, 499.99999999971, 0, 0
9 6, 999.99999999888, 0, 0
10 7, 1499.9999999988, 0, 0
11 8, 1999.9999999989, 0, 0
12 9, 2499.9999999969, 0, 0
13 10, 2999.9999999991, 0, 0
14 11, 3500.0000000014, 0, 0
15 12, 4000.0000000037, 0, 0
16 13, 4500.0000000022, 0, 0
17 14, 5000, 499.99999999971, 0
(continues on next page)
The text file containing the problem definition is a merge of the ones employed for first the
force density method and then the stress density method. The arguments are stacked cleverly
to allow the parser to process first the cable elements and then the membrane elems.
The deck resembles the skeleton to be used by the form-finding algorithm to assign the fictitious
thermal loads to simulate numerically the prestressing. An elastic isotropic material for both
the membrane and the steel is used, and nonlinear geometric effects are turned on.
Code block 10.37: Calculix deck for the analysis of the plane
membrane with cables.
1 ** Include all the mesh produced by the form-finding algorithm.
2 *INCLUDE, INPUT=all.msh
3 ** Include also the sets defining the membrane and cables.
4 *INCLUDE, INPUT=memb.nam
5 *INCLUDE, INPUT=cables.nam
6
7 ** Material for the membrane and steel; simple linear elastic isotropic
8 ** model.
9 *MATERIAL, NAME=membr
10 *ELASTIC
11 8e3, 0.4
12 *EXPANSION, ZERO=0.0
13 1.2e-5, 0.0
14 *MATERIAL, NAME=steel
15 *ELASTIC
16 160e3, 0.3
17 * EXPANSION, TYPE=ISO, ZERO=0.0
18 1.2e-5, 0.0
19
20 ** The membrane is 2mm thick.
21 *SHELL SECTION, ELSET=Ememb, MATERIAL=membr
22 2
23
24 ** Simulate the cable as a two-node beam element with negligible flexural
25 ** stifness. Approximated to a rectangular beam with same cross-section.
26 ** First order circular beam elements are not possible in Calculix,
27 ** only second order.
28 *BEAM SECTION, ELSET=Ecables, MATERIAL=steel, SECTION=RECT
29 17.72, 17.72
30
31 ** Impose a null temperature field on the whole model at T0.
32 *INITIAL CONDITIONS, TYPE=TEMPERATURE
33 Nall, 0.0
34
35 ** Dirichlet boundary conditions for the four side pins.
36 *BOUNDARY
37 1, 1, 3
38 2, 1, 3
39 3, 1, 3
40 4, 1, 3
41
42 ** Nonlinear static analysis: request a fixed time incrementation;
43 ** apply the pretensioning to the cables at T0, then ask for
44 ** the nodal displacements and reaction forces, but also element stresses
45 ** and strains.
46 *STEP, NLGEOM
47 *STATIC, SOLVER=SPOOLES, DIRECT
(continues on next page)
The pre-post phase and calculation steps are fully automated by a script written in Calculix
using a sort of bash language. The script in gmsh is meshed, the physical groups are written to
file to be retrieved by ccx and finally the countour plots are drawn.
Code block 10.38: Calculix file for the pre-post phase of the
plane membrane with cables.
1 # change element type because gmsh by default uses plane stress elems.
2 sys sed -i "s/CPS3/S3/g" form_finding_resu.inp
3
4 # read the mesh produced by the form-finding algorithm.
5 read form_finding_resu.inp
6
7 # save an image of the mesh as filled geometry.
8 view fill
9 rot z
10 hcpy png
11 sys sleep 1
12 sys mv hcpy_1.png form_finding_membrane_plane_resu_mesh_fill.png
13
˓→resu_uxy.png
57
58 ds -4 e 3
59 plot fv all
60 sys sleep 1
61 hcpy png
62 sys mv hcpy_6.png form_finding_membrane_plane_resu_uz.png
63
64 ds -3 e 1
65 plot fv all
66 min 0
67 max 10
68 sys sleep 1
69 hcpy png
70 sys mv hcpy_7.png form_finding_membrane_plane_resu_sxx.png
71
72 ds -3 e 2
73 plot fv all
74 min 0
75 max 10
76 sys sleep 1
77 hcpy png
78 sys mv hcpy_8.png form_finding_membrane_plane_resu_syy.png
79
80 sys montage form_finding_membrane_plane_resu_sxx.png form_finding_membrane_
˓→plane_resu_syy.png -tile 2x1 -geometry +0+0 form_finding_membrane_plane_
˓→resu_sxxyy.png
81
82 ds -3 e 7
83 plot fv all
84 min 0
85 max 10
86 sys sleep 1
87 hcpy png
88 sys mv hcpy_9.png form_finding_membrane_plane_resu_sm.png
(continues on next page)
˓→resu_exxyy.png
103
104 ds -2 e 7
105 plot fv all
106 sys sleep 1
107 hcpy png
108 sys mv hcpy_12.png form_finding_membrane_plane_resu_em.png
109
110 sys montage form_finding_membrane_plane_resu_sm.png form_finding_membrane_
˓→plane_resu_em.png -tile 2x1 -geometry +0+0 form_finding_membrane_plane_
˓→resu_sem.png
111
112 # switch back to the deformed shape
113 rot z
114 ds -4 e 4
115 plot fv all
The updated geometry after the form-finding procedure is automatically written to a new mesh
file in Calculix to be later processed by the ccx solver.
The modified deck in Calculix contains the nodal temperatures to simulate the initial prestress-
ing; as already explained in the section concerning the analysis of the code, the pretensioning
is due to the taut cables, which in turn stress the membrane.
7 ** Material for the membrane and steel; simple linear elastic isotropic
8 ** model.
9 *MATERIAL, NAME=membr
10 *ELASTIC
11 8e3, 0.4
12 *EXPANSION, ZERO=0.0
13 1.2e-5, 0.0
14 *MATERIAL, NAME=steel
15 *ELASTIC
16 160e3, 0.3
17 *EXPANSION, TYPE=ISO, ZERO=0.0
18 1.2e-5, 0.0
19
20 ** The membrane is 2mm thick.
21 *SHELL SECTION, ELSET=Ememb, MATERIAL=membr
22 2
23
10.4.9 Post-processing
We report here the results in terms of displacements (magnitude, x and y-direction), stresses and
strains (magnitude, x and y-directions). We also plot the same quantities for some characteristic
points versus the applied load multiplier, in this case resembled by the time incrementation.
Even though showing a similar trend for the displacements, we shall immediately point out that
the solution implemented in Dlubal is different: the stress in the cables is a force density but not
10.23: Displacements in the x and y-direction of the plane membrane with cables.
10.24: Normal stresses in the x and y-direction of the plane membrane with cables.
10.25: Normal strains in the x and y-direction of the plane membrane with cables.
10.26: Von-Mises stresses and strains of the plane membrane with cables.
per finite element whereas for geometrical entity, and the target stress in an actual stress. For
simple plane and 3D structures this does not affect the reasoning, instead for complex shapes
Dlubal introduces the projected algorithm, which returns results similar to the ones discussed
in the thesis.
10.4.11 Conclusions
The collaboration between the cables and the membrane is essential to guarantee the structural
stability of the system before, but especially under, the application of the live loads. As already
seen for the nonlinear behaviour of cable structures, the pretension should not however be
excessive, causing the membrane to enter the plastic field under the application of the additional
loads; this would certainly lead to total loss of the aesthetic properties of the material, but at the
same time of the structural properties, because the membrane would manifest a shape similar
to the one of a sagged cable. In conclusion, compressive zones are absolutely to be avoided,
but also plasticity should never be reached; in addition, viscoelasticity under severe loads may
lead to increased further deformation in time.
10.27: Check of the displacements and the 𝜎x of the plane membrane with cables.
10.28: Check of 𝜎y and the 𝜎mises of the plane membrane with cables.
10.5.1 Introduction
The test suite includes a case for a hypar membrane with cables on the borders. pytest collects
iteratively the tests and executes them (eventually parallely with xdist).
test_ff_clc.test_hypar()
The form-finding procedure is tested against a simple square membrane, hinged at its
corners, but with a couple of hinges located 2m higher than the other two.
The size of the membrane is 5x5m, and its thickness is 2mm. The cables are laid at the
edges of the membrane. The target force in the cables is set to 30 kN (which will be
adjusted according to the force density method) and the isotropic stress in the membrane
to 3 MPa (it will also be corrected according to the stress density method).
In this example, the combined method of the stress density and force density method is
tested.
Returns asserts the positiveness of the test.
Return type None
10.5.3 Pre-processing
The mesh is created as usual with the bottom-up approach already seen in gmsh (top-down with
OpenCASCADE would be counterproductive).
17 // Geometric lines defining the edges of the membrane, which are also
18 // the cable themselves.
19 Line(1) = {1, 2};
20 Line(2) = {2, 3};
21 Line(3) = {3, 4};
22 Line(4) = {4, 1};
23
24 // Line looo enclosing the membrane to define the surface.
25 Line Loop(1) = {1, 2, 3, 4};
26
27 // Geometric surface representing the membrane.
28 Surface(1) = {1};
29
30 // Ask for a transfinite surface if possible.
31 Transfinite Surface {1} = {};
32
33 // Physical lines for the cables and membranes; mark this entities for
34 // export.
35 Physical Line("cables", 200) = {1, 2, 3, 4};
36 Physical Surface("memb", 300) = {1};
37
38 // Mesh first order for faster converence.
39 Mesh.ElementOrder = 1;
40 // Create also the sets for the nodes in Calculix.
41 Mesh.SaveGroupsOfNodes = 1;
42
43 // Mesh only 2D; we do not have any 3D elems.
44 Mesh 2;
45
The definition of the problem is a merge of the files used for the force and stress density method.
The file is sorted according to the type of element to allow the parser to trasverse the file only
once.
The solution deck plays a fundamental role in the analysis of the hypar membrane with cables;
not only it acts as the skeleton for the insertion of the thermal loads simulating the initial
inelastic prestressing, but it also drives the nonlinear Newton solver. Greater attention needs
to be payed to the solver params due to the high nonlinearity induced by the interaction of the
cables and membrane.
Code block 10.48: Calculix deck for the analysis of the hypar
membrane with cables.
1 ** Include all the mesh generated by the form-finding algorithm.
2 *INCLUDE, INPUT=all.msh
3 ** Include also the physical sets defining the cables and the membrane.
4 *INCLUDE, INPUT=memb.nam
5 *INCLUDE, INPUT=cables.nam
6
7 ** Assign a simple linear elastic isotropic material to both the membrane
8 ** and the cables.
9 *MATERIAL, NAME=membr
10 *ELASTIC
11 8e3, 0.4
12 *EXPANSION, ZERO=0.0
13 1.2e-5, 0.0
14 *MATERIAL, NAME=steel
15 *ELASTIC
16 160e3, 0.3
17 * EXPANSION, TYPE=ISO, ZERO=0.0
18 1.2e-5, 0.0
19
20 ** The membrane is 2mm thick; assign the material too.
21 *SHELL SECTION, ELSET=Ememb, MATERIAL=membr
22 2
23
24 ** Equivalent beam section for the cable, sharing the same area.
25 ** Beams with circular section are not possible for first order elems.
26 *BEAM SECTION, ELSET=Ecables, MATERIAL=steel, SECTION=RECT
27 17.72, 17.72
28
The pre-post and solution phases are carefully automated using cgx from Calculix. The param-
eters can be changed in the problem definition file, and the printed results for the report can be
straightforwardly generated by cgx which takes care of entire solution workflow until the final
print of the contour plots.
Code block 10.49: Calculix file for the pre-post phase of the
hypar membrane with cables.
1 # change element type since gmsh uses plane stress elems by default.
2 sys sed -i "s/CPS3/S3/g" form_finding_resu.inp
3
4 # read the mesh generated by the form-finding algorithm.
5 read form_finding_resu.inp
6 view line
7
8 # save an image of the mesh as filled geometry.
9 view fill
10 rot -z
11 rot u -50
12 rot c -10
13 hcpy png
14 sys sleep 1
15 sys mv hcpy_1.png form_finding_membrane_plane_resu_mesh_fill.png
16
17 # save an image of the mesh with wireframe geometry.
18 view line
19 rot -z
20 rot u -50
(continues on next page)
˓→resu_uxy.png
64
65 ds -4 e 3
66 plot fv all
67 sys sleep 1
68 hcpy png
69 sys mv hcpy_6.png form_finding_membrane_hypar_resu_uz.png
70
71 ds -3 e 1
72 plot fv all
73 min 0
74 max 7
(continues on next page)
88
89 ds -3 e 7
90 plot fv all
91 min 0
92 max 7
93 sys sleep 1
94 hcpy png
95 sys mv hcpy_9.png form_finding_membrane_hypar_resu_sm.png
96
97 ds -2 e 1
98 plot fv all
99 sys sleep 1
100 hcpy png
101 sys mv hcpy_10.png form_finding_membrane_hypar_resu_exx.png
102
103 ds -2 e 2
104 plot fv all
105 sys sleep 1
106 hcpy png
107 sys mv hcpy_11.png form_finding_membrane_hypar_resu_eyy.png
108
109 sys montage form_finding_membrane_hypar_resu_exx.png form_finding_membrane_
˓→hypar_resu_eyy.png -tile 2x1 -geometry +0+0 form_finding_membrane_hypar_
˓→resu_exxyy.png
110
111 ds -2 e 7
112 plot fv all
113 sys sleep 1
114 hcpy png
115 sys mv hcpy_12.png form_finding_membrane_hypar_resu_em.png
116
117 sys montage form_finding_membrane_hypar_resu_sm.png form_finding_membrane_
˓→hypar_resu_em.png -tile 2x1 -geometry +0+0 form_finding_membrane_hypar_
˓→resu_sem.png
118
119 # switch back to the deformed shape
120 rot -z
121 rot u -50
122 rot c -10
123 ds -4 e 4
124 plot fv all
The new geometry obtained thanks to the form-finding procedure is, as far as the central part
of the mesh is concerned, very similar to the input one; in fact, the convergence has been fast
thanks to the clever initial guess provided by the geometry drawn in gmsh.
The modified deck is again the exact copy of the input one, except for the fact that the thermal
loads to simulate the prestressing have been added. We here suggest to edit the solution deck
before submitting another analysis for live loads to tweak the convergence parameters in case
the automatic incrementation were not enough to fulfill the desired convergence requirements.
24 ** Equivalent beam section for the cable, sharing the same area.
25 ** Beams with circular section are not possible for first order elems.
26 *BEAM SECTION, ELSET=Ecables, MATERIAL=steel, SECTION=RECT
27 17.72, 17.72
28
10.5.9 Post-processing
We here report the final contour plots automated by cgx for the displacements (magnitude, x
and y-direction), stresses and strains (x and y-dir, Mises).
The shapes produced by the Dlubal commercial software are slightly different from the ones
generated by the algorithm discussed in this thesis, because:
1. The force density for the cables is not assigned as a density per finite element, but per
geometrical entitiy, which means for the entire cable between two anchored pins.
10.33: Displacements in the x and y-direction of the hypar membrane with cables.
10.34: Normal stresses in the x and y-direction of the hypar membrane with cables.
10.35: Normal strains in the x and y-direction of the hypar membrane with cables.
10.36: Von-Mises stresses and strains of the hypar membrane with cables.
2. The stress is imposed as a target fixed stress; the software does not allow the definition
of a stress density on different finite elements.
Especially for 3D structures, this difference is generally negligible, especially when the addi-
tional live loads are considered. However, this aspect may impose a limitation on the possible
obtainable shapes since the geometry must be partitioned in advanced to assign different coef-
ficients.
10.5.11 Conclusions
The example of the hypar membrane with cables is the first real-world application against
which the code is tested. Even though only the topology of the problem is important as far as
the form-finding procedure is concerned, we want to stress the fact that a solid guess of the
solution should be generally given to the solver. This can be easily achieved by OpenCAS-
CADE kernel in gmsh, or in case of more complex geometries, an industrial alternative could
be the employment of CAD packages such as Solid Edge, or even the tools in Shape Studio of
Siemens NX Nastran.
10.37: Check of the displacement and the 𝜎x of the hypar membrane with cables.
10.38: Check of 𝜎y and the 𝜎mises of the hypar membrane with cables.
10.40: Check of 𝜀mises and of the internal action in the cable of the hypar membrane with
cables.
The most important goal of the form-finding is the establishment of a mostly uniform pre-
tensioning in the membrane in order to avoid compressive zones after the application of the
external live loads. In addition, a serious pre-tensioning ensures that the membrane is nowhere
slack, thus the exhibited behaviour, even though geometrically nonlinear, is only slightly su-
perlinear.
We therefore verify an hexagonal membrane pre-tensioned by means of cables on both the
borders and from the pins at ground level to the top of the middle columns.
11.1.1 Pre-processing
The pre-processing for the production of the mesh is done in gmsh, requesting a mesh with a
frontal algorithm to increase the regularity of the triangular discretizazion, which in turn allows
a faster convergence of the NK method.
391
Form-Finding and Mechanical Behaviour of Cable and Membrane Structures
53 Mesh 2;
54
55 Save "membrana_hex.msh";
56 Mesh.Format = 39;
57 Save "membrana_hex_gmsh.inp";
58 Mesh.Format = 33;
59 Save "membrana_hex.med";
60
61 Print "membrane_hex.png";
11.1.2 Solution
Uniform pretensioning is imposed to both the membrane and the cables. The pins are located
on the borders and at the top of the columns. The coefficient of thermal expansion, as already
outlined, is employed to simulate the prestressing of the cables.
Code block 11.6: ASCII file for the definition of the form-
finding problem for the hexagonal membrane.
1 CHARACTERISTICLENGTH
2 1000.0
3 BCC
4 1
5 2
6 3
7 4
8 5
9 6
10 7
11 8
12 PRECABLE
13 60.0e3
14 1 60.0e3
15 CABLEAREA
16 314.0
17 1 314.0
18 CABLEYOUNG
19 160E3
20 1 160E3
21 CABLEEXPANSION
22 1.2E-5
23 1 1.2E-5
24 PREMEMBRANE
(continues on next page)
Code block 11.8: Calculix file for the pre-post and analysis
phase of the hexagonal membrane.
1 # change element type since gmsh by default outputs plane stress
2 # elements.
3 sys sed -i "s/CPS3/S3/g" form_finding_resu.inp
4
5 # read the mesh produced by the form-finding algorithm.
6 read form_finding_resu.inp
7 view line
8
9 # save an image of the filled geometry.
(continues on next page)
74 ds -4 e 3
75 plot fv all
76 sys sleep 1
77 hcpy png
78 sys mv hcpy_6.png membrane_hex_displacement_z.png
79
80 ds -3 e 1
81 plot fv all
82 max 20
83 sys sleep 1
84 hcpy png
85 sys mv hcpy_7.png membrane_hex_stress_xx.png
86
87 ds -3 e 2
88 plot fv all
89 max 15
90 sys sleep 1
91 hcpy png
92 sys mv hcpy_8.png membrane_hex_stress_yy.png
93
94 sys montage membrane_hex_stress_xx.png membrane_hex_stress_yy.png -tile
˓→2x1 -geometry +0+0 membrane_hex_stress_xxyy.png
95
96 ds -3 e 7
97 plot fv all
98 max 10
99 sys sleep 1
100 hcpy png
101 sys mv hcpy_9.png membrane_hex_stress_m.png
102
103 ds -2 e 1
104 plot fv all
105 sys sleep 1
106 hcpy png
107 sys mv hcpy_10.png membrane_hex_strain_xx.png
108
109 ds -2 e 2
110 plot fv all
111 sys sleep 1
112 hcpy png
113 sys mv hcpy_11.png membrane_hex_strain_yy.png
114
115 sys montage membrane_hex_strain_xx.png membrane_hex_strain_yy.png -tile
˓→2x1 -geometry +0+0 membrane_hex_strain_xxyy.png
116
117 ds -2 e 7
118 plot fv all
(continues on next page)
137
138 graph mid time STRESS Mises
139 sys mv graph_3.ps membrane_hex_stress_mises_mid.eps
140 graph mid time TOSTRAIN EXX
141 sys mv graph_4.ps membrane_hex_strain_xx_mid.eps
142 graph mid time TOSTRAIN EYY
143 sys mv graph_5.ps membrane_hex_strain_yy_mid.eps
144
145 sys ps2pdf membrane_hex_strain_xx_mid.eps
146 sys ps2pdf membrane_hex_strain_yy_mid.eps
147 sys pdfnup --no-landscape --a4paper --nup 1x2 --outfile membrane_hex_
˓→strain_xxyy_mid.pdf membrane_hex_strain_xx_mid.pdf membrane_hex_strain_
˓→yy_mid.pdf
148
149 graph mid time TOSTRAIN Mises
150 sys mv graph_6.ps membrane_hex_strain_mises_mid.eps
151
˓→front.pdf
185
186 sys ps2pdf membrane_hex_displacement_mid.eps
187 sys ps2pdf membrane_hex_displacement_front.eps
188 sys pdfnup --no-landscape --a4paper --nup 1x2 --outfile membrane_hex_
˓→displacements.pdf membrane_hex_displacement_mid.pdf membrane_hex_
˓→displacement_front.pdf
189
190 # switch back to the deformed shape
191 rot -z
192 rot u -50
193 rot c -15
194 view disp
195 ds -4 e 4
196 plot fv all
Code block 11.9: Bash script to start Calculix for the solu-
tion of the hexagonal membrane.
1 #!/usr/bin/env bash
2 cgx_2.14 -b membrane-hex.fbd
11.1.3 Post-processing
The post-processing phase is automated with a script in Calculix to produce both the countour
plots and the graphs. Each analysis is fully automated from the pre to the post-processing.
11.2: Total displacements of the hexagonal membrane under the live loads.
11.6: Von-Mises stresses and strains of the hexagonal membrane under the live loads.
-2
-4
-6
-8
-10
D3
-12
-14
-16
-18
-20
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
-0.5
-1
-1.5
-2
-2.5
D3
-3
-3.5
-4
-4.5
-5
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
11.7: Displacements for the front and middle membrane of the hex structure.
1.6
1.4
1.2
1
SXX
0.8
0.6
0.4
0.2
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
3.5
2.5
SYY
1.5
0.5
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
4
SXX
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
3.5
2.5
SYY
1.5
0.5
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
-2x10-5
-4x10-5
-6x10-5
-8x10-5
EXX
-0.0001
-0.00012
-0.00014
-0.00016
-0.00018
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
0.0003
0.00025
0.0002
EYY
0.00015
0.0001
5x10-5
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
0.0006
0.0005
0.0004
EXX
0.0003
0.0002
0.0001
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
0.00014
0.00012
0.0001
8x10-5
EYY
6x10-5
4x10-5
-5
2x10
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
5
Mises
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
0.0009
0.0008
0.0007
0.0006
Mises
0.0005
0.0004
0.0003
0.0002
0.0001
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
11.12: Von-Mises stresses and strains of the front membrane of the hexagonal structure.
4
Mises
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
0.0007
0.0006
0.0005
Mises
0.0004
0.0003
0.0002
0.0001
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Time
11.13: Von-Mises stresses and strains of the middle membrane of the hexagonal structure.
The analysis in Calculix is compared to the same structural case solved in the commercial soft-
ware Dlubal with the use of the standard algorithm (the projected algorithm is different from the
formulation proposed in this thesis). Both the form-finding step conducted with the developed
algorithm and the step for the additional live loads show an extremely good adherence.
11.1.5 Conclusions
A serious form-finding, as shown by the countour plots produced by Calculix, allows the engi-
neer to avoid the presence of zones where the membrane is slack or even compressed. However,
the most notable aspect is that all the studied static quantities, stresses, strains and displace-
ments, vary linearly with the intensity of the load. This aspect is particularly important be-
cause, having run two analysis in an interval, for example for two loads producing a maximum
Von-Mises stress of the 10% and 80% of the yield stress, we could extrapolate the intermediate
values by simple linear interpolation. As for cables, the behaviour of membranes is mostly
linear for taut membranes, whereas strongly not linear for slack elements.
11.17: Check of the form-finding: 𝜀mises and the tension in the cables
11.21: Check after the application of a live load: 𝜀mises and the tension in the cables
The formulations described in Chapter 2 and 3 do not consider the possibility to perform an
actual structural optimization in the case of nonlinear constitutive laws or in presence of a
significant interaction with the soil or substructure.
We therefore introduce a new approach based on the constrained optimization of a convex
function which defines the target for one or more specific static or geometric quantities. The
problem can be mathematically formulated as:
where the inequality constraints are employed to establish either the domain for a specific
variable (an area must be strictly positive) or the interval where the algorithm is likely to find
the target value (limit the search to a sensible interval).
As we will later discuss, the problem can be greatly simplified if we consider 𝑛 optimization
parameters and exactly 𝑛 targets so that, even though the solution of the problem may still not
be unique, the solution can be approached by the solution of a nonlinear system of 𝑛 unknowns
with the Newton-Krylov algorithm.
This approach must be preferred for thourough and precise investigations of the mechanical
behaviour of membranes but limited to problems which do not involve a large number of un-
knowns unless we are aware of the computational cost of the operations.
We have already discussed the example of a plane square membrane with fixed edges (let us
image infinitely rigid cables on the border) loaded perpendicular to its plane. We want now to
compute the external distributed load, which is able to produce a displacement of approximately
93 mm. We already know that the desired load is 1𝐾𝑁/𝑚2 .
427
Form-Finding and Mechanical Behaviour of Cable and Membrane Structures
The deck in code_aster is written in pure Python, but the file is later translated in Fortran
code to speed up the execution. For this reason, in order to pass the parameters to the deck
written in code_aster, we must first serialize an object to file which represents the params of the
optimization, which will be later read by code_aster before the translation to Fortran. Equally,
at each iteration, code_aster will save the value of the target parameters to the same file which
will be retrieved by the main function which is responsible for the guidance of the Netwon-
Krylov method.
class aster_params.AsterParams(*args, **kwargs)
Wrapper class to encapsulate the parameters to be passed to code_aster for the solution
of the optimization problem.
Parameters cannot be passed directly to the solver since the deck of the analysis is moved
to a temporary folder and especially a new deck is created at every iteration. The class,
thanks to the module pickle, takes care of the serialization of the parameters to be passed
to solver which are written on disk and thus accessible by the python deck before the
execution of the C and Fortran routines for the actual solution.
save(dirpath)
Saves the parameters in the pickle file in serialized binary format.
Dumps the parameters in a neutral serialized binary file to be retrieved by code_aster
before the beginning of the analysis. The number of parameters is generally limited,
so pickle is a simple but straightforward choice.
Parameters dirpath (str) – path to save the pickle file containing the
parameters used in the optimization.
Returns None
Return type None
Raises IOError
load(dirpath)
Loads the parameters from the picke file de-serializing the binary format and storing
them in the corresponding class.
The params store in the binary file are de-serialised and mapped to the AsterParams
class and stored in memory for modification and reuse by the solver deck.
Parameters dirpath (str) – path to save the pickle file containing the
parameters used in the optimization.
Returns None
Return type None
Raises FileNotFoundError
add_params(params)
Adds new parameters to the dictionary of the class to be later saved by the dedicated
function.
Parameters are saved in the dictionary but not immediately written to file until the
class is instructed to do so. A dictionary must be passed, so as many parameters as
desired can be inserted at the same time.
Parameters params (dict) – dictionary containing the parameters to
be stored in self.__dict__.
Returns None
Return type None
Code_aster is invoked with a simple wrapper function in Python which opens a new dedicated
shell to execute the deck and eventually supresses the output on the terminal, redirecting it only
to the log file with astk.
membrane_iterative_run.solve_aster()
Passes control to code_aster for the solution at each iteration of the optimization.
The command executes a separate shell to run the analysis in code_aster; the system may
throw an exception similar to FileNotFoundError in case that the hardcoded path for the
.export file is not found.
Returns None
Return type None
Raises FileNotFoundError
Code block 12.2: Code_aster input file for the statical opti-
mization problem.
1 def solve_aster():
2 # Run code_aster in a separate shell session and suppress output.
3 subprocess.call("as_run membrane_iterative.export > /dev/null",
4 shell=True)
The optimization is guided by the following code which references the objective function defin-
ing the target of the optimization: in this case, the maximum displacement must be equal to
93mm.
membrane_iterative_run.optimize_aster(x, apar, st_path)
Objective function for the optimization problem in the FEA analysis software
code_aster.
The unknowns are initially mapped to the corresponding parameters in the AsterParams
class instance. The latter are saved on disk by the serialization operated by pickle. The
problem is then submitted for analysis where the target parameters are updated according
to the results of the latest iteration. Finally, the error function is evaluated to check the
convergence of the procedure.
Parameters
• x – unknowns of the problem for which the optimization problem
must be solved.
• apar (AsterParams) – instance of the class containing the pa-
rameters for the optimization problem.
• st_path (str) – path where the serialized binary pickle file is
stored to retrieve the parameters for the optimization.
Returns equations to be solved by the Newton-Krylov algorithm.
Return type ndarray(float)
The optimization is called by the setup prepared by main, which, due to the simplicity of the
analysis hardcodes the starting guess for the solution:
membrane_iterative_run.main()
Main entry point for the solution of the optimization problem of the maximum displace-
ment of a plane membrane pinned along its edges under distributed loading.
The target of the optimization is the value of the distributed load per unit area which can
produce a displacement of 93𝑚𝑚.
The main function supports also a debug feature, which even if hardcoded, allows the
user to stop the process at the first iteration to check if the code_aster deck has been
successfully written. The core of the function gravitates around the Newton-Krylov al-
gorithm which is responsible for the modification of the parameters at each iteration or
line search.
Returns distributed load on the membrane solving the optimization problem.
Return type float
45 # Stop the timer of the basic profiler and print the elapsed time.
46 end_time = time.time()
47 elapsed_time = end_time - start_time
48 print("Total elapsed time: %.2f sec" % elapsed_time)
49
50 # Print and return the solution to add this case to the test suite.
51 print(solution)
52 return solution.x[0]
The mesh for the plane membrane is created with gmsh as already observed previously:
Code block 12.5: Gmsh script to produce the mesh for the
statical optimization.
1 // Characteristic length for the mesh of the membrane.
2 lcm = 0.25;
3
4 // Half length of the edge of the square membrane.
5 l = 2.5;
6
7 // Point entities defining the pins where the membrane is anchored.
8 Point(1) = {-l, -l, 0, lcm};
9 Point(2) = {l, -l, 0, lcm};
10 Point(3) = {l, l, 0, lcm};
11 Point(4) = {-l, l , 0, lcm};
(continues on next page)
The deck in code_aster is rather complex and thorougly commented in the following algorithm:
Code block 12.6: Code_aster input deck file for the statical
optimization problem.
1 ADIR = "/home/anava/Documents/Projects/tension-structures/src/code_aster" \
2 "/membrane_iterative"
(continues on next page)
12 # Load the params that have been saved by the main function to the
13 # pickle file.
14 aparams = AsterParams()
15 aparams.load(ADIR)
16
17 # -------------------------------------------------------------------------
18
19 # Read the mesh created with gmsh and exported to the compatible binary
20 # MED format.
21 mesh = LIRE_MAILLAGE(
22 FORMAT='MED',
23 UNITE=20,
24 )
25
26 # Define a nodal group to allow the possibility to extract the values
27 # with POST_RELEVE_T.
28 mesh = DEFI_GROUP(
29 reuse=mesh,
30 MAILLAGE=mesh,
31 CREA_GROUP_NO=(_F(
32 GROUP_MA='membrane',
33 NOM='membrane',
34 ),),
35 )
36
37 # -------------------------------------------------------------------------
38
39 # Create the finite element model and assign to the membrane the
40 # 'membrane' finite element. Partition the mesh for parallel solving
41 # with Metis.
42 model = AFFE_MODELE(
43 AFFE=_F(
44 GROUP_MA='membrane',
45 PHENOMENE='MECANIQUE',
46 MODELISATION='MEMBRANE',
47 ),
48 DISTRIBUTION=_F(PARTITIONNEUR='METIS',
49 METHODE='SOUS_DOMAINE',
50 ),
51 MAILLAGE=mesh,
52 )
53
54 # -------------------------------------------------------------------------
55
56 # Create the material for the membrane, assign the Young modulus,
57 # the Poisson ratio, the density and the thermal expansion coefficient.
58 materM = DEFI_MATERIAU(
(continues on next page)
91 # -------------------------------------------------------------------------
92
93 # Define the kinematic boundary conditions for the pin supports.
94 pins = AFFE_CHAR_CINE(
95 MECA_IMPO=_F(
96 DX=0,
97 DY=0,
98 DZ=0,
99 GROUP_MA=(
100 'pins',
101 'cables',
102 ),
103 ),
104 MODELE=model,
105 )
106
107 # -------------------------------------------------------------------------
108
109 # Define the function to gradually apply the load to grant a slightly
110 # faster convergence.
111 ramp = DEFI_FONCTION(
112 NOM_PARA='INST',
113 VALE=(0.0, 0.0, 1E-3, 1E-7, 1E-2, 1E-3, 0.1, 0.1, 1.0, 1.0),
114 INTERPOL=('LIN',),
(continues on next page)
226 # -------------------------------------------------------------------------
(continues on next page)
An example run of the program is here reported to show how fast convergence can be achieved:
12.1.8 Post-processing
12.1.9 Tests
Tests are run with a test function which in turn calls the main method of the current folder and
tests for the positiveness of the test:
test_membrane_iterative.test_membrane_iterative()
Tests that the load to obtain a displacement of 93mm for a square membrane 2mm thick
with sides 5m long is 1kN/m2.
The initial guess for the solution is extremely far from the target value, so that the effec-
tiveness of the Newton-Krylov is highlighted.
Returns asserts the positiveness of the test.
Return type None
6
7 def test_membrane_iterative():
8 optimized_load = main()
9
10 assert optimized_load == pytest.approx(1000.0, 0.02)
Tests are invoked for code_aster by the bash script that follows, since pytest has troubles in
employing the correct python version when invoking code_aster:
In this example, we show the topological optimization of a cable under a concentrated load,
as already studied at the beginning of the thesis. The target displacement to be reached is
365mm and the corresponding target function is the minimization of the error function describ-
ing the discrepancy between the maximum displacement for the current iteration and the target
displacement imposed for the analysis.
The optimization concerns two parameters which are chosen arbitrarly by the Truncated New-
ton Algorithm to satisfy the objective function of the minimization but also the constraints
applied on the parameters: the length of is limited to satisfy the requirements of an hypothetic
span, whereas the cross section of the cable satisfies two inequality constraints which limit the
interval to values with a commercial meaning (the cable cannot be too thin or on the other hand
thick).
Returns None
Return type None
Raises FileNotFoundError
add_params(params)
Adds new parameters to the dictionary of the class to be later saved by the dedicated
function.
Parameters are saved in the dictionary but not immediately written to file until the
class is instructed to do so. A dictionary must be passed, so as many parameters as
desired can be inserted at the same time.
Parameters params (dict) – dictionary containing the parameters to
be stored in self.__dict__.
Returns None
Return type None
5 class AsterParams:
6
7 def __init__(self, *args, **kwargs):
8 self.__dict__.update(kwargs)
9
10 def save(self, dirpath):
11 with open(os.path.join(dirpath, 'aster_params.pkl'), 'wb') as fout:
12 pickle.dump(self, fout, pickle.HIGHEST_PROTOCOL)
13
14 def load(self, dirpath):
15 with open(os.path.join(dirpath, 'aster_params.pkl'), 'rb') as fin:
16 tparams = pickle.load(fin)
17 self.__dict__.update(tparams.__dict__)
18
19 def add_params(self, params):
20 self.__dict__.update(params)
Code_aster is invoked with the python wrapper function which opens two serial shells for the
execution of gmsh to update the mesh at every iteration since the length of the cable has been
modified by the TNC algo, and to launch code_aster for the solution of the current iteration.
cable_iterative_run.solve_aster(apar)
Passes control to code_aster for the solution at each iteration of the optimization.
The command executes a separate shell to run the analysis in code_aster; the system may
throw an exception similar to FileNotFoundError in case that the hardcoded path for the
.export file is not found.
Warning: Pay attention to gmsh! No specific error codes are returned in case that
the wrong parameters are passed to the .geo file. Check mesh at least once before
running the optimization with the debug flag.
The optimization targets a single error function for the maximum displacement of the cable:
cable_iterative_run.optimize_aster(x, apar, st_path)
Objective function for the optimization problem in the FEA analysis software
code_aster.
Being a sort of topology optimization problem, the objective function takes care also of
meshing the new geometry at every step.
The unknowns are initially mapped to the corresponding parameters in the AsterParams
class instance. The latter are saved on disk by the serialization operated by pickle. The
problem is then submitted for analysis where the target parameters are updated according
to the results of the latest iteration. Finally, the error function is evaluated to check the
convergence of the procedure.
Parameters
• x – unknowns of the problem for which the optimization problem
must be solved.
• apar (AsterParams) – instance of the class containing the pa-
rameters for the optimization problem.
• st_path (str) – path where the serialized binary pickle file is
stored to retrieve the parameters for the optimization.
Returns equations to be solved by the Newton-Krylov algorithm.
Return type ndarray(float)
The optimization is called by the setup prepared by main, which, due to the simplicity of the
problem, hardcodes the starting guesses of the solution:
cable_iterative_run.main()
Main entry point for the solution of the optimization problem concerning the combination
of maximum length and minimum cross-section of the cable to sustain a concentrated
load located at midspan under a certain pretension.
The known parameters are the elastic modulus of the cable, the coefficient of thermal
expansion to apply the fictitious elastic thermal load equivalent to the pretension, the
Poisson ratio and the prestressing.
The main function supports also a debug feature, which even if hardcoded, allows the
user to stop the process at the first iteration to check if the code_aster deck has been
successfully written.
The core of the function gravitates around the TNC Truncated Newton Algorithm which
is responsible for the iterative procedure and line search during the minimization of the
objective function.
Returns expansion of two floats representing the cross-section of the cable
and the horizontal distance between the supports.
Return type float
The mesh for the cable is created with gmsh, which accepts the length of the cable as a param-
eter on the command line at every iteration:
Code block 12.14: Gmsh script to produce the mesh for the
topological optimization.
1 // Characteristic length of the cable; for simple 1D elements,
2 // it equals the actual discretization operated by the meshing algorithm.
3 lc = 0.1;
4
5 // Point geometrical entities; the variable 'l' must be either suppòied by
6 // the user on the command line or by the script which is invoking gmsh.
7 Point(1) = {0.0, 0.0, 0.0, lc};
8 Point(2) = {l, 0.0, 0.0, lc};
9 Point(3) = {l/2, 0.0, 0.0, lc};
10
11 // Physical point entities defining the nodes representing the pins and
12 // one for the point of application of the load.
13 Physical Point("pins", 100) = {1, 2};
14 Physical Point("vload", 101) = {3};
15
16 // Geometrical lines for the two parts of the cable, since the latter has
17 // been split to allow the insertion of the physical point representing
18 // the point of application of the load.
19 Line(1) = {1, 3};
20 Line(2) = {3, 2};
21
The deck in code_aster is rather complex, and thorougly commented to show the capabilities
of the solver:
17 # ------------------------------------------------------------------------
18
19 # Read the mesh created with gmsh and exported to the compatible binary
20 # MED format.
21 mesh = LIRE_MAILLAGE(
22 FORMAT='MED',
23 UNITE=20,
24 )
25
26 # Define a nodal group to allow the possibility to extract the values
27 # with POST_RELEVE_T.
28 mesh = DEFI_GROUP(
(continues on next page)
144 # -------------------------------------------------------------------------
145
146 # Define the sectional properties of the cable, among which the
147 # section of the cable and the initial fictitious tension to allow
148 # the calculation of the tangent matrix during the first iteration.
149 cara = AFFE_CARA_ELEM(
150 MODELE=model,
151 CABLE=_F(
152 GROUP_MA='cable',
153 SECTION=aparams.section,
154 N_INIT=1.0,
155 ),
156 )
157
158 # -------------------------------------------------------------------------
159
160 # Define the kinematic boundary conditions for the pin supports.
161 pins = AFFE_CHAR_CINE(
162 MECA_IMPO=_F(
163 DX=0,
164 DY=0,
165 DZ=0,
166 GROUP_MA=(
167 'pins',
168 ),
169 ),
170 MODELE=model,
171 )
172
173 # -------------------------------------------------------------------------
174
175 # Define the function to gradually apply the load to grant a slightly
176 # faster convergence.
177 ramp = DEFI_FONCTION(
178 NOM_PARA='INST',
179 VALE=(0.0, 0.0, 1.0, 1.0),
180 )
181
182 # Concentrated load to be applied in the middle of the span following
183 # the ramp previously defined.
184 grav = AFFE_CHAR_MECA(
185 FORCE_NODALE=_F(
186 GROUP_NO='vload',
187 FZ=aparams.vload,
188 ),
189 MODELE=model,
190 )
191
192 # -------------------------------------------------------------------------
193
194 # Define the time increments at which the problem must be solved.
195 lisnon = DEFI_LIST_REEL(VALE=(0.0, 0.1, 0.25, 0.5, 0.75, 1.0), )
196
(continues on next page)
229 # -------------------------------------------------------------------------
230
231 # Compute the generalised efforts for the cable and also the nodal
232 # and reaction forces.
233 reslin = CALC_CHAMP(
234 RESULTAT=reslin,
235 reuse=reslin,
236 CONTRAINTE=('EFGE_ELNO',),
237 FORCE=(
238 'REAC_NODA',
239 'FORC_NODA',
240 ),
241 )
242
243 # -------------------------------------------------------------------------
244
245 # Write the results in ASCII format to file.
246 IMPR_RESU(
247 FORMAT='RESULTAT',
248 RESU=_F(
249 FORMAT_R='1PE21.3',
250 RESULTAT=reslin,
251 ),
252 UNITE=8,
(continues on next page)
268 # Knowing that the maximum displacement is obtained for the point of
269 # application of the load, extract DX, DY and DZ for that location.
270 tresu = POST_RELEVE_T(
271 ACTION=_F(
272 INTITULE="Ottimizzazione",
273 OPERATION="EXTRACTION",
274 GROUP_NO='vload',
275 RESULTAT=reslin,
276 NOM_CHAM='DEPL',
277 RESULTANTE=('DX', 'DY', 'DZ'),
278 ), )
279
280 # -------------------------------------------------------------------------
281
282 # Convert the extracted values to a list and pass them to the objective
283 # function of the TNC algorithm by saving the binary object to
284 # file with pickle.
285 optiresu = tresu.EXTR_TABLE()
286 aparams.add_params({"frex": optiresu.values()['DZ'][-1]})
287 aparams.save(ADIR)
288
289 # -------------------------------------------------------------------------
290
291 # Clean the allocated resources for the execution and exit cleanly.
292 FIN()
An example run of the program is the following, reporting not only the number of iterations but
also the number of line searches:
12.2.8 Post-processing
12.2.9 Tests
Tests are conducted with the test function operated by pytest, which is invoked in a separate
shell in order to choose the correct python interpreter:
test_membrane_iterative.test_membrane_iterative()
Tests that the load to obtain a displacement of 93mm for a square membrane 2mm thick
with sides 5m long is 1kN/m2.
The initial guess for the solution is extremely far from the target value, so that the effec-
tiveness of the Newton-Krylov is highlighted.
Returns asserts the positiveness of the test.
Return type None
5 def test_cable_iterative():
6 length, section = main()
7
8 assert [length, section] == pytest.approx([9.277, 0.0001], 1e-2)
457
Form-Finding and Mechanical Behaviour of Cable and Membrane Structures
458 Bibliography
Form-Finding and Mechanical Behaviour of Cable and Membrane Structures
[DaBu205c] A.Day, J.Bunce (no date) - “The analysis of hanging roofs”. ARUP Journal ,
30-31;
[GaLe206] S.Gale, W.Lewis - Patterning of tensile fabric structures with a discrete element
modelusing dynamic relaxation;
[GBZ207] P.D.Gosling, B.N.Bridgens, L.Zhang - Adoption of a reliability approach for
membrane structure analysis;
[Lin208] T.Linthout - Form-finding and patterning of fabric structures using shape opti-
mization techniques;
[LoLS209] S.Lopez, G.LaSala - A finite element approach to static and dynamic analysis of
geometrically nonlinear structures;
[KnKe210] D.A.Knoll, D.E.Keyes (2004) - Jacobian-free Newton–Krylov methods: a survey
of approaches and applications
[MaMo211] B.Maurin, R.Motro (1998) - The Surface Stress Density Method as a Form-
Finding Tool for Tensile Membrane – Engineering Structures, Volume 20, Issue
8, pages 712-719;
[MoTo212] E.Moncrieff, B.H.V.Topping - Computer methods for the generation of mem-
brane cutting patterns
[Pri213] D.Priour (2005) - FEM modeling of flexible structures made of cables, bars and
nets;
[PWB214] B.Philippa, R.Wüchnera, K.U.Bletzinger - Advances in the form-finding of
structural membranes;
[Vee215] D.Veenendaal, P.Block - An overview and comparison of structural form finding
methodsfor general networks;
[WuBl216] R. Wüchner, K.U.Bletzinger - Stress-adapted numerical form finding of pre-
stressed surfaces by the updated reference strategy
[XuCo217] W.Xu, T.F. Coleman (2011) - Solving Nonlinear Equations with Newton-Krylov
Method Based on Automatic Differentiation
[Bar301] M.R.Barnes (1977) - Form finding and analysis of tension space structures by
dynamic relaxation - Unpublished Doctoral thesis, City University London
[Bia302] M.Biasielli (2016) - Reti di Funi e Membrane: Ricerca della Forma e Risposta
Strutturale - Politecnico di Milano;
[Col303] D.Colmerales (2016) – Numerical Analysis of cabled-truss structures – Czech
Technical University;
[Hen304] E.Henrysson (2012) - Conceptual Design and Analysis of Membrane Structures;
[AGC401] AGC - www.agcchem.com/blog/europes-largest-infrastructure-project-relies-
fluon-etfe-film/
[ArJ402] Architect’s Journal - www.architectsjournal.co.uk/home/roof-in-place-on-
hopkins-olympic-velodrome/5216064.article
Bibliography 459
Form-Finding and Mechanical Behaviour of Cable and Membrane Structures
[AS403] AS - en.as.com/en/2017/05/17/football/1495044694_565468.html
[EnBr404] Encyclopedia Brittanica - www.britannica.com/technology/membrane-structure
[Hol405] C.Hollander (2007) - Modeling of Thermal Oxidation and Stress Effects -
www.iue.tuwien.ac.at/phd/hollauer/
[Maf406] Maffeis - www.maffeis.it/en/projects/#/
[MaTe407] Maco Technology - www.macotechnology.com/blog/la-progettazione-delle-
tensostrutture/
[Meh408] Mehgies - www.mehgies.com/mta/Wiki-Textile-Architecture/Form-Finding-of-
Textiel-Architecture-.php
[Tri409] TriPyramid - www.tripyramid.com/?q=tension-elements
[Pin410] Pintrest - www.pinterest.com/pin/361062095119373910/
[Wiki411] Wikipedia - en.wikipedia.org/
460 Bibliography
Index
L T
test_catenary() (in module test_catenary), 70
load() (aster_params.AsterParams method),
test_hypar() (in module test_ff_clc), 311, 365
428, 441
test_membrane_iterative() (in module
M test_membrane_iterative), 438,
453
main() (in module cable_iterative_run), 444
test_net() (in module test_ff_clc), 278
main() (in module catenary), 68
test_net_anti() (in module test_ff_clc), 294
main() (in module ff_clc), 248
test_parabola() (in module test_parabola), 78
main() (in module membrane_iterative_run),
test_plane() (in module test_ff_clc), 336
431
main() (in module parabola), 76 W
O write_test() (in module code_ff.parser), 259
write_topology() (in module code_ff.parser),
optimize_aster() (in module ca-
255
ble_iterative_run), 443
461