Structural Mechanics Module Users Guide
Structural Mechanics Module Users Guide
User’s Guide
Structural Mechanics Module User’s Guide
© 1998–2022 COMSOL
Protected by patents listed on www.comsol.com/patents, or see Help>About COMSOL Multiphysics on the
File menu in the COMSOL Desktop for less detailed lists of U.S. Patents that may apply. Patents pending.
This Documentation and the Programs described herein are furnished under the COMSOL Software License
Agreement (www.comsol.com/sla) and may be used or copied only under the terms of the license
agreement.
COMSOL, the COMSOL logo, COMSOL Multiphysics, COMSOL Desktop, COMSOL Compiler,
COMSOL Server, and LiveLink are either registered trademarks or trademarks of COMSOL AB. All other
trademarks are the property of their respective owners, and COMSOL AB and its subsidiaries and products
are not affiliated with, endorsed by, sponsored by, or supported by those trademark owners. For a list of such
trademark owners, see www.comsol.com/trademarks.
Version: COMSOL 6.1
Contact Information
Visit the Contact COMSOL page at www.comsol.com/contact to submit general inquiries
or search for an address and phone number. You can also visit the Worldwide Sales Offices
page at www.comsol.com/contact/offices for address and contact information.
If you need to contact Support, an online request form is located on the COMSOL Access
page at www.comsol.com/support/case. Other useful links include:
Chapter 1: Introduction
Study Types 55
Introduction . . . . . . . . . . . . . . . . . . . . . . . . 55
Stationary Analysis . . . . . . . . . . . . . . . . . . . . . . 55
Eigenfrequency Analysis . . . . . . . . . . . . . . . . . . . . 57
Mode Analysis . . . . . . . . . . . . . . . . . . . . . . . . 62
Time-Domain Analysis . . . . . . . . . . . . . . . . . . . . . 63
Frequency-Domain Analysis . . . . . . . . . . . . . . . . . . . 64
Mode Superposition. . . . . . . . . . . . . . . . . . . . . . 67
Harmonic Perturbation . . . . . . . . . . . . . . . . . . . . 69
Modal Reduced-Order Models. . . . . . . . . . . . . . . . . . 72
Linearized Buckling Analysis. . . . . . . . . . . . . . . . . . . 77
Bolt Pretension Study . . . . . . . . . . . . . . . . . . . . . 78
Random Vibration (PSD) Study . . . . . . . . . . . . . . . . . 78
Response Spectrum Analysis Study . . . . . . . . . . . . . . . . 79
CONTENTS |3
Membrane . . . . . . . . . . . . . . . . . . . . . . . . . 86
Beam . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
Pipe Mechanics . . . . . . . . . . . . . . . . . . . . . . . 87
Truss . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
Wire . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
Selecting Discretization 89
Shape Function Order . . . . . . . . . . . . . . . . . . . . . 89
Lagrange and Serendipity Shape Functions . . . . . . . . . . . . . 90
Choosing Shape Functions in Multiphysics Models . . . . . . . . . . 90
Implicit Shape Function Orders . . . . . . . . . . . . . . . . . 91
4 | CONTENTS
Attachments . . . . . . . . . . . . . . . . . . . . . . . 122
CONTENTS |5
Modeling Magnetostrictive Materials 176
Magnetostriction Coupling . . . . . . . . . . . . . . . . . . 176
Linear vs. Nonlinear Magnetostriction . . . . . . . . . . . . . . 176
Create the Magnetostriction Interface and Define Domains. . . . . . 177
6 | CONTENTS
Solver Settings for Contact Analysis. . . . . . . . . . . . . . . 235
Monitoring the Solution . . . . . . . . . . . . . . . . . . . 239
Dependent Variables in Contact Analysis. . . . . . . . . . . . . 240
Important Contact Variables . . . . . . . . . . . . . . . . . 243
References for Contact Modeling . . . . . . . . . . . . . . . 245
CONTENTS |7
Modeling Pretensioned Bolts 290
8 | CONTENTS
Using Reduced Integration 355
Stresses 409
Defining Stress. . . . . . . . . . . . . . . . . . . . . . . 409
Invariants of the Stress Tensor . . . . . . . . . . . . . . . . 410
Plane Stress . . . . . . . . . . . . . . . . . . . . . . . . 413
CONTENTS |9
Initial Stresses and Strains . . . . . . . . . . . . . . . . . . 414
External Stress. . . . . . . . . . . . . . . . . . . . . . . 415
External Strain . . . . . . . . . . . . . . . . . . . . . . . 417
Axial Symmetry and Stresses . . . . . . . . . . . . . . . . . 418
10 | C O N T E N T S
Linear Buckling . . . . . . . . . . . . . . . . . . . . . . 707
Damping 709
Rayleigh Damping . . . . . . . . . . . . . . . . . . . . . 709
Loss Factor Damping . . . . . . . . . . . . . . . . . . . . 710
Viscous Damping . . . . . . . . . . . . . . . . . . . . . . 711
Maximum Loss Factor . . . . . . . . . . . . . . . . . . . . 712
Wave Attenuation . . . . . . . . . . . . . . . . . . . . . 713
CONTENTS | 11
Bolt Modeling Theory 783
Tightening Torque . . . . . . . . . . . . . . . . . . . . . 783
References 815
12 | C O N T E N T S
Change Cross Section . . . . . . . . . . . . . . . . . . . . 844
Linear Elastic Material . . . . . . . . . . . . . . . . . . . . 845
Nonlinear Elastic Material . . . . . . . . . . . . . . . . . . 853
Elastoplastic Soil Material. . . . . . . . . . . . . . . . . . . 859
Hyperelastic Material . . . . . . . . . . . . . . . . . . . . 866
Shape Memory Alloy . . . . . . . . . . . . . . . . . . . . 878
Phase Transformation Direction . . . . . . . . . . . . . . . . 883
Piezoelectric Material . . . . . . . . . . . . . . . . . . . . 883
Piezomagnetic Material . . . . . . . . . . . . . . . . . . . 887
Magnetostrictive Material. . . . . . . . . . . . . . . . . . . 891
Viscoelasticity . . . . . . . . . . . . . . . . . . . . . . . 891
Mullins Effect . . . . . . . . . . . . . . . . . . . . . . . 898
Plasticity . . . . . . . . . . . . . . . . . . . . . . . . . 899
Set Variables . . . . . . . . . . . . . . . . . . . . . . . 907
Creep . . . . . . . . . . . . . . . . . . . . . . . . . . 909
Additional Creep . . . . . . . . . . . . . . . . . . . . . . 914
Viscoplasticity . . . . . . . . . . . . . . . . . . . . . . . 915
Porous Plasticity . . . . . . . . . . . . . . . . . . . . . . 923
Soil Plasticity . . . . . . . . . . . . . . . . . . . . . . . 929
Concrete . . . . . . . . . . . . . . . . . . . . . . . . . 933
Rocks . . . . . . . . . . . . . . . . . . . . . . . . . . 935
Fiber . . . . . . . . . . . . . . . . . . . . . . . . . . 937
Thermal Expansion (for Fibers) . . . . . . . . . . . . . . . . 939
Thermal Expansion (for Materials) . . . . . . . . . . . . . . . 940
Hygroscopic Swelling . . . . . . . . . . . . . . . . . . . . 943
Intercalation Strain . . . . . . . . . . . . . . . . . . . . . 945
Initial Stress and Strain. . . . . . . . . . . . . . . . . . . . 947
External Stress. . . . . . . . . . . . . . . . . . . . . . . 948
External Strain . . . . . . . . . . . . . . . . . . . . . . . 952
Inelastic Strain Rate . . . . . . . . . . . . . . . . . . . . . 955
Damage . . . . . . . . . . . . . . . . . . . . . . . . . 958
Activation . . . . . . . . . . . . . . . . . . . . . . . . 963
Safety . . . . . . . . . . . . . . . . . . . . . . . . . . 964
Damping . . . . . . . . . . . . . . . . . . . . . . . . . 970
Mechanical Damping . . . . . . . . . . . . . . . . . . . . 976
Coupling Loss . . . . . . . . . . . . . . . . . . . . . . . 978
Dielectric Loss. . . . . . . . . . . . . . . . . . . . . . . 979
Conduction Loss (Time-Harmonic) . . . . . . . . . . . . . . . 981
CONTENTS | 13
External Stress-Strain Relation. . . . . . . . . . . . . . . . . 982
Rigid Material . . . . . . . . . . . . . . . . . . . . . . . 984
Initial Values (Rigid Material) . . . . . . . . . . . . . . . . . 988
Fixed Constraint (Rigid Material) . . . . . . . . . . . . . . . . 989
Prescribed Displacement/Rotation . . . . . . . . . . . . . . . 990
Applied Force (Rigid Material) . . . . . . . . . . . . . . . . . 993
Location Nodes . . . . . . . . . . . . . . . . . . . . . . 995
Applied Moment (Rigid Material) . . . . . . . . . . . . . . . . 995
Mass and Moment of Inertia (Rigid Material) . . . . . . . . . . . 997
Center of Mass Nodes. . . . . . . . . . . . . . . . . . . . 999
Spring Foundation (Rigid Material) . . . . . . . . . . . . . . . 999
Free. . . . . . . . . . . . . . . . . . . . . . . . . . . 1002
Prescribed Displacement . . . . . . . . . . . . . . . . . . . 1003
Prescribed Velocity . . . . . . . . . . . . . . . . . . . . . 1006
Prescribed Acceleration . . . . . . . . . . . . . . . . . . . 1009
Fixed Constraint . . . . . . . . . . . . . . . . . . . . . . 1011
Thermal Expansion (for Constraints) . . . . . . . . . . . . . . 1013
Roller . . . . . . . . . . . . . . . . . . . . . . . . . . 1016
Symmetry . . . . . . . . . . . . . . . . . . . . . . . . 1018
Symmetry Plane . . . . . . . . . . . . . . . . . . . . . . 1021
Antisymmetry . . . . . . . . . . . . . . . . . . . . . . . 1023
Rigid Motion Suppression . . . . . . . . . . . . . . . . . . 1024
Body Load . . . . . . . . . . . . . . . . . . . . . . . . 1026
Gravity . . . . . . . . . . . . . . . . . . . . . . . . . 1028
Base Excitation . . . . . . . . . . . . . . . . . . . . . . 1030
Rotating Frame . . . . . . . . . . . . . . . . . . . . . . 1031
Linearly Accelerated Frame . . . . . . . . . . . . . . . . . . 1034
Boundary Load . . . . . . . . . . . . . . . . . . . . . . 1036
Edge Load . . . . . . . . . . . . . . . . . . . . . . . . 1040
Point Load . . . . . . . . . . . . . . . . . . . . . . . . 1041
Point Load, Free . . . . . . . . . . . . . . . . . . . . . . 1045
Ring Load . . . . . . . . . . . . . . . . . . . . . . . . 1047
Ring Load, Free . . . . . . . . . . . . . . . . . . . . . . 1048
Point Load (on Axis) . . . . . . . . . . . . . . . . . . . . 1050
Thin Layer . . . . . . . . . . . . . . . . . . . . . . . . 1052
Prescribed Displacement (Thin Layer) . . . . . . . . . . . . . . 1053
Fixed Constraint (Thin Layer) . . . . . . . . . . . . . . . . . 1055
Roller (Thin Layer) . . . . . . . . . . . . . . . . . . . . . 1056
14 | C O N T E N T S
Face Load (Thin Layer) . . . . . . . . . . . . . . . . . . . 1057
Boundary Load (Thin Layer) . . . . . . . . . . . . . . . . . 1059
Spring Material. . . . . . . . . . . . . . . . . . . . . . . 1061
Spring Foundation . . . . . . . . . . . . . . . . . . . . . 1062
Thin Elastic Layer. . . . . . . . . . . . . . . . . . . . . . 1068
Predeformation . . . . . . . . . . . . . . . . . . . . . . 1072
Spring-Damper . . . . . . . . . . . . . . . . . . . . . . 1073
Source Point (for Spring-Damper) . . . . . . . . . . . . . . . 1079
Destination Point (for Spring-Damper) . . . . . . . . . . . . . 1079
Source Point Nodes. . . . . . . . . . . . . . . . . . . . . 1079
Destination Point Nodes . . . . . . . . . . . . . . . . . . . 1080
Added Mass . . . . . . . . . . . . . . . . . . . . . . . . 1080
Periodic Condition . . . . . . . . . . . . . . . . . . . . . 1083
Adiabatic Heating. . . . . . . . . . . . . . . . . . . . . . 1086
Cell Periodicity . . . . . . . . . . . . . . . . . . . . . . 1088
Boundary Pair . . . . . . . . . . . . . . . . . . . . . . . 1092
Low-Reflecting Boundary . . . . . . . . . . . . . . . . . . . 1093
Rigid Connector . . . . . . . . . . . . . . . . . . . . . . 1093
Center of Rotation Nodes . . . . . . . . . . . . . . . . . . 1104
Thermal Expansion (Rigid Connector) . . . . . . . . . . . . . . 1104
Applied Force (Rigid Connector) . . . . . . . . . . . . . . . . 1106
Applied Moment (Rigid Connector) . . . . . . . . . . . . . . . 1108
Mass and Moment of Inertia (Rigid Connector) . . . . . . . . . . 1110
Spring Foundation (Rigid Connector) . . . . . . . . . . . . . . 1111
Attachment . . . . . . . . . . . . . . . . . . . . . . . . 1113
Thermal Expansion (Attachment). . . . . . . . . . . . . . . . 1115
Reduced Flexible Components . . . . . . . . . . . . . . . . 1117
Component Definition. . . . . . . . . . . . . . . . . . . . 1119
Fixed Joint . . . . . . . . . . . . . . . . . . . . . . . . 1120
Source Filter . . . . . . . . . . . . . . . . . . . . . . . 1121
Destination Filter. . . . . . . . . . . . . . . . . . . . . . 1121
Contact . . . . . . . . . . . . . . . . . . . . . . . . . 1122
Friction . . . . . . . . . . . . . . . . . . . . . . . . . 1134
Slip Velocity. . . . . . . . . . . . . . . . . . . . . . . . 1138
Adhesion . . . . . . . . . . . . . . . . . . . . . . . . . 1141
Decohesion . . . . . . . . . . . . . . . . . . . . . . . . 1143
Wear . . . . . . . . . . . . . . . . . . . . . . . . . . 1146
Continuity . . . . . . . . . . . . . . . . . . . . . . . . 1149
CONTENTS | 15
Bolt Pretension . . . . . . . . . . . . . . . . . . . . . . 1150
Bolt Selection . . . . . . . . . . . . . . . . . . . . . . . 1153
Bolt Thread Contact . . . . . . . . . . . . . . . . . . . . 1156
Crack . . . . . . . . . . . . . . . . . . . . . . . . . . 1158
Crack Closure . . . . . . . . . . . . . . . . . . . . . . . 1160
Face Load (Crack) . . . . . . . . . . . . . . . . . . . . . 1162
J-Integral . . . . . . . . . . . . . . . . . . . . . . . . . 1163
Reverse Crack Front . . . . . . . . . . . . . . . . . . . . 1165
Port . . . . . . . . . . . . . . . . . . . . . . . . . . . 1166
Elastic Predeformation. . . . . . . . . . . . . . . . . . . . 1170
Phase . . . . . . . . . . . . . . . . . . . . . . . . . . 1172
Harmonic Perturbation . . . . . . . . . . . . . . . . . . . 1174
Stress Linearization . . . . . . . . . . . . . . . . . . . . . 1175
Section Forces . . . . . . . . . . . . . . . . . . . . . . . 1178
Test Material . . . . . . . . . . . . . . . . . . . . . . . 1181
Local System Results . . . . . . . . . . . . . . . . . . . . 1183
Average Rotation. . . . . . . . . . . . . . . . . . . . . . 1185
Wave Speeds . . . . . . . . . . . . . . . . . . . . . . . 1187
Important Variables in the Solid Mechanics Interface . . . . . . . . 1188
16 | C O N T E N T S
The Shell and Plate Interfaces 1242
Domain, Boundary, Edge, Point, and Pair Nodes for the Shell and
Plate Interfaces . . . . . . . . . . . . . . . . . . . . . 1247
Initial Values . . . . . . . . . . . . . . . . . . . . . . . 1249
Thickness and Offset . . . . . . . . . . . . . . . . . . . . 1250
Thickness Change . . . . . . . . . . . . . . . . . . . . . 1252
Linear Elastic Material . . . . . . . . . . . . . . . . . . . . 1253
Layered Linear Elastic Material. . . . . . . . . . . . . . . . . 1255
Shell Local System . . . . . . . . . . . . . . . . . . . . . 1259
Layered Hyperelastic Material . . . . . . . . . . . . . . . . . 1260
Viscoelasticity . . . . . . . . . . . . . . . . . . . . . . . 1262
Plasticity . . . . . . . . . . . . . . . . . . . . . . . . . 1263
Set Variables . . . . . . . . . . . . . . . . . . . . . . . 1265
Creep . . . . . . . . . . . . . . . . . . . . . . . . . . 1265
Additional Creep . . . . . . . . . . . . . . . . . . . . . . 1267
Viscoplasticity . . . . . . . . . . . . . . . . . . . . . . . 1268
Thermal Expansion (for Materials) . . . . . . . . . . . . . . . 1269
Hygroscopic Swelling . . . . . . . . . . . . . . . . . . . . 1273
Initial Stress and Strain. . . . . . . . . . . . . . . . . . . . 1276
External Stress. . . . . . . . . . . . . . . . . . . . . . . 1279
External Strain . . . . . . . . . . . . . . . . . . . . . . . 1281
Inelastic Strain Rate . . . . . . . . . . . . . . . . . . . . . 1283
Damage . . . . . . . . . . . . . . . . . . . . . . . . . 1284
Mullins Effect . . . . . . . . . . . . . . . . . . . . . . . 1285
Safety . . . . . . . . . . . . . . . . . . . . . . . . . . 1286
Damping . . . . . . . . . . . . . . . . . . . . . . . . . 1287
Section Stiffness . . . . . . . . . . . . . . . . . . . . . . 1289
Fixed Constraint . . . . . . . . . . . . . . . . . . . . . . 1291
Prescribed Displacement/Rotation . . . . . . . . . . . . . . . 1293
Thermal Expansion (for Constraints) . . . . . . . . . . . . . . 1296
Prescribed Velocity . . . . . . . . . . . . . . . . . . . . . 1298
Prescribed Acceleration . . . . . . . . . . . . . . . . . . . 1300
Pinned . . . . . . . . . . . . . . . . . . . . . . . . . . 1303
No Rotation . . . . . . . . . . . . . . . . . . . . . . . 1304
Simply Supported. . . . . . . . . . . . . . . . . . . . . . 1306
Symmetry . . . . . . . . . . . . . . . . . . . . . . . . 1307
Symmetry Plane . . . . . . . . . . . . . . . . . . . . . . 1310
Antisymmetry . . . . . . . . . . . . . . . . . . . . . . . 1312
CONTENTS | 17
Body Load . . . . . . . . . . . . . . . . . . . . . . . . 1313
Face Load . . . . . . . . . . . . . . . . . . . . . . . . 1315
Edge Load . . . . . . . . . . . . . . . . . . . . . . . . 1317
Point Load . . . . . . . . . . . . . . . . . . . . . . . . 1320
Point Load, Free . . . . . . . . . . . . . . . . . . . . . . 1322
Ring Load . . . . . . . . . . . . . . . . . . . . . . . . 1325
Ring Load, Free . . . . . . . . . . . . . . . . . . . . . . 1327
Spring Foundation . . . . . . . . . . . . . . . . . . . . . 1329
Predeformation . . . . . . . . . . . . . . . . . . . . . . 1333
Added Mass . . . . . . . . . . . . . . . . . . . . . . . . 1333
Point Mass . . . . . . . . . . . . . . . . . . . . . . . . 1335
Point Mass Damping . . . . . . . . . . . . . . . . . . . . 1336
Periodic Condition . . . . . . . . . . . . . . . . . . . . . 1337
Layered Adiabatic Heating . . . . . . . . . . . . . . . . . . 1341
Boundary to Boundary. . . . . . . . . . . . . . . . . . . . 1344
Edge to Boundary . . . . . . . . . . . . . . . . . . . . . 1345
Edge to Edge . . . . . . . . . . . . . . . . . . . . . . . 1348
Rigid Connector . . . . . . . . . . . . . . . . . . . . . . 1351
Thermal Expansion (Rigid Connector) . . . . . . . . . . . . . . 1359
Attachment . . . . . . . . . . . . . . . . . . . . . . . . 1360
Thermal Expansion (Attachment). . . . . . . . . . . . . . . . 1361
Phase . . . . . . . . . . . . . . . . . . . . . . . . . . 1362
Harmonic Perturbation . . . . . . . . . . . . . . . . . . . 1363
18 | C O N T E N T S
Hyperelastic Material . . . . . . . . . . . . . . . . . . . . 1393
Piezoelectric Material . . . . . . . . . . . . . . . . . . . . 1395
Viscoelasticity . . . . . . . . . . . . . . . . . . . . . . . 1399
Mullins Effect . . . . . . . . . . . . . . . . . . . . . . . 1400
Plasticity . . . . . . . . . . . . . . . . . . . . . . . . . 1401
Set Variables . . . . . . . . . . . . . . . . . . . . . . . 1403
Creep . . . . . . . . . . . . . . . . . . . . . . . . . . 1404
Additional Creep . . . . . . . . . . . . . . . . . . . . . . 1405
Viscoplasticity . . . . . . . . . . . . . . . . . . . . . . . 1407
Thermal Expansion (for Materials) . . . . . . . . . . . . . . . 1408
Hygroscopic Swelling . . . . . . . . . . . . . . . . . . . . 1409
Initial Stress and Strain. . . . . . . . . . . . . . . . . . . . 1411
External Stress. . . . . . . . . . . . . . . . . . . . . . . 1412
External Strain . . . . . . . . . . . . . . . . . . . . . . . 1414
Inelastic Strain Rate . . . . . . . . . . . . . . . . . . . . . 1415
Damage . . . . . . . . . . . . . . . . . . . . . . . . . 1416
Activation . . . . . . . . . . . . . . . . . . . . . . . . 1417
Delamination . . . . . . . . . . . . . . . . . . . . . . . 1419
Safety . . . . . . . . . . . . . . . . . . . . . . . . . . 1422
Damping . . . . . . . . . . . . . . . . . . . . . . . . . 1423
Mechanical Damping . . . . . . . . . . . . . . . . . . . . 1425
Coupling Loss . . . . . . . . . . . . . . . . . . . . . . . 1426
Dielectric Loss. . . . . . . . . . . . . . . . . . . . . . . 1427
Rigid Material . . . . . . . . . . . . . . . . . . . . . . . 1428
Free. . . . . . . . . . . . . . . . . . . . . . . . . . . 1429
Prescribed Displacement . . . . . . . . . . . . . . . . . . . 1430
Prescribed Displacement, Interface . . . . . . . . . . . . . . . 1431
Prescribed Velocity . . . . . . . . . . . . . . . . . . . . . 1432
Prescribed Velocity, Interface . . . . . . . . . . . . . . . . . 1433
Prescribed Acceleration . . . . . . . . . . . . . . . . . . . 1434
Prescribed Acceleration, Interface . . . . . . . . . . . . . . . 1436
Fixed Constraint . . . . . . . . . . . . . . . . . . . . . . 1437
Fixed Constraint, Interface . . . . . . . . . . . . . . . . . . 1438
Thermal Expansion (for Constraints) . . . . . . . . . . . . . . 1439
Roller . . . . . . . . . . . . . . . . . . . . . . . . . . 1440
Roller, Interface . . . . . . . . . . . . . . . . . . . . . . 1441
Symmetry . . . . . . . . . . . . . . . . . . . . . . . . 1442
Antisymmetry . . . . . . . . . . . . . . . . . . . . . . . 1442
CONTENTS | 19
Rigid Motion Suppression . . . . . . . . . . . . . . . . . . 1443
Body Load . . . . . . . . . . . . . . . . . . . . . . . . 1444
Face Load . . . . . . . . . . . . . . . . . . . . . . . . 1446
Rotating Frame . . . . . . . . . . . . . . . . . . . . . . 1447
Linearly Accelerated Frame . . . . . . . . . . . . . . . . . . 1448
Boundary Load . . . . . . . . . . . . . . . . . . . . . . 1449
Edge Load . . . . . . . . . . . . . . . . . . . . . . . . 1451
Line Load. . . . . . . . . . . . . . . . . . . . . . . . . 1452
Point Load . . . . . . . . . . . . . . . . . . . . . . . . 1453
Phase . . . . . . . . . . . . . . . . . . . . . . . . . . 1454
Spring Foundation . . . . . . . . . . . . . . . . . . . . . 1455
Spring Foundation, Interface . . . . . . . . . . . . . . . . . 1456
Thin Elastic Layer. . . . . . . . . . . . . . . . . . . . . . 1458
Thin Elastic Layer, Interface . . . . . . . . . . . . . . . . . . 1459
Predeformation . . . . . . . . . . . . . . . . . . . . . . 1460
Added Mass . . . . . . . . . . . . . . . . . . . . . . . . 1461
Added Mass, Interface . . . . . . . . . . . . . . . . . . . . 1462
Adiabatic Heating. . . . . . . . . . . . . . . . . . . . . . 1463
Rigid Connector . . . . . . . . . . . . . . . . . . . . . . 1466
Rigid Connector, Interface . . . . . . . . . . . . . . . . . . 1467
Thermal Expansion (Rigid Connector) . . . . . . . . . . . . . . 1468
Attachment . . . . . . . . . . . . . . . . . . . . . . . . 1469
Thermal Expansion (Attachment). . . . . . . . . . . . . . . . 1470
Continuity . . . . . . . . . . . . . . . . . . . . . . . . 1471
Chapter 7: Membrane
20 | C O N T E N T S
Plasticity . . . . . . . . . . . . . . . . . . . . . . . . . 1500
Set Variables . . . . . . . . . . . . . . . . . . . . . . . 1501
Creep . . . . . . . . . . . . . . . . . . . . . . . . . . 1502
Additional Creep . . . . . . . . . . . . . . . . . . . . . . 1503
Viscoplasticity . . . . . . . . . . . . . . . . . . . . . . . 1505
Wrinkling . . . . . . . . . . . . . . . . . . . . . . . . 1506
Thermal Expansion (for Materials) . . . . . . . . . . . . . . . 1507
Hygroscopic Swelling . . . . . . . . . . . . . . . . . . . . 1509
Initial Stress and Strain. . . . . . . . . . . . . . . . . . . . 1511
External Stress. . . . . . . . . . . . . . . . . . . . . . . 1513
Inelastic Strain Rate . . . . . . . . . . . . . . . . . . . . . 1515
Safety . . . . . . . . . . . . . . . . . . . . . . . . . . 1516
Damping . . . . . . . . . . . . . . . . . . . . . . . . . 1517
Symmetry . . . . . . . . . . . . . . . . . . . . . . . . 1519
Antisymmetry . . . . . . . . . . . . . . . . . . . . . . . 1521
Layered Adiabatic Heating . . . . . . . . . . . . . . . . . . 1522
Chapter 8: Beam
CONTENTS | 21
Hygroscopic Swelling . . . . . . . . . . . . . . . . . . . . 1546
Initial Load and Strain . . . . . . . . . . . . . . . . . . . . 1547
Implementation . . . . . . . . . . . . . . . . . . . . . . 1548
Theory for Section Stiffness . . . . . . . . . . . . . . . . . . 1550
Stress Evaluation . . . . . . . . . . . . . . . . . . . . . . 1552
Common Cross Sections . . . . . . . . . . . . . . . . . . . 1554
22 | C O N T E N T S
Chapter 9: Beam Cross Section
CONTENTS | 23
Pinned . . . . . . . . . . . . . . . . . . . . . . . . . . 1691
Thermal Expansion (for Constraints) . . . . . . . . . . . . . . 1691
Symmetry . . . . . . . . . . . . . . . . . . . . . . . . 1693
Antisymmetry . . . . . . . . . . . . . . . . . . . . . . . 1695
Edge Load . . . . . . . . . . . . . . . . . . . . . . . . 1696
Point Mass . . . . . . . . . . . . . . . . . . . . . . . . 1698
Point Mass Damping . . . . . . . . . . . . . . . . . . . . 1699
24 | C O N T E N T S
Connection Between Pipes and Structures . . . . . . . . . . . . 1736
Cross Sections. . . . . . . . . . . . . . . . . . . . . . . 1738
CONTENTS | 25
Chapter 13: Multiphysics Interfaces and Couplings
26 | C O N T E N T S
Interface 1805
CONTENTS | 27
Layered Piezoelectric Effect . . . . . . . . . . . . . . . . . . 1847
Electrostriction . . . . . . . . . . . . . . . . . . . . . . 1848
Magnetomechanical Forces . . . . . . . . . . . . . . . . . . 1852
Magnetic Forces . . . . . . . . . . . . . . . . . . . . . . 1854
Piezomagnetic Effect . . . . . . . . . . . . . . . . . . . . 1856
Nonlinear Magnetostriction . . . . . . . . . . . . . . . . . . 1858
Magnetostriction . . . . . . . . . . . . . . . . . . . . . . 1861
Lorentz Coupling 1861
Fluid-Structure Interaction 1863
Fluid-Structure Interaction, Pair . . . . . . . . . . . . . . . . 1865
Fluid-Structure Interaction, Fixed Geometry . . . . . . . . . . . 1868
Fluid-Pipe Interaction . . . . . . . . . . . . . . . . . . . . 1868
Structure Thin-Film Flow Interaction 1870
Shell Thin-Film Flow Interaction . . . . . . . . . . . . . . . . 1871
Index 1903
28 | C O N T E N T S
1
Introduction
This chapter introduces you to the capabilities of this module and includes a
summary of the physics interfaces as well as information about where you can find
additional documentation and model examples. The last section is a brief overview
with links to each chapter in this guide.
29
About the Structural Mechanics
Module
In this section:
The physics interfaces in this module are fully multiphysics enabled, making it possible
to couple them to any other physics interfaces in COMSOL Multiphysics or the other
modules. Available physics interfaces include:
• Solid mechanics for 1D and 2D plane stress, plane strain, and generalized plane
strain, 2D axial symmetry, and 3D solids
• Beams in 2D and 3D, Euler and Timoshenko theory
• Pipes in 2D and 3D
• Truss elements in 2D and 3D
• Wires in 2D and 3D
• Shells and plates, Mindlin theory, 3D and 2D axial symmetry
• Membranes, 3D and 2D axial symmetry
• With the Composite Materials Module, the Layered Shell interface is also available.
30 | CHAPTER 1: INTRODUCTION
The module’s study capabilities include static, eigenfrequency, time dependent
(transient), frequency response, buckling, response spectrum, random vibration, and
parametric studies.
• Linear Elastic Materials can be isotropic, orthotropic, or fully anisotropic, and you
can use local coordinate systems to specify material properties.
• Linear Viscoelastic Materials
• Piezoelectric Materials
• Magnetostrictive Materials are available when used together with the AC/DC
module.
• Material models for hyperelasticity, metal plasticity, porous plasticity, creep,
viscoplasticity, nonlinear elasticity, soil plasticity, concrete, rocks, and clay are
available with the optional Nonlinear Structural Materials Module and
Geomechanics Module.
Coupling structural analysis with thermal analysis is one example of multiphysics easily
implemented with the module, which provides predefined multiphysics couplings for
thermal stress and other types of multiphysics. Piezoelectric materials, coupling the
electric field and strain in both directions are fully supported inside the module
through special multiphysics couplings solving for both the electric potential and the
displacements. Structural mechanics couplings are common in simulations done with
COMSOL Multiphysics and occur in interaction with, for example, fluid flow (fluid–
structure interaction, FSI), chemical reactions, acoustics, electric fields, magnetic
fields, and optical wave propagation.
EIGENFREQUENCY ANALYSIS
An eigenfrequency analysis finds the damped or undamped eigenfrequencies and
mode shapes of a structure, sometimes referred to as the free vibration of a structure.
Prestress effects and damping can also be taken into account.
TRANSIENT ANALYSIS
A transient analysis finds the transient response for a time-dependent model, taking
into account mass and mass moment of inertia. The transient analysis can be either
direct or using a modal solution.
PARAMETRIC ANALYSIS
A parametric analysis finds the solution dependence due to the variation of a specific
parameter, which could be, for instance, a material property or the position of a load.
32 | CHAPTER 1: INTRODUCTION
THERMAL STRESS
In a transient thermal stress study, the program neglects mass effects, assuming that
the time scale in the structural mechanics problem is much smaller than the time scale
in the thermal problem.
LARGE DEFORMATIONS
You can also enable geometric nonlinearity for all structural mechanics interfaces. The
engineering strain is then replaced with the Green–Lagrange strain and the stress with
the second Piola–Kirchhoff stress. To solve the problem, the program uses a total
Lagrangian formulation.
ELASTOPLASTIC MATERIALS
An elastoplastic analysis involves a nonlinear material with or without hardening.
Several isotropic and kinematic hardening models are available.
The elastoplastic material models are available in the Solid Mechanics, Shell, Layered
Shell, Membrane, and Truss interfaces.
HYPERELASTIC MATERIALS
In hyperelastic materials, the stresses are computed from a strain energy density
function. They are often used to model rubber and biological tissue, but are also used
in acoustic elasticity. Many different models are available.
The hyperelastic materials are available in the Solid Mechanics, Shell, Layered Shell,
and Membrane interfaces.
CONTACT MODELING
You can model contact between parts of a structure. The Solid Mechanics, Shell,
Layered Shell, and Membrane interfaces support contact with or without friction.
Three contact algorithms are available: penalty, augmented Lagrangian, and Nitsche
methods. The contact models can be augmented with adhesion and decohesion.
• Bolt pretension
• Bolt thread modeling
• Stress linearization
• Rigid connectors
• Transitions between solid, shell, and beam elements
34 | CHAPTER 1: INTRODUCTION
At any time, a new model can be created or physics interfaces added. Right-click the
Root (top) node and select Add Component or right-click a Component node and select
Add Physics.
Acoustics
Elastic Waves
Fluid-Structure Interaction
36 | CHAPTER 1: INTRODUCTION
PHYSICS INTERFACE ICON TAG SPACE AVAILABLE STUDY TYPE
DIMENSION
Structural Mechanics
38 | CHAPTER 1: INTRODUCTION
PHYSICS INTERFACE ICON TAG SPACE AVAILABLE STUDY TYPE
DIMENSION
40 | CHAPTER 1: INTRODUCTION
PHYSICS INTERFACE ICON TAG SPACE AVAILABLE STUDY TYPE
DIMENSION
42 | CHAPTER 1: INTRODUCTION
PHYSICS INTERFACE ICON TAG SPACE AVAILABLE STUDY TYPE
DIMENSION
5
Requires the addition of the Composite Materials Module.
6 Requires the addition of the CFD Module, or the Polymer Flow, or the Microfluidics
Module.
7
Requires the addition of the Pipe Flow Module.
8
Requires the addition of the Porous Media Flow Module.
9
Requires the addition of the AC/DC Module or the MEMS Module.
10
Requires the addition of the Polymer Flow Module.
44 | CHAPTER 1: INTRODUCTION
displays information about that feature (or click a node in the Model Builder followed
by the Help button ( ). This is called topic-based (or context) help.
• In the Model Builder or Physics Builder, click a node or window and then
press F1.
• On the main toolbar, click the Help ( ) button.
• From the main menu, select Help>Help.
• Press Ctrl+F1.
• From the File menu, select Help>Documentation ( ).
• Press Ctrl+F1.
• On the main toolbar, click the Documentation ( ) button.
• From the main menu, select Help>Documentation.
Once the Application Libraries window is opened, you can search by name or browse
under a module folder name. Click to view a summary of the model or application and
its properties, including options to open it or its associated PDF document.
To include the latest versions of model examples, from the Help menu,
select ( ) Update COMSOL Application Library.
46 | CHAPTER 1: INTRODUCTION
by email. You can also access technical support, software updates, license information,
and other resources by registering for a COMSOL Access account.
48 | CHAPTER 1: INTRODUCTION
THE BEAM CROSS SECTION INTERFACE
The Beam Cross Section chapter describes The Beam Cross Section Interface, which
is used for computing cross section properties for beams. It can also be used for a
detailed evaluation of stresses in a beam when the section forces to which it is subjected
are known. The first section discusses Using the Beam Cross Section Interface, and the
underlying theory is described in Theory for the Beam Cross Section Interface.
• The Thermal Stress, Solid Interface combines a Solid Mechanics interface with a
Heat Transfer interface. The coupling appears on the domain level, where the
temperature from the Heat Transfer interface acts as a thermal load for the Solid
Mechanics interface, causing thermal expansion.
• The Thermal Stress, Shell Interface combines a Shell interface with a Heat Transfer
in Shells interface. The coupling appears on the boundary level, where the
temperature from the Heat Transfer in Shells interface acts as a thermal load for the
Shell interface, causing thermal expansion.
• The Thermal Stress, Layered Shell Interface combines a Layered Shell interface with
a Heat Transfer in Shells interface. The coupling appears on the boundary level,
where the temperature from the Heat Transfer in Shells interface acts as a thermal
load for the Layered Shell interface, causing thermal expansion.
• The Thermal Stress, Membrane Interface combines a Membrane interface with a
Heat Transfer in Shells interface. The coupling appears on the boundary level,
50 | CHAPTER 1: INTRODUCTION
in a situation where the fluid domain can be considered to be nondeforming. The
solid material exists on domains which are adjacent to the fluid.
• The Fluid-Shell Interaction, Fixed Geometry Interface combines fluid flow with the
Shell interface to capture the interaction between the fluid and the solid in a
situation where the fluid domain can be considered to be nondeforming. The shell
is modeled on the boundary of the fluid
• The Fluid-Membrane Interaction, Fixed Geometry Interface combines fluid flow
with the Membrane interface to capture the interaction between the fluid and the
membrane in a situation where the fluid domain can be considered to be
nondeforming. The membrane is modeled on the boundary of the fluid.
• The Fluid-Pipe Interaction, Fixed Geometry Interface combines flow computed
using the Pipe Flow interface with structural analysis in the Pipe Mechanics
interface. Different types of fluid loads are transferred to the structural analysis.
• The Fluid-Solid Interaction, Conjugate Heat Transfer Interface combines fluid flow
with the Solid Mechanics interface and the Heat Transfer in Solids and Fluids
interface. It combines fluid-structure interaction modeling with a nonisothermal
flow. Heat transfer is considered both in the fluid and in the solid in order to capture
thermal expansion effects.
• The Fluid-Solid Interaction, Two-Phase Flow, Phase Field Interface combines
two-phase fluid flow with the Solid Mechanics interface to capture the interaction
between the fluid and the solid in a situation where the fluid domain has significant
deformation. The solid material exists on domains which are adjacent to the fluid.
• The Fluid-Solid Interaction, Two-Phase Flow, Phase Field, Fixed Geometry
Interface combines two-phase fluid flow with the Solid Mechanics interface to
capture the interaction between the fluid and the solid in a situation where the fluid
domain can be considered to be nondeforming. The solid material exists on
domains which are adjacent to the fluid.
The goal of this chapter is to give you an insight on how to approach the modeling
of various structural mechanics problems.
Some physics interfaces and features discussed in this chapter are only available with
certain products. For a detailed overview of the functionality available in each
product, visit https://www.comsol.com/products/specifications/
53
In this chapter:
• Stationary Analysis
• Eigenfrequency Analysis
• Mode Analysis
• Time-Domain Analysis
• Frequency-Domain Analysis
• Mode Superposition
• Harmonic Perturbation
• Modal Reduced-Order Models
• Linearized Buckling Analysis
• Bolt Pretension Study
• Random Vibration (PSD) Study
• Response Spectrum Analysis Study
For general information about study types and solvers, see Studies and
Solvers in the COMSOL Multiphysics Reference Manual
Stationary Analysis
You can consider a structural mechanics problem as stationary if the following two
criteria are fulfilled:
• The loads vary so slowly that inertial forces are negligible. Problems of this type are
referred to as quasi static.
• There are no explicit time dependencies in the material model. Viscoelasticity and
creep have such time dependencies.
STUDY TYPES | 55
In many cases, there is a variation in the load, even though the solution for each value
of the load can be considered as stationary. There are three conceptually different
cases:
• The load values are independent; it is just a number of different load cases you want
to compute. The load case handling functionality described in Load Cases is well
suited for this purpose.
• You want to study a nonlinear problem where the solution is path dependent, or
where the load must be increased in small increments in order to obtain a converged
solution. In this case you should use the parametric continuation solver. Create a
parameter under Global Definitions>Parameters, which you use to control the
variation of the load. Then select Auxiliary sweep under Study Extensions in the
settings for the Stationary solver. In the table for the auxiliary sweep parameters, add
the load controlling parameter, and define its range of variation.
• In a multiphysics problem, another physical quantity might be truly time dependent
but on a time scale that is “slow” from the structural mechanics point of view. This
is usually the case with, for example, problems coupled to heat transfer or diffusion.
If the problem also is such that the structural deformations do not affect the other
physics, it will be unnecessarily expensive to solve also the structural problem in the
time domain, irrespective of whether it is linear or nonlinear. In this situation, you
should first solve the other physics in a time-dependent study and then the structural
mechanics problem in a subsequent stationary study step using the time t as the
parameter in the auxiliary sweep.
CONSTRAINTS
A stationary problem is solvable only if the structure is sufficiently constrained. There
must not be any possible rigid body modes. Thus, no stress-free deformation states are
allowed.
When performing an eigenfrequency analysis, you can specify whether to look at the
mathematically more fundamental eigenvalue, , or the eigenfrequency, f, which is
more commonly used in a structural mechanics context. The relation between the two
is
f = – ---------
2i
2
K – M u = 0
where K is the stiffness matrix, M is the mass matrix, u is the eigenmode displacement
vector, and f is the angular frequency. If damping is present, the eigenvalue
equation is expanded to
2
K + iC – M u = 0 (2-1)
Because only the shape and not the size of the modes (eigenvectors) have physical
significance, the computed modes can be scaled arbitrarily. You can select the method
for scaling in the Eigenvalue Solver node of the solver sequence. If Scaling of eigenvectors
is set to Mass matrix, the eigenmodes u are orthogonalized with respect to the mass
matrix M so that
T
u i Mu i = 1 (2-2)
This is a common choice for the scaling of eigenvectors within the structural mechanics
field. The choice of eigenvector scaling does not affect for example the results of a
subsequent mode superposition analysis, but it will affect the interpretation of an
exported modal representation of the system.
STUDY TYPES | 57
MODAL PARTICIPATION FACTORS
Modal (or mass) participation factors are useful tools when working with the modal
representation of a structure. Through them, you can get the following information:
• The fraction of the total mass of a structure that a certain number of modes
represent is a result. This can be important when judging if a set of modes forms a
good enough base for a mode superposition.
• The main direction of vibration for a certain mode can be seen from the relation
between the participation factors.
• When you have a large set of modes, an examination of the participation factors can
give information about the dominant modes.
To compute modal participation factors, a Participation Factors node must be present
under Definitions in the current component. When you add an Eigenfrequency study
from the Add Study window, such a node is automatically created.
You can also add it manually under Definitions>Physics Utilities. If you do that after an
eigenfrequency study has been run, you need to do an Update Solution in order to get
access to the variables containing the participation factors.
The modal participation factors are available as global variables, and these can for
example be displayed in a table using a Global Evaluation node under Derived Values in
the Results branch. The participation factor results are available as predefined variables
in the Definitions submenu for the component. In Table 2-1, the variables created from
a Participation Factors node is listed (assuming the default tag mpf1).
TABLE 2-1: PARTICIPATION FACTOR VARIABLES
VARIABLE DESCRIPTION
The normalized participation factors are those that would be obtained if mass matrix
scaled eigenmodes would have been used.
If you would compute all eigenmodes of a structure, and sum all modal
masses, they will usually not exactly match the total mass of the structure.
The reason is that any mass which is associated with constrained degrees
of freedom is lost. This effect is discretization dependent. The mass lost
is a fraction of the mass of the elements having constrained nodes.
• Eigenvalue Solver
• Studies and Solvers
• Postprocessing of Eigenmodes
• Derived Values, Evaluation Groups, and Tables
3D 6 (3 translations + 3 rotations)
2D axisymmetric 1 (Z-direction translation)
STUDY TYPES | 59
TABLE 2-2: NUMBER OF POSSIBLE RIGID BODY MODES
In a piezoelectric model, one more zero eigenfrequency could appear if you have not
set a reference value for the electric potential.
In practice, the natural frequencies of the rigid body modes are not computed as
exactly zero, but can appear as small numbers which may even be negative or complex.
If rigid body modes are present in the model, then it is important to use a nonzero
value in the Search for eigenfrequencies around text field in the settings for the
Eigenfrequency study step. The value should reflect the order of magnitude of the first
important nonzero eigenfrequency.
DAMPING
If any type of damping is included in the model, an eigenfrequency solution
automatically returns the damped eigenvalues. The eigenfrequencies and, in general,
also the mode shapes are complex in this case. A complex-valued eigenfrequency can
be interpreted so that the real part represents the actual frequency, and the imaginary
part represents the damping. The damping ratio of the corresponding eigenmode can
be defined as
imag i imag i
i = ------------------------- -------------------------
i real i
where the approximate expression is valid with high accuracy (within 2%) as long as the
damping is less than 0.2.
Some damping types will still give real-valued eigenmodes, this is the case for Rayleigh
damping and loss factor damping.
PRESTRESSED ANALYSIS
In a loaded structure, the natural frequencies may be shifted due to stress stiffening.
The prestress loading can include a contact analysis, in which case the subsequent
eigenfrequency analysis provide as linearization around the current contact state.
Prestressed Structures
2
K + iC – M u = 0 (2-3)
The eigenvalue solver as such assumes that the matrices involved are constant, so they
must be evaluated at a certain frequency, the linearization point.
2
K L + iC L – M u = 0 (2-4)
STUDY TYPES | 61
In order to get a correct solution to Equation 2-3, the linearization point L must be
close to the actual eigenvalue . This is in general possible only for one single
eigenfrequency at a time. You must solve this problem either by manual iteration, or
by using some type of scripting, for example through a model method.
Mode Analysis
The Mode Analysis study type ( ) is available with the Solid Mechanics interface in
2D plane strain.
Elastic waves can propagate over large distances in structures like rails and pipes, with
a generic name referred to as waveguides. After some distance of propagation in a
waveguide of uniform cross section, such guided waves can be described as a sum of
just a few discrete propagating modes, each with its own shape and phase speed. The
equation governing these modes can be obtained as a spatial Fourier transform of the
linearized time-harmonic equation in the waveguide axial z direction or by inserting
the assumption that the mode is harmonic in space,
u = ue – ikz z
Similar to the full time-harmonic equation, the transformed equation can be solved at
a given frequency with a nonzero excitation for most axial wave numbers kz. But at
certain discrete values the equation breaks down. These values are the propagation
constants or wave numbers of the propagating or evanescent waveguide modes. The
The most common use for mode analysis is to define sources for a subsequent
time-harmonic simulation. If there is a component with one or more waveguide
connections, its behavior can be described by simulating its response to the discrete set
of propagating modes on the waveguide opening cross sections.
Time-Domain Analysis
There are two classes of problems where a stationary solution cannot be used:
• When the inertial forces no longer are negligible, the full problem as given by
Newton’s first law must be solved.
• When there are time dependencies in the material model, as for creep or
viscoelasticity.
The most general way of handling time-dependent problems is to use a Time Dependent
study. In this type of analysis, you can incorporate any type of nonlinearity, and there
are no limitations on the time dependence of the loads.
A time-domain solution can be preceded by a stationary study, if, for example, prestress
effects are needed.
For a linear problem including inertia, using the mode superposition method is often
much more efficient than using the standard direct method.
SOLVER SELECTION
The two classes of dynamic problems presented above have quite different properties.
The inertial forces in the full structural dynamics problem contain second-order time
derivatives of the displacements, whereas creep and viscoelasticity only have first-order
STUDY TYPES | 63
derivatives. The physical and numerical properties of these equations differ
significantly.
• Time-Dependent Solver
• Studies and Solvers
Frequency-Domain Analysis
In a frequency-domain analysis, you study the response to a harmonic steady state
excitation for certain frequencies. Such a steady state can prevail once all transient
effects have been damped out.
The response must be linear, so that the single frequency harmonic excitation gives a
pure harmonic response with the same frequency. The model may, however, contain
nonlinearities. The harmonic response is computed around a certain linearization
point. In such a case, the frequency-domain analysis can be considered as a very small
perturbation around that linearization point.
Harmonic Perturbation
• Add a Phase subnode to the load, in which you directly give the phase angle.
• Enter the load as a complex value, for example as
100[N]*(1+0.3*i)/sqrt(1+0.3^2).
Most results from a frequency domain analysis are complex-valued. In many results
evaluation nodes, the real value of any result quantity will be shown. Assuming that
you want to display for example the displacement in the x direction, u, you have the
following options:
• Plot u or real(u). This gives the displacement at current (default zero) phase angle.
• Plot imag(u). This gives the displacement at a phase angle shifted 90 degrees from
the current value.
• Plot abs(u). This gives the amplitude of the displacement.
• Plot arg(u). This gives the phase angle of the displacement.
The reference phase, with respect to which the results above are reported can be
entered in the settings for the dataset.
STUDY TYPES | 65
Result quantities that are nonlinear in terms of the displacements, such as principal
stresses, should be interpreted with great care in frequency domain. They will in
general not be harmonic, so the information about amplitude and phase is not reliable.
PRESTRESSED ANALYSIS
The shift in the natural frequencies in a prestressed structure may have a significant
effect on the frequency response. This is particularly important when the frequencies
of the load are close to any of the natural frequencies of the structure.
• Prestressed Structures
• Harmonic Perturbation
Mode Superposition
Analyzing forced dynamic response for large models can be very time-consuming. You
can often improve the performance dramatically by using the mode superposition
technique. The following requirements must be met for a modal solution to be
possible:
• The analysis is linear. It is possible, however, that the structure has been subjected
to a preceding nonlinear history. The modal response can then be a linear
perturbation around that state.
• There are no nonzero prescribed displacements.
• The important frequency content of the load is limited to a range that is small when
compared to all the eigenfrequencies of the model, so that its response can be
STUDY TYPES | 67
approximated with a small number of eigenmodes. In practice, this excludes wave
and shock type problems.
• If the modal solution is performed in the time domain, all loads must have the same
dependency on the time. This requirement can be relaxed if you use a reduced-order
model to represent the system, rather than using one of the mode superposition
studies.
When using the Structural Mechanics Module, there are three predefined study types
for mode superposition:
The two first of these study types consist of two study steps: One step for computing
the eigenfrequencies and one step for the modal response. The last one has three study
steps. Before the eigenfrequency step, you solve a static load case in order to get a
prestress state used in the eigenfrequency computation.
In practice, you have often computed the eigenfrequencies already, and then want to
use them in a mode superposition. In this case, start by adding an empty study, and
then add a Time Dependent, Modal or Frequency Domain, Modal study step to it. After
having added the study step this way, you must point the modal solver to the solution
containing the eigenfrequencies and eigenmodes. You do this by first selecting Show
default solver at the study level, and then selecting the eigenfrequency solution to be
used in the Eigenpairs section of the generated modal solver.
For many common cases, the mode superposition analysis is not sensitive to whether
the eigenmodes were computed using damping or not. The reason is that the
eigenmodes of problems with Rayleigh damping and loss factor damping can be shown
to be identical to those of the undamped problem, so that the projection to the
subspace spanned by the eigenmodes is the same in both cases. For more general
damping, it is however recommended that you suppress all contributions to the
damping during the eigenfrequency step, and thus base the mode superposition on the
solution to the undamped eigenfrequency problem.
TIME-DEPENDENT ANALYSIS
Only the factor of the load which is independent of time should be specified in the load
features. The dependency on time is specified as Load factor under the Advanced section
of the modal solver. This factor is then applied to all loads.
Harmonic Perturbation
Analyses in the frequency domain assume that the problem your study is linear, at least
with respect to the response to the harmonic excitation. There may be other
nonlinearities, such that the structure has responded nonlinearly to a previous loading.
STUDY TYPES | 69
This loading could, for example, have caused a large rotations or prestress of a rubber
membrane.
The default settings for the different structural mechanics study types in the frequency
domain are summarized in Table 2-3.
TABLE 2-3: DEFAULT PERTURBATION SETTINGS FOR STRUCTURAL MECHANICS STUDY TYPES
• With the default settings you cannot use the same set of loads for a Frequency Domain
and a Frequency Domain, Modal study because only the latter responds to
perturbation loads.
• You can change the behavior of a Frequency Domain study to be of the perturbation
type by modifying the solver sequence. In the General section of the settings for the
Stationary Solver, change Linearity to Linear perturbation.
• A solver that does not have Linearity set to either Linear perturbation or Linear may
respond to nonlinear effects. There are multiphysics problems where this is wanted
because there may be a nonlinearity in another physics, even though the harmonic
solution within structural mechanics is linear. But if there are nonlinearities within
• Pressure loads
• Loads defined in coordinate systems with deformation dependent axis
orientation
• User-defined expressions containing spatial (“lowercase”) coordinates
For most load types, the use of Harmonic Perturbation is straightforward, but some
cases need a more detailed discussion:
• A Rigid Connector can be assigned a Harmonic Perturbation subnode in which you can
prescribe harmonic perturbation values to constrained degrees of freedom. If you
have added Applied Force or Applied Moment nodes under a Rigid Connector, you can
independently assign Harmonic Perturbation to these nodes, so that the loads are
considered as being of the perturbation type.
STUDY TYPES | 71
• Even though initial stresses and strains are not usually considered as loads, you can
assign Harmonic Perturbation also to the Initial Stress and Strain nodes.
• Rather assigning Harmonic Perturbation to a load, you can write the load value
enclosed in the linper() operator. This is particularly useful when the feature that
provide the loading does not have the Harmonic Perturbation option.
This can be employed in two different ways: Either you can use the built-in modal
solvers for the time or frequency domain, or you can export the small equivalent system
and analyze it outside COMSOL Multiphysics, for example, as a component in a larger
system simulation.
ꞏꞏ ꞏ
Mu + Du + Ku = F (2-5)
where u is the displacement vector (size: n-by-1), K is the stiffness matrix (size:
n-by-n), D is the damping matrix (size: n-by-n), and M is the mass matrix (size:
n-by-n). In the frequency domain the problem takes the form
where u = u0eit.
Initially consider the system in the absence of damping and forces. The undamped
system has n eigenvalues i, which satisfy the equation
ˆ = 2 Mu
Ku ˆ (2-6)
i i i
ˆ T Mu
u ˆ = 0 i j , i j (2-7)
j i
ˆ T Ku
u ˆ = 0 i j , i j (2-8)
j i
Next the following n-by-n matrix is constructed, with columns taken from the n
eigenvectors:
ˆ ,u
U = u ˆ u
ˆ
1 2 n
u ˆ u
ˆ T Mu ˆT ˆ
1 1 1 Mu 2
ˆ T Mu
u ˆ u ˆT ˆ
2 1 2 Mu 2
T
U MU =
ˆ T Mu
u ˆ ˆT ˆ
n–1 n – 1 u n – 1 Mu n
ˆ T Mu
u ˆ ˆ T Mu
u ˆ
n n–1 n n
From Equation 2-7 it is clear that this is a diagonal matrix. Similarly, from
Equation 2-8 it is clear that UTKU is also diagonal.
From the properties of the eigenvectors it is possible to expand any function in terms
of the eigenvectors. Thus, the displacement u can be written as:
n
ˆ
u = ai ui
i=1
STUDY TYPES | 73
u = Ua (2-9)
Now consider the original equation: Equation 2-5. First substitute for u using
Equation 2-9. Then transform the equation to the modal coordinate system by
premultiplying by UT. This gives:
T ꞏꞏ T ꞏ T T
U MUa + U DUa + U KUa = U F (2-10)
It has already been established that the matrices UTMU and UTKU are diagonal and
frequently a damping model is chosen that results in a diagonal damping matrix. For
example, in Rayleigh damping D M + K, where and are constants. For a
general damping, the transformed damping matrix is however not diagonal. As an
alternative, a damping ratio, i, can be assigned to each mode.
EIGENVALUE SCALING
The precise form of Equation 2-10 is determined by the normalization adopted for the
eigenfunctions. In structural applications the eigenfunctions are often normalized such
that UTMU I. This is referred to as mass matrix scaling in the eigenvalue solver. In
this case Equation 2-6 gives
ˆ T Ku
u ˆ = 2u
ˆ T Mu
ˆ = 2
i i i i i i
so that
T 2
U KU = diag i
where diag(i2) is the diagonal matrix with diagonal elements i2. Similarly, if
damping ratios for each mode are defined, the damping matrix can be expressed in the
form
T
U DU = diag 2 i i
Thus, if mass matrix scaling is used Equation 2-10 takes the form
ꞏꞏ ꞏ 2 T
a + diag 2 i i a + diag i a = U F (2-11)
It is also possible to scale the eigenvectors so that the point of maximum displacement
has given displacement. This is referred to as max scaling in the eigenvalue solver. For
ꞏꞏ ꞏ T
diag m eff , i a + diag c eff , i a + diag k eff , i a = U F (2-12)
Here meff,i is the effective mass of the i:th mode, ceff,i = 2meff,iii is the effective
damping parameter for the mode, and keff,i is the effective spring constant. Each
element of the vector UTF gives the force component that acts on each of the
respective modes.
REDUCED-ORDER MODELS
The preceding discussion did not consider how to reduce the number of degrees of
freedom in the system. For systems in which the vector UTF has only a few significant
components (for example, components i = 1, …, m where m « n) the following
approximation can be made:
m
ˆ
u ai ui
i=1
u = U'a'
where U' is now an m-by-n and a' is a vector with m components. The equation
system in modal coordinates now takes the form
T ꞏꞏ T ꞏ T T
U' MU'a' + U' DU'a' + U' KU'a' = U' F (2-13)
STUDY TYPES | 75
The matrices U'TMU', U'TDU', and U'TKU' now have dimensions m-by-m.
Similarly, the vector U'TF has m components. This results in a significant reduction in
the system complexity.
2 ˆ + i Du
–i M u ˆ + Kuˆ
r ,i i r ,i r ,i = F
u = U' r a'
where U'r is the n-by-m matrix containing the right eigenvectors chosen for the modal
analysis. Once again a' is a vector with m components. The system in modal
coordinates takes the form
T ꞏꞏ T ꞏ T T
U' l MU' r a' + U' l DU' r a' + U' l KU' r a' = U' l F
where U'l is the n-by-m matrix containing the left eigenvectors chosen for the modal
analysis.
The matrices U'lTMU'r, U'lTDU'r, and U'lTKU'r are no longer necessarily diagonal.
The modal solver accepts any linearly independent set of vectors to project the solution
vector and equations onto and constructs the reduced-order system accordingly.
After the model has solved, right-click the Results>Derived Values node and select
System Matrices. In the output section choose the Matrix to display in the list. The mass
matrix corresponds to the matrix U'lTMU'r the stiffness matrix corresponds to
U'lTKU'r, and the damping matrix corresponds to U'lTDU'r. The vector U'lTF is
available as the load vector. These matrices are given in a format that respects the
COMSOL reports a critical load factor, ,which is the multiplier to the initial load at
which the structure becomes unstable. The corresponding eigenmode is the shape of
the structure in its buckled state.
The level of the initial load used is immaterial since a linear problem is solved. If the
initial load actually was larger than the buckling load, then the critical value of is
smaller than 1. It is also possible that the computed value of is negative. This signifies
that a reversed load will give the critical case.
The buckling computed buckling modes can be used to provide an initial imperfection
for a subsequent nonlinear buckling analysis,
• For more details about how to model buckling, see Buckling Analysis.
• The numerical formulation is described in the section Linear Buckling
in the theory chapter.
• Settings for the solvers are described in Studies and Solvers, Linear
Buckling, and Buckling Imperfection in the COMSOL Multiphysics
Reference Manual.
STUDY TYPES | 77
• Bracket — Linear Buckling Analysis: Application Library path
Structural_Mechanics_Module/Tutorials/bracket_linear_buckling
• Buckling Analysis of a Truss Tower: Application Library path
Structural_Mechanics_Module/Buckling_and_Wrinkling/
truss_tower_buckling
The input to a random vibration analysis is given in terms of power spectral densities
(PSD) and, in the case of several loading sources, the load cross-correlations.
The analysis is based on a mode superposition and the reduced-order model (ROM)
functionality. Except from the computation of eigenfrequencies and corresponding
eigenmodes, and the creation of the reduced model, the core of the computation is
performed during result evaluation.
The Random Vibration (PSD) study is mainly an entry point when adding studies. When
you select it, you actually get two studies and a number of nodes under Global
Definitions added to the model.
The Response Spectrum Analysis study is mainly an entry point when adding studies.
What you actually get when you add such study is an Eigenfrequency study step together
with a Response Spectrum node under Definitions.
The actual response is computed on demand during result evaluation, using the
computed eigenfrequencies and modes. The settings for the evaluation are done in the
Response Spectrum 2D and Response Spectrum 3D datasets.
STUDY TYPES | 79
If your response spectrum evaluation requires inclusion of missing mass correction,
you need also to compute a set of stationary load cases. To set up that analysis, use the
Create missing mass correction study button ( ) on the header of the Response
Spectrum section the in the Response Spectrum node settings.
Solid Mechanics
The Solid Mechanics interface offers the most general modeling of structural
mechanics problems and is based on general principles of continuum mechanics. It is
the interface which contains the largest number of material models, and the most
advanced boundary conditions.
3D SOLID GEOMETRY
The degrees of freedom (dependent variables) in 3D are the global displacements u, v,
and w in the global x, y, and z directions, respectively.
Figure 2-1: Loads and constraints applied to a 3D solid using the Solid Mechanics
interface.
2D GEOMETRY
Plane Stress
The plane stress variant of the 2D physics interface is useful for analyzing thin in-plane
loaded plates. For a state of plane stress, the out-of-plane components of the stress
tensor are zero.
Figure 2-2: Plane stress is used to model plates where the loads are only in the plane; it does
not include any out-of-plane stress components.
The 2D physics interface for plane stress allows loads in the x and y directions, and
assumes that these are constant throughout the material’s thickness, which can vary
with x and y. The plane stress condition prevails in a thin (compared to the in-plane
Plane Strain
The plane strain variant of the 2D physics interface that assumes that all out-of-plane
strain components of the total strain tensor z, yz, and xz are zero.
Loads in the x and y directions are allowed. The loads are assumed to be constant
throughout the thickness of the material, but the thickness can vary with x and y.
Formally, the plane strain condition requires that the ends of the object are constrained
in the z direction, but it is often also used for central parts of an object that is long in
the z direction (compared to the in-plane dimensions). One example is a long tunnel
along the z-axis where it is sufficient to study a unit-depth slice in the xy-plane.
In the default version of the interface, displacements occur only in the r-z plane, and
there are two degrees of freedom, u and w. By selecting the Include circumferential
displacement option, you can model also torsion around the axis of rotational
symmetry. The azimuthal rotation degree of freedom v is then included. In addition,
many features, such as load features, allow values to be entered in the direction.
In general, the interface is suited for modeling the propagation of waves over large
distances relative to the wavelength, for example, ultrasound propagation for
nondestructive testing (NDT), or seismic waves. The interface exists in 2D
(generalized plane strain) and 3D.
The Plate interface is a specialization of the Shell interface, used for 2D modeling in
the XY-plane. A plate model has its main action in bending out of the plane, but can
optionally also treat in-plane forces. If the loads act only in the plane, using Solid
Mechanics with the Plane Stress option is a better choice.
Shells are modeled on boundaries, and the transverse direction is represented only by
the mathematical model. The degrees of freedom consist of displacements and
rotations at the modeled boundary. This results in an assumption where the in-plane
strains vary linearly through the thickness, and the stress in the thickness direction is
zero. The thickness of a shell does not have to be constant, although this is by far the
most common case.
For nonlinear material models, a layered approach with a single layer is used. There is
a virtual mesh in the thickness direction, in order to accommodate the potential
variation of the material properties in the thickness direction.
Rather than computing the shell stiffness from material properties and thickness, you
can also directly enter that stiffness properties in tension, bending, and shear.
The Shell and Plate interfaces can be used both for “thin” and “thick” shells. Shear
deformations are taken into account; this is usually called Mindlin theory. The material
model is linear elastic.
The in-plane stiffness of an elastic shell is proportional to the thickness h, while the
bending stiffness is proportional to h3. The difference in stiffness along different
directions can thus become very large. When an object is very thin, a shell model may
be numerically ill-posed due to the negligible bending stiffness. It is then better to use
the Membrane interface.
Membrane
The Membrane interface can be used for very thin objects, like cloth, where only
in-plane forces are important. Membranes can be considered as plane stress elements
but with an arbitrary, possibly curved configuration in space. The Membrane interface
is available in 3D and 2D axisymmetry.
You can study configurations when there is local wrinkling in a membrane by adding
a special nonlinear material model.
In the Membrane interface a large number of different material models can be used.
When the Composite Materials Module is available, it is possible to model also
multilayered membranes.
Beam
A beam is an abstract model where only the extension in the axial direction is modeled
explicitly on an edge. The cross section is usually specified in terms of geometrical
properties such as area and moments of inertia. Several predefined cross-section types
are also available. The Beam interface is available in 3D and 2D.
The exact stress distribution in the beam is not explicitly modeled. It is actually not
even fully determined by the cross-sectional properties. Instead, six (in 3D) resultant
section forces are used: axial force, shear forces in two perpendicular directions, two
bending moments, and one twisting moment.
• The classical Euler-Bernoulli beam theory, which is applicable for slender beams.
• Timoshenko theory, where also shear deformations are considered. This theory
makes it possible to use the Beam interface to model rather thick beams.
Pipe Mechanics
Pipes are similar to beams, and many properties of the Pipe Mechanics interface are
shared with the Beam interface. The most distinguishing feature is that the internal
pressure usually causes a significant part of the stresses in a pipe. Also, temperature
gradients usually occur through the pipe wall, rather than across the entire section.
The loads from internal pressure and drag forces can be taken directly from results in
the Pipe Flow interface. Similarly, the temperature in the pipe walls can be taken from
the Heat Transfer in Pipes interface.
Truss
The Truss interface has two main purposes:
For a truss model, only one geometrical property is needed, the cross-section area. The
material model can be linear elastic, elastoplastic, or a shape memory alloy. There is also
a special material model for creating spring/damper data.
The truss element has no stiffness in the directions perpendicular to its extension. For
trusses, this is usually not a problem since they are designed such that each member is
stabilized by its neighbors.
The main difference between the Truss interface and the Wire interface is that the
wires cannot sustain any compressive forces. In reality the wire will wrinkle.
Numerically, this is handled by using a very low stiffness in compression.
It is well known that using first-order shape functions in solid mechanics will give an
overly stiff solution, unless a very fine mesh is used. This is especially noticeable for
triangular and tetrahedral elements. This can, to some extent, be counteracted by using
reduced integration, see Using Reduced Integration.
If the purpose of the analysis is only to compute stiffness, rather than stresses, the use
of linear shape functions can still be justified. This is the default choice in the
Multibody Dynamics interface, available with the Multibody Dynamics Module.
If the solution contains discontinuities, for example when some type of front is moving
through the material, first-order elements and a fine mesh is often a good choice, since
the advantage of the higher-order elements lies in their ability to represent smooth
gradients.
TRUSS ELEMENTS
In the Truss interface the default is to use first-order shape functions, since the
elements are mainly used in a context where the axial force in each element is constant.
When truss elements share an edge with other structural elements, you should choose
the same discretization in both interfaces, usually quadratic.
BEAM ELEMENTS
The beam elements have only one set of shape functions, which cannot be changed.
The axial displacement and the twist are represented by first-order shape functions,
while the bending is represented by cubic Hermitian shape functions. This element can
then represent a constant axial force, a constant twisting torque, a linear bending
moment, and a constant shear force. This is the exact solution for a beam having no
distributed loads.
SELECTING DISCRETIZATION | 89
A consequence of this formulation is that it may not possible to obtain a perfectly
conforming approximation if a beam shares an edge with elements from another
physics interface.
The serendipity elements have the advantage of generating significantly fewer degrees
of freedom for structured meshes. The accuracy is in most cases almost as good as for
the Lagrange elements. The Lagrange elements are however less sensitive to strong
mesh distortions.
The serendipity shape functions differ from the Lagrange shape functions only for the
following element shapes:
When coupling two structural mechanics physics interfaces, the same type of shape
functions should be used in both interfaces to ensure conformity in displacement shape
functions. Since there is no difference between the two families of shape functions in
1D, this is not an issue when connecting edges.
Another type of coupling appears on the boundary between two domains having
different physics, as in fluid-structure interaction and acoustic-structure interaction.
When, for example, Thermoviscous Acoustics is coupled to Solid Mechanics, then the
time derivative of the displacement in the solid is set equal to the velocity in the
acoustic medium on the shared boundary. In this case, it makes sense to have the same
shape function order for these two fields.
The discontinuous Lagrange shape functions will have an order that is one below what
is used for the displacement shape functions.
If Gauss point data is used, the same integration points as used for the numerical
integration of the stiffness matrix are used. This order depends on the selected
displacement discretization order and whether reduced integration is used or not.
INELASTIC STRAINS
For material models like plasticity and creep, the inelastic strains are formally degrees
of freedom. They will be allocated at the same integration points as used for the
SELECTING DISCRETIZATION | 91
numerical integration of the stiffness matrix. This order depends on the selected
displacement discretization order and whether reduced integration is used or not.
• Structures that are thin in large regions, but more three-dimensional at certain
locations. A mixture of solids and shells can then significantly reduce the model size.
• Plates or shells having beams as stiffeners.
• Truss elements acting as reinforcement bars in a concrete structure.
• A thin layer of one material on top of another material. In this case, an idealization
with shells or membranes covering the boundary of a solid can be useful.
When several physics interfaces are added in COMSOL Multiphysics, the default is
always that each physics interface has individual and unique degrees of freedom. In
structural mechanics, the first physics interface has the displacement variables (u, v, w),
then the second physics interface has (u2, v2, w2), and so on. This means that the
physics interfaces initially are independent, even when defined on the same geometrical
part. To get the intended interaction requires that a coupling is established between
the physics interfaces.
Various methods to couple different element types are discussed in this section.
Coupling Techniques
The following basic techniques to connect physics interfaces with displacement
degrees of freedom are discussed in this section:
In the Beam, Pipe Mechanics, Shell, and Plate interfaces, the deformation is described
also by rotational degrees of freedom. In the general case, these degrees of freedom
interact with the translational degrees of freedom in a connection.
In some special cases — for example, when a thin shell acts as cladding on a solid — it
is sufficient to make the degree of freedom names for the displacements common; the
rotational degrees of freedom are not important. If, however, a shell edge is connected
to a solid, it acts as a “hinge”, which in most cases is not the intended behavior. You
then need to use the more sophisticated techniques described next.
The default shape functions in the Solid Mechanics interface are of the
serendipity type, whereas Lagrange shape functions are used in the Shell
interface. If you are placing a shell element on the boundary of a solid
element, you must select Lagrange shape functions also in the Solid
Mechanics interface so that the two physics interfaces share the same node
points.
The shape functions used in the Beam and Pipe Mechanics interfaces have
special properties, and a beam cannot have the same degrees of freedom
as another physics interface if the same edge is shared.
Also, the representation of rotations differs between the Shell and Plate
interfaces (displacement of normal) and the Beam and Pipe Mechanics
interfaces (rotation angle). It is therefore not possible to use common
degree of freedom names for the rotational degrees of freedom.
You can choose between two different formulations, by setting Method to either Rigid
or Flexible. The flexible version is significantly more accurate locally at the connected
solid boundary, but it comes with a cost in terms of some extra degrees of freedom.
Also, this method requires a large enough number of degrees of freedom in the
thickness direction of the solid. For second-order elements, typically three elements are
required.
You can choose between two different formulations, by setting Method to either Rigid
or Flexible. The flexible version is significantly more accurate locally at the connected
solid boundary, but it comes with a cost in terms of some extra degrees of freedom.
Also, this method requires a large enough number of degrees of freedom in the
thickness direction of the solid. For second-order elements, typically three elements are
required.
The Solid boundaries to beam points, general connection is used for modeling a beam
with one end “welded” to the face of the solid. You can specify the size of the area on
the solid boundary that is connected to the endpoint of the beam in several different
ways.
The Solid boundaries to beam points, transition coupling is intended for modeling a
transition from a beam to a solid where beam assumptions are valid on both sides of
the connection. Thus, the geometry of the solid at the transition should match the
cross-section data given to the beam.
This connection type can include warping of the solid cross section. In order to
compute the warping properties, an extra PDE is solved over the cross-section
boundaries. To improve the performance, you should preferably solve for these
variables once in a separate stationary study or study step. In that study step, clear all
physics interfaces except the Solid-Beam Connection multiphysics coupling in the Physics
and Variables Selection section.
You can manually control whether to include warping or not. If not included, the setup
of the solver sequence is simplified, but there will be significant stress disturbances
close to the connection boundaries if the cross section is susceptible to warping.
There are four warping variables: one named Warping function and three named
Warping constant. In the successive study steps, you need to manually suppress
The connection can be considered an extension of the Solid boundaries to beam points,
transition coupling in Solid-Beam Connection to also account for radial deformation
of the pipe caused by the fluid pressure and the temperature difference over the cross
section.
The connection can be considered an extension of the Shell edges to beam points
coupling in Shell-Beam Connection to also account for radial deformation of the pipe
caused by the fluid pressure and the temperature difference over the cross section.
Embedded Reinforcement
Lower dimension structural elements can be connected to a solid domain by adding
an Embedded Reinforcement multiphysics coupling. This connection supports
coupling truss, beam, and membrane elements to a Solid Mechanics interface. The
connection can either be rigid, or made by attaching springs between points on the
embedded structure and points in the solid. A more detailed discussion about this type
of modeling is given in Modeling Embedded Structures and Reinforcements.
The underlying theory and more details about the built-in couplings can
be found in
An example could be a shell stiffened by beams. In practice, you would probably use
the built-in coupling described in Beam Edge to Shell Edge (3D) for this case, but the
example shows the principles.
In structure like this, the beam is usually placed at one side of the shell, so that the
centerline of the beam and the midsurface of the shell do not coincide. This difference
Beam centerline
Mathematically, the connection between the beam and the shell can be expressed as
u beam = u shell + X beam – X shell
beam = shell
or equivalently as
u beam = u shell + X beam – X shell n a
beam = shell
Here, is the rotation vector, which contains the rotational degrees of freedom in the
Beam interface. The rotation vector is also available as a variable in the Shell interface,
where it is derived from the rotational degrees of freedom a. The shell normal is
denoted by n. Small rotations are assumed.
1 Add a General Extrusion node under Definitions>Nonlocal Couplings. Select the line on
the shell midsurface as source. Enter data in the Destination Map.
Because a shell does not have a valid rotation degree of freedom around
its normal, the rotation of the beam should not be connected in that
direction.
Loads can be applied in the structural mechanics interfaces on the body, face, edge, or
point levels. You can also apply loads to special features like Rigid Material or Rigid
Connector. There is also an option to apply point loads to given coordinates, which do
not have to coincide with a geometrical point or a mesh node.
In this section:
USING UNITS
Enter loads in any unit, independently of the base SI unit system in the model, because
COMSOL automatically converts any unit to the base SI unit system. To use the
feature for automatic unit conversion, enter the unit in square brackets, for example,
100[lbf/in^2].
The exception is random vibration analysis. In that case, no automatic unit conversions
are available, so you must enter loads in the base units of the model.
Custom coordinate systems are also available and are useful, for example, to specify a
load in any direction without splitting it into components. From the Definitions
toolbar, select a Coordinate System ( ) from the menu.
VISUALIZATION
If you have switched on the physics symbols (see Displaying Physics Symbols in the
Graphics Window — An Example in the COMSOL Multiphysics Reference Manual),
then an applied load is indicated by a symbol together with a coordinate system
indicator displaying the definition directions for the load. The actual direction or
magnitude of the load you enter is not, however, reflected by the symbol. As a load in
COMSOL Multiphysics can be a function of parameters, variables, the solution, or
results from other physics interfaces, it is not possible to display it with only the
information available in the individual load feature.
Once you have turned on the physics symbols for a certain physics interface, you can
fine-tune the display. Every feature which has associated physics symbols will now have
a check box Show physics symbols, by which you can control the display of the symbols
for that specific feature.
Loads are among the results for which predefined plots are generated, so you will
always have access to a visual feedback of the loads after the solution. How to work
with the default load plots is described in detail in the Plotting Applied Loads section.
• Physics Symbols
• Using Units
• Coordinate Systems
• Plotting Applied Loads
Load Cases
For a Stationary or Frequency Domain study, you can define load cases and constraint
cases. Any load or constraint can be assigned to a load or constraint group, and then
be used conditionally.
For most load types, the load case acts as a simple multiplier, but some cases need a
more detailed discussion:
Any expression that acts purely as a load, that is, contributes only to the
right-hand side of the system of equations, can be part of the load case
handling. This is true even if, for example, the corresponding feature does
not have a setting for load cases, or if it is a contribution you have created
using equation based modeling.
For example, if you have a boundary load which partially is always active,
and partially is conditional, you can write 20[MPa]+group.lg1*10[MPa]
in the input field for the pressure.
In particular, this approach is useful for features that override each other,
like Thermal Expansion, since you can then accommodate several load
cases in a single node in the Model Builder.
Singular Loads
In reality, loads always act on a finite area or over a volume. However, in a model loads
are sometimes defined on points or edges, which leads to a singularity. The reason for
this is that points and lines have no area, so the stress becomes infinite. Because of the
stress singularity, there are high stress values in the area surrounding the applied load.
The size of this area and the magnitude of the stress depend on both the mesh and the
material properties. The stress distribution at locations far from these singularities is
unaffected according a to a well-known principle in solid mechanics, the St. Venant’s
principle. It states that for an elastic body, statically equivalent systems of forces
produce the same stresses in the body, except in the immediate region where the loads
are applied.
Figure 2-5 shows a plate with a hole in plane stress loaded with a distributed load and
a point load of the same magnitude. The mesh consists of triangular elements with
quadratic shape functions. The high stress around the point load is dissipated within
the length of a few elements for both mesh cases. The stresses in the middle of the plate
and around the hole agree for the distributed load and the point load. The problem is
that due to the high stress around the singular load it is easy to overlook the high stress
region around the hole. When the point load is applied, the range must be manually
set for the stress plot to get the same visual feedback of the high stress region around
the hole in the two cases. This is because the default plot settings automatically set the
range based on the extreme values of the expression that is plotted.
Despite these findings it is good modeling practice to avoid singular loads because it is
difficult to estimate the size of the singular region. In the Structural Mechanics
The Plasticity and Creep nodes are available as a subnode to Linear Elastic
Material nodes with the Nonlinear Structural Materials Module or the
Geomechanics Module.
Figure 2-5: A plate with a hole subject to a distributed load (left) and a point load (right).
https://www.comsol.com/blogs/
singularities-in-finite-element-models-dealing-with-red-spots/
https://www.comsol.com/blogs/
applying-and-interpreting-saint-venants-principle/
Pressure
A pressure is a load acting toward the normal of a face of the structure. If there are
large deformations in the model and the Include geometric nonlinearity check box is
selected under the Study settings section of the current study step, the pressure acts as
a follower load. The pressure is then defined with respect to the geometry and, as the
geometry deforms locally, the orientation of the load changes. The size of the loaded
area can also change as an effect of straining.
Acceleration Loads
Within the structural mechanics interfaces, you will find four different types of loads
to describe acceleration loads:
• Gravity
• Base Excitation
• Rotating Frame
• Linearly Accelerated Frame
The two first nodes have are of a global type. They will be applied to all features in the
physics interface, and cannot have a spatial variation.
The Rotating Frame and Linearly Accelerated Frame nodes have a domain selection. They
can be applied selectively, and also depend on the coordinates. When applied to a set
of selected domains, such loads are applied also to lower dimensional objects, for
example, a point mass or an added mass on an edge.
All acceleration loads share the property that they are not applied to mass contributions
that belong to global features such as rigid connectors. There, you must add loads
explicitly.
When a separate physics interface is used to model heat transfer in the material, the
entry for the temperature is the dependent variable for the temperature from that
physics interface, typically T. In most cases, possible temperature variables from other
physics interfaces can be directly selected from a list.
• For more information about how to couple heat transfer analysis with
structural mechanics analysis, see Thermal-Structure Interaction. This
module also includes The Thermal Stress, Solid Interface.
• For a detailed discussion about thermal effects in structural mechanics
models, see Thermally Coupled Problems.
Hygroscopic Swelling
Some materials have the capability to absorb significant amounts of moisture through
diffusion processes. Changes in the moisture content may then cause volume changes.
When a separate physics interface is used to model the moisture diffusion in the
material, the entry for the concentration is the dependent variable for the
concentration from that physics interface, typically c. In most cases, possible
concentration variables from other physics interfaces can be directly selected from a
list.
The diffusion of the moisture into the material also adds to the mass density. You can
choose to automatically include this effect in a dynamic analysis, and also in mass
proportional loads, such as gravity and rotating frame loads.
Total Loads
You can specify a load either as a distributed load per unit length, area, or volume, or
as a total force to be distributed on a boundary. In the case of a total load, the applied
distributed load is the given load divided by the area (or length, or volume) on which
its acts. Thus, entering a total load is usually only meaningful when its orientation is
given by a Cartesian coordinate system.
HYDROSTATIC LOAD
Hydrostatic loading is a common special case of spatial variation. In this case, there is
often a fluid surface, above which there is no load. Such a load you can describe with
an expression like if(Z<ZSurf,rhoFluid*g_const*(ZSurf-Z),0). Here, ZSurf
and rhoFluid are assumed to be parameters containing the Z-coordinate of the fluid
surface and the mass density of the fluid respectively.
The local stress state within the loaded element is still limited by what can be described
by the shape functions, but the total load applied on the structure will be more
accurate it you increase the integration order.
TRAVELING LOADS
Loads that are moving along the structure with time can be modeled using an
expression X-v*t, where v is the velocity of the load. The mesh independent point
load of the type Point Load, Free is particularly well suited for this type of modeling. If
a distributed load is modeled using this approach, it is often necessary to increase the
integration order as discussed in the previous section, since the load patch will typically
cover partial element faces.
DESCRIPTION EXAMPLES
In this section:
The number of possible rigid body modes for different geometrical dimensions is
shown in the table below.
TABLE 2-5: NUMBER OF POSSIBLE RIGID BODY MODES
3D 6 (3 translations + 3 rotations)
2D axisymmetric 1 (Z-direction translation)
2D (solid, beam, truss) 3 (2 translations + 1 rotation)
2D (plate) 3 (1 translation + 2 rotations)
For a single body, it is seldom difficult to see whether it is fully constrained or not, but
for a more complex assembly, including several physics interfaces, or advanced
couplings and boundary conditions, it may not be trivial. If you suspect that rigid body
modes are a problem in your model, you can run an eigenfrequency analysis, and check
for modes with zero eigenfrequency as described in Eigenfrequency Analysis.
If there are no constraints which are dictated by the physical boundary conditions of
the structure, you can use the Rigid Motion Suppression feature to automatically
remove the rigid body motions. When you do this, the assumption is that the external
loads are in equilibrium. If not, reaction forces and stress concentrations will appear at
seemingly arbitrary points where the automatic constraints were placed.
As an alternative to applying constraints, you can also add elastic supports through a
Spring Foundation node to suppress rigid body motion.
Orientation
You can specify constraints in global as well as in any previously defined local
coordinate system.
For dynamic analysis, you can also directly prescribe the velocity or acceleration. The
conditions for prescribing displacements, velocities, or accelerations are mutually
exclusive for the same geometrical object since they prescribe the same degree of
freedom.
FREQUENCY DOMAIN
In frequency domain, a prescribed velocity vp or prescribed acceleration ap can be
directly interpreted as a prescribed displacement up:
vp
u p = ------
i
–ap
u p = --------
2
-
A prescribed velocity with zero phase is assumed to have its peak at the
reference phase. As an effect, the corresponding peak displacement is
shifted by 90°. Similarly, a positive prescribed acceleration with zero phase
corresponds to a negative value of the displacement.
TIME DOMAIN
In the case of a time-dependent analysis, the prescribed displacement is obtained as
u p t = u 0 t 0 + v p d
t0
or
where u0 and v0 are is given by the initial conditions. It is not possible to set explicit
initial conditions, but if initial values are taken from a previous study, they will be
respected. In order to compute the integrals, up is introduced as a separate degree of
freedom which is solved for by adding an extra ODE.
As prescribing the velocity or acceleration in time domain comes with an extra cost,
you should always consider using a prescribed displacement instead. As long as the
time history of the velocity or acceleration is a known a priori and does not depend on
the solution itself, this is possible.
• When the velocity or acceleration has a simple time dependence, you can integrate
it analytically one or two times to obtain the displacement, and directly prescribe the
displacement instead.
• When you have complicated known velocity or acceleration histories, for example
from measurements, you can use the integrate() operator. In this case, you enter
the prescribed displacement as integrate(my_data(tau),tau,0,t). Here
my_data is the measured data as function of time, and tau is a dummy integration
variable
STATIONARY ANALYSIS
In a stationary analysis, the prescribed velocity and acceleration nodes can have two
different behaviors. As a default, they are ignored, but you can also select that the
degrees of freedom having a prescribed velocity or acceleration in a dynamic analysis
should be constrained to zero in a static analysis.
Symmetry Constraints
In many cases symmetry of the geometry and loads can be used to your advantage in
modeling. Symmetries can often greatly reduce the size of a model and hence reduce
the memory requirements and solution time. When a structure exhibits axial
symmetry, use the axisymmetric physics interfaces. A solid that is generated by rotating
a planar shape about an axis is said to have axial symmetry. In order to make use of
For other types of symmetry, use the predefined symmetry and antisymmetry
constraints. This means that no expressions need to be entered — instead just add the
type of constraint to apply to the model.
If the geometry exhibits two symmetry planes (Figure 2-6), model a quarter of the
geometry by using the Symmetry node for the two selected surfaces.
Figure 2-6: If the geometry exhibits two symmetry planes, model a quarter of the geometry
by using the Symmetry feature for the two selected surfaces.
Both geometric symmetry and loads are important when selecting the
correct constraints for a model.
SYMMETRY IN 2D AXISYMMETRY
In an axisymmetric model, the only possible symmetry is when the symmetry plane is
normal to the Z-axis. For models in 2d axisymmetry, the Symmetry Plane node is used
for prescribing this type of symmetry.
You can modify the symmetry condition, so that it can translate in various ways by
using the controls in the Normal Direction Condition section of the settings for the
Symmetry constraint. You can model the following cases:
When using nodal constraints, one constraint is generated for each node within the
selection a certain constraint feature. With elemental constraints, the number of
constraints added at a node equals the number of elements connected to that node.
This means that if some values used in the constraints differ between the elements,
then different constraints will be generated by the elemental method, whereas with the
nodal method an average is computed at the node before adding the constraint.
When several constraints are present at a node, the internal constraint elimination
algorithm is responsible for reducing them to a minimum unique set. Using elemental
constraints will clearly put an extra burden on this algorithm, so whenever possible you
should use nodal constraints.
The two different options exist, since under some circumstances the actual constraints
can differ between the two methods. Consider for example a symmetry constraint,
where the displacement in the direction normal to the boundary is constrained by the
equation
un = 0
• If both boundaries are selected in the same Symmetry node, then only a single
constraint is applied for each node along the common edge, while you actually want
constraints along the normals of both planes. The normal used would be pointing
somewhere between the two planes, since a nodal constraint uses averaging of the
values from the adjacent elements.
• If two Symmetry nodes are used, so that the selection in any one of them only
contains boundaries without a normal direction discontinuity, the intended
constraints are added. On the common edge, there will be two contributions, one
from each Symmetry node, and each using the normal direction of its boundary. If
you want to use nodal constraints, you must set up your model in this way if the
constraints are orientation dependent.
Elemental constraints, on the other hand, can cause problems if the constraints added
by adjacent elements are not exactly the same. This could for example happen if the
normal orientation differs between neighboring elements. In such a case, a boundary
could behave as if it were fixed when a Symmetry, Antisymmetry, or Roller constraint is
applied. Such a situation could occur when the component consists of an imported
mesh, so that no underlying geometry exists.
The default type of the constraint, nodal or elemental, differs between different
constraint features. A nodal formulation is the default whenever it is considered safe,
like for a Fixed Constraint. Whenever the constraint can have a dependency on the
surface orientation, the default value is elemental.
For most constraints in the structural mechanics interfaces, you have the possibility to
select that certain objects of lower dimensions should be excluded from the main
selection. To do this, you must first select Advanced Physics Options. In the settings for
a constraint, like for example Prescribed Displacement, new sections named Excluded
In the structural mechanics interfaces, there are many types of complex constraints,
and sometimes you may get conflicts or duplicates which makes the model either
overconstrained, or problematic for the automatic constraint elimination algorithm. If
you are aware of such situations, it is good practice to remove one of the potentially
conflicting constraints. One example of such a situation is when you have a Solid-Shell
Connection meeting a symmetry plane, as shown in Figure 2-8.
Here you would add a Symmetry condition on a boundary in the Solid Mechanics
interface, as well as a Symmetry condition on an edge in the Shell interface. But at the
same time, the displacements on whole boundary where the solid meets the shell are
controlled by shell degrees of freedom as an effect of the Solid-Shell Connection. As a
result, on the edge marked with Conflict in the sketch, the displacements will be
controlled both by the symmetry condition is Solid Mechanics, and implicitly through
the coupling, by the symmetry condition in the Shell interface. Particularly if the
geometry is curved, there is a risk that these constraints are not identical from a
numerical point of view. In this case, excluding the conflicting edge from the selected
boundary in the Solid Mechanics interface will make the behavior unique and fully
predictable.
Another example where constraints will come in conflict is if you want to constrain the
displacement on parts of the geometry using weak constraints, while keeping the
default pointwise constraints on other parts. If the same mesh node has both types of
See also Excluded Surfaces, Excluded Edges, and Excluded Points in the
COMSOL Multiphysics Reference Manual.
Kinematic Constraints
Kinematic constraints are equations that control the motion of solids, faces, edges, or
points. Add a Prescribed Displacement constraint to enter expressions for constraints.
You can define the equations using predefined coordinate systems as well as custom
coordinate systems. Special constraints, for instance to keep an edge of body straight
or to make a boundary rotate, require such constraint equations.
In the 3D and 2D Solid Mechanics interfaces and in the Shell and Beam
interface there is a special constraint called a Rigid Connector. A rigid
connector is applied to one or more boundaries, edges, or points and force
them to behave as connected to a common rigid body. The rigid
connector can be given prescribed displacements and rotations and thus
simplifies the realization of some constraints.
Rotational Joints
Joints between elements in The Truss Interface are automatically rotational joints
because the truss elements have no rotational degrees of freedom. For beams, however,
the rotational degrees of freedom are by default coupled between elements. To create
a rotational joint between two beam elements, add one additional Beam interface to a
geometry. Make sure that it is only active for the edge that includes the point where
the joint is positioned and that no other physics interface is active here. Couple the
translational degrees of freedom and leave the rotational degrees of freedom
uncoupled at the joint.
Attachments
An Attachment is a set of boundaries, edges, or points on a flexible or rigid component
used to connect it to another flexible or rigid component through a joint or spring. An
When an attachment is connected to a flexible body, you can use two different
formulations: rigid or flexible.
In the rigid attachment formulation, all selected boundaries or edges behave as if they
were connected to a common rigid body. This may cause an unrealistic stiffening and
local stress concentrations.
In the flexible version, the boundaries are allowed to deform, and the rigid body
constraints are enforced only in an average sense.
The attachment formulation is similar to that of a rigid connector. In the rigid case,
the only degrees of freedom needed to represent this assembly are the ones needed to
describe the movement of a rigid body. In 2D this is just two in-plane translations, and
the rotation around the out-of-plane axis.
In 3D the situation is more complex. Six degrees of freedom are necessary, usually
selected as three translations and three parameters for the rotation. For finite rotations
any choice of three rotation parameters is singular at some specific set of angles. For
this reason, a four-parameter quaternion representation is used.
Some extra degrees of freedom are added for each attachment where the flexible
formulation is used.
When an attachment is defined on a rigid component, it does not create any degrees
of freedom of its own and directly picks the degrees of freedom of the rigid
component.
The following sections describe the merits and costs of these methods.
Reaction forces are computed as the sum of the nodal values over the selected volume,
face, or edge. Reaction moments are calculated as the sum of the moment from the
reaction forces with respect to a reference point, and any explicit reaction moments (if
Since the reaction force variables are added to the solution components, the number
of DOFs for the model increases slightly, depending on the mesh size for the
boundaries in question. Boundaries that are adjacent to each other must have the same
constraint settings. The reason for this is that adjacent boundaries share a common
node.
Using weak constraints affects the structure of the equation system to be solved, and
is not suitable for all types of equation solvers.
The first type, contained in the variables interface.Tax, is computed from the
stresses. It is always available. Since the surface traction vector is based on computed
stress results, this method is less accurate for computing reactions than the other
methods.
You can also add a material model which you have coded yourself and made available
as a binary library file using an External Stress-Strain Relation.
Many of the material models can be augmented by effects like thermal expansion,
hygroscopic swelling, initial stresses and strains, external stress, and activation.
1 An elastic strain is computed by removing all inelastic strains (for example, plastic or
thermal strains) from the total strain.
2 An “elastic stress” is computed from the elastic strains.
3 Any additional stresses (for example viscous stresses, or initial stresses) are added to
form the total stress.
This concept will give you a great freedom in combining effects. Some such useful
combinations are
Constitutive matrices, such as the elasticity matrix for an anisotropic material, are in
many cases per definition symmetric. Only the upper diagonal of the matrix given as
input is used for forming the matrix used, so you need not enter the lower diagonal
part.
For isotropic linear elasticity, two parameters are enough to describe the material
behavior. The number of parameters increases to (at most) 21 for the fully anisotropic
case in 3D. When setting up a model, make sure that the material parameters are
defined in agreement with the type of relationship used. If necessary, transform the
material data before entering it in the physics interface. For example, for orthotropic
materials calculate the Poisson’s ratio xy by
Ex
xy = yx ------
Ey
• Linear Viscoelasticity
• Viscoelasticity
In several material models you will find settings named Use mixed formulation or
Compressibility, by which you can introduce a mixed formulation.
Use a mixed formulation when the material data is such that the deformation is close
to being incompressible. For an isotropic elastic material, this happens when Poisson’s
ratio approaches 0.5.
Note that in the case of a mixed mesh, you may have to sub-divide the
domain by element type, to be able to control the shape functions for the
auxiliary variable.
DISCONTINUOUS LAGRANGE
If this shape function type is selected, the auxiliary variable shape function is
discontinuous across element boundaries, and it is one order lower than the shape
function order of the displacements.
LINEAR
If this shape function type is selected, the auxiliary variable is regarded as a linearly
interpolated field in the element. In the case of an auxiliary pressure, the interpolation
in a 3D isoparametric element is
p = p0 + 1 p1 – p0 + 2 p2 – p0 + 3 p3 – p0
where p0, p1, p2, and p3 are auxiliary pressure coefficients, and 1, 2, and 3 are
isoparametric coordinates. Note that the field is linear in the local element coordinates,
and that it is not continuous across element boundaries.
The linear shape function type is not available for layered material
features.
CONSTANT
This shape function type represents a constant auxiliary variable in the element, p = p0.
The mixed formulation is useful not only for linear elastic materials but
also for nonlinear elastic materials, elastoplastic materials, hyperelastic
materials, and viscoelastic materials.
Note that some iterative solvers do not work well together with mixed
formulation because the stiffness matrix becomes indefinite.
When using the crack band method or no regularization at all, the following steps are
recommended:
• The size of the biggest mesh element h should not exceed 2EGf/ts2, where E is
the Young’s modulus, Gf is the fracture energy per unit area, and ts is the tensile
strength. Larger values of h will cause a snap-back of the stress-strain curve at the
material point level.
• Use linear shape functions for the displacement field. When using higher-order
shape functions, strains may localize in either rows of Gauss points, or entire
elements, depending on the stress state and numerical rounding errors.
• If cracks are located on a symmetry plane, the model parameters should be modified
so that the amount of dissipated energy is reduced by one half in the elements
adjacent to the symmetry plane. This can be achieved, for example, by reducing the
fracture energy in that row of elements.
When using the implicit gradient method, the element size should be sufficiently small
to resolve damaged zones. The same applies to the Phase field damage model, where it
is recommended that size of the mesh elements in the expected crack path follows
h < lint/2 for a linear displacement field, otherwise h < lint.
1 Initialization. The crack phase field, displacement field and state variables are
known at step n.
2 Update state variables. Update internal state variables used by the phase field
model with values from step n.
3 Solve for the Crack Phase Field. Compute the crack phase field variable in a Newton
step, with the displacement field frozen at step n.
4 Solve for the Displacement field. Compute displacement field variables in a Newton
step with the updated crack phase field.
This leads to a single-pass algorithm that is accurate only for sufficiently small
parameter or time steps.
The default solver will suggest the above single-pass algorithm for the
Phase field damage model when it is feasible to perform the operator split.
Cases where this is not possible include when some multiphysics
couplings are present in the model and when a segregated contact
algorithm has to be used.
Before moving to implementing your own material model, there are however two
other options to consider:
• Many material models like hyperelasticity, creep and plasticity have User defined as
one of the options in addition to the standard models. Any material model which
you can describe using built-in variables is most conveniently described here.
• A material model which can partially described in terms of a PDE can often be
implemented using one of the mathematical interfaces. Stresses or strains computed
in that interface are then injected in an existing material model using the External
Stress and External Strain subnodes.
There are two basic types of external material functions: those which completely
replace other material definitions in a domain, and those that just compute an inelastic
strain contribution to be used as part of an existing material model. The former is
During the solution, an external material routine is always called for each Gauss point
during evaluation of stiffness matrices and computation of residuals. During result
presentation, the external material can be called from any location in the geometry, as
requested by for example graphs and point evaluations.
Almost invariably, you need to store state variables in the external material, such as for
example plastic strains. The state variables are stored at the Gauss points. If an external
material is called at another location, the state variables will be interpolated to that
location. This means that the state of the material may not be exactly consistent there,
which can lead to some artifacts during result presentation. You can avoid this problem
by using the gpeval operator.
References
1. D. Chapelle and K.J. Bathe, The Inf-Sup test, Comp. Struct., Vol. 47, No. 4/5, pp.
537–545, 1993.
Within a piezoelectric, there is a coupling between the strain and the electric field,
which is determined by the constitutive relation:
T
S = sET + d E
(2-14)
D = dT + T E
Here, the naming convention used in piezoelectricity theory is assumed: S is the strain,
T is the stress, E is the electric field, and D is the electric displacement field. The
material parameters sE, d, and T, correspond to the material compliance, the coupling
properties and the permittivity. These quantities are tensors of rank 4, 3, and 2,
Equation 2-14 will, using the notation from structural mechanics, then read
T
= sE + d E
(2-15)
D = d + 0 rT E
Equation 2-14 (or Equation 2-15) is known as the strain-charge form of the
constitutive relations. The equation can be re-arranged into the stress-charge form,
which relates the material stresses to the electric field:
T
= cE – e E
(2-16)
D = e + 0 rS E
The material properties, cE, e, and S are related to sE, d, and T. It is possible to use
either form of the constitutive relations. In addition to Equation 2-14 or
Equation 2-16, the equations of solid mechanics and electrostatics must also be solved
within the material.
• Piezoelectric Coupling
• Modeling Piezoelectric Problems
• Piezoelectricity in the theory section
The stiffness, compliance, coupling, and dielectric material property matrices are
defined with the crystal axes aligned with the local coordinate axes. Note that the signs
of several matrix components differ between the 1949 and the 1978 standards (see
Table 2-8). In the absence of a user-defined coordinate system, the local system
corresponds to the global X, Y, and Z coordinate axes. When an alternative coordinate
system is selected this system defines the orientation of the crystal axes. This is the
mechanism used in COMSOL to define a particular crystal cut, and typically it is
necessary to calculate the appropriate Euler angles for the cut (given the thickness
orientation for the wafer). All piezoelectric material properties are defined using the
Voigt form of the abbreviated subscript notation, which is universally employed in the
literature (this differs from the standard notation used for the Solid Mechanics
interface material properties). The material properties are defined in the material
frame, so that if the solid rotates during deformation the material properties rotate
with the solid. See Modeling Geometric Nonlinearity.
TABLE 2-8: SIGNS FOR THE MATERIAL PROPERTIES OF QUARTZ, WITHIN THE TWO STANDARDS COMMONLY
EMPLOYED.
s14 + + - -
c14 - - + +
d11 - + + -
d14 - + - +
e11 - + + -
e14 + - + -
TABLE 2-9: CRYSTAL CUT DEFINITIONS FOR QUARTZ CUTS WITHIN THE TWO STANDARDS COMMONLY
EMPLOYED AND THE CORRESPONDING EULER ANGLES FOR DIFFERENT ORIENTATIONS OF THE CRYSTAL
THICKNESS.
There are also a number of methods to define the local coordinate system with respect
to the global system. Usually, it is most convenient to define the local coordinates with
a Rotated System node, which defines three Euler angles according to the ZXZ
convention (rotation about Z, then X, then Z again). Note that these Euler angles
define the local (crystal) axes with respect to the global axes — this is distinct from the
approach of defining the cut (global) axes with respect to the crystal (local) axes.
The hysteretic electrical losses are usually used to represent high frequency electrical
losses that occur as a result of friction impeding the rotation of the microscopic dipoles
that produce the material permittivity.
• Using the Rayleigh Damping option in the Mechanical Damping and Coupling Loss
subnodes to the Piezoelectric Material,.
Rayleigh Damping
• Using either the Dispersion, or Complex permittivity, or Maximum Loss Tangent option
in the Dielectric Loss subnode to the Piezoelectric Material.
Rayleigh Damping
HYSTERETIC LOSS
In the frequency domain, the dissipative behavior of the material can be modeled using
complex-valued material properties, irrespective of the loss mechanism. Such hysteretic
losses can be applied to model both electrical and mechanical losses. For more
information about hysteretic losses, see Ref. 1 to Ref. 4.
T = c˜ E S – e˜ E
T
D = e˜ T + ̃ E
0 rS
S = s˜ E T + d˜ E
T
D = d˜ S + ̃ E
0 rT
All constitutive matrices in the above equation can be complex-valued matrices, where
the imaginary parts define the dissipative function of the material.
Both the real and complex parts of the material data must be defined so as to respect
the symmetry properties of the material being modeled and with restrictions imposed
by the laws of physics.
A key requirement is that the dissipation density is positive; that is, there
is no power gain from the passive material. This requirement sets rules for
the relative magnitudes for all material parameters. This is important
when defining the coupling losses.
In COMSOL, you can enter the complex-valued data directly or by means of loss
˜
factors. When loss factors are used, the complex data X is represented as pairs of a
real-valued parameter
˜
X = real X
˜ ˜
X = imag X real X
the ratio of the imaginary and real part, and the complex data is then:
˜
X = X 1 j X
The loss factors are defined so that a positive loss factor value usually corresponds to a
positive loss. The complex-valued data is then based on sign rules.
By default, there is no damping until at least one of the damping and losses related
subnodes is added.
For the Piezoelectric Material node, the following equations apply via the
corresponding three subnodes:
Mechanical Damping
m n m n m n
c˜ E = c E 1 + j cE
m n m n m n
s˜ E = s E 1 – j sE
Coupling Loss
m n m n m n
e˜ =e 1 + j e
m n m n m n
d˜ =d 1 + j d
Dielectric Loss
m n m n m n
̃ rS = rS 1 – j S
m n m n m n
̃ rT = rT 1 – j T
The loss factors can also be entered as scalar isotropic factors independently of the
material and the other coefficients.
Two important particular cases are available in COMSOL, for which the
dielectric loss data can be entered based on experimentally measurable
quantities:
• Complex permittivity
• Maximum Loss Tangent
DIELECTRIC DISPERSION
The Dielectric Loss subnode can be set to use the Dispersion option. In such case, the
following equations need to be solved in the time domain:
D
+ J p = 0 (2-17)
t
J p E
Jp + d = 0 rS (2-18)
t t
where you can specify two material parameters: the relaxation time d and the relative
permittivity increment rS. The latter can be either a matrix or a scalar quantity. This
model is a one-term version of the more general Debye dispersion model, Ref. 13.
D = eS + 0 E (2-19)
where S is the strain tensor, and is the relative permittivity in the high frequency
limit (that is, for excitations with a characteristic time much shorter than the relaxation
time d).
The parent Piezoelectric Material node has an input for the relative
permittivity, rS, which is used in stationary study. You can chose how this
input will be interpreted in the dispersion computations. The options are:
D
D + d + 0 rS E = 0
t
This is the equation form used in COMSOL Multiphysics for time-dependent analysis.
For the eigenfrequency and frequency domain analyses, the corresponding equation is:
1 + j d jD + 0 rS jE = 0
In most cases, i can be factored out, so that the following equation is solved:
D + j d D + 0 rS E = 0
This equation, together with the constitutive relation Equation 2-19 gives
where
rS
' = + --------------------------2-
1 + d
and
d 0
rS
'' = ---------------------------2-
1 + d
Equation 2-20 shows how the dispersion parameters contribute to the polarization
and losses. Thus, the effective relative permittivity decreases with the excitation
frequency from the low frequency limit + rS down to the high frequency limit
. The damping effect vanishes for both large and small frequencies, and it reaches
the maximum for 1d.
The following two sections present two cases, for which the dielectric dispersion data
can be related to other experimentally measurable quantities. Both cases, can be used
Complex permittivity
In this case, the complex relative permittivity is known at a certain reference frequency
–1 –1
d = ref '' ' –
2 –1
rS = ' – + '' ' –
If the relative permittivity rS (input on the Piezoelectric Material parent node) is
selected to represent the low frequency limit, one has
2 –1
= ' – '' rS – '
If rS is selected to represent the high frequency limit, one can simply use = rS
instead.
r f = ' I – j
For many materials, the loss tangent reaches a maximum at certain frequency fref
within the frequency range of interested
–1 2 12
d = ref max + max + I
If the relative permittivity rS (input on the Piezoelectric Material parent node) is
selected to represent the low frequency limit, the relative permittivity contribution is
computes as
rS = rS –
–2
where = d ref rS .
If the relative permittivity rS is selected to represent the high frequency limit, it is
computed as
2
rS = d ref – I
where = rS .
jD + J e = 0
Je = e E
where e is the material electrical conductivity, and E is the electric field. The above
form of the equation is used for the eigenfrequency analysis in COMSOL
Multiphysics.
When the conduction loss is applied, the default eigenvalue analysis will
in most cases return a number of pure imaginary eigenfrequencies. To
avoid this, you can configure the solver to search for eigenvalues with real
part larger than zero.
Do not use any dielectric loss factor damping together with the
conduction loss in an eigenfrequency analysis.
For the frequency domain analysis, the angular frequency is just a parameter, and the
equation can be transformed into
–1
D – j e E = 0
which allows you to use both a dielectric loss factor and electrical conductivity in a
frequency response study. In such case, ensure that the loss factor refers to the
alternating current loss tangent, which dominates at high frequencies, where the effect
of ohmic conductivity vanishes (Ref. 7).
Conduction loss can be combined with Dielectric Dispersion for both eigenfrequency
and frequency domain analyses. The following equation forms are used, respectively,
in the frequency domain:
–1
D + j d D + 0 rS E + d – j J e = 0
1 + j d jD + 0 rS jE + 1 + j d J e = 0
6. P.C.Y. Lee and N.H. Liu, “Plane Harmonic Waves in an Infinite Piezoelectric Plate
With Dissipation,” Frequency Control Symposium and PDA Exhibition, pp. 162–
169, IEEE International, 2002.
11. B.A. Auld, Acoustic Fields and Waves in Solids, Krieger Publishing, 1990.
13. N.S. Stoykov, T.A. Kuiken, M.M. Lowery, and A. Taflove, “Finite-element
time-domain algorithms for modeling linear Debye and Lorentz dielectric dispersions
at low frequencies,” IEEE Transactions on Biomedical Engineering, vol. 50, no. 9,
pp. 1100–1107, 2003.
Piezoelectric Coupling
The piezoelectric effect is an interaction between the mechanical and electrical physics,
where a stress applied on a piezoelectric material generates a voltage (direct effect) or
a voltage applied on it generates the deformation of the material (inverse effect). In
COMSOL Multiphysics, the Piezoelectricity interface is constituted of one Solid
Mechanics and one Electrostatics interface, which are coupled together by a Piezoelectric
Effect multiphysics feature. Hence a piezoelectric problem contains solid and
electrostatic domains, with at least one domain shared by the two physics interfaces
and with the piezoelectric coupling defined on it.
In the first two cases, by default all the domains in the model are assumed to be
piezoelectric materials.
When setting up the problem manually (that is, by adding single physics interfaces one
at a time) both Solid Mechanics and Electrostatics interfaces should be added. Then, you
have to specify which domains are in each physics, and which domains are to be
modeled as piezoelectric materials.
1 On the Solid Mechanics interface Settings window, locate the Domain Selection
section. Select the domains which undergo structural deformation, including the
piezoelectric material domains.
2 Go to the Solid Mechanics>Piezoelectric Material node (if it is not yet available, add
it). Select the domains where the piezoelectric effect applies. Domains which are not
piezoelectric can be modeled using other available material models.
3 Go to the Electrostatics node and under Domain Selection select the domains the
electrostatics equations must be solved. These domains include all the piezoelectric
domains, and any other insulating domains.
• Coordinate System Selection section: The material is poled in the x3 direction of the
coordinate system (x1, x2, x3) specified in this section. By default, it is set to the
global coordinate system. If the piezoelectric material is poled along another
direction, you need to define a coordinate system so that its third direction is aligned
• Mechanical Damping: Specify the domains of application, then choose if you want to
define a loss factor for cE, a loss factor for sE (in Strain-charge form), an isotropic
loss factor, or a Rayleigh damping.
• Dielectric Loss: Specify the domains of application, then choose if you want to define
a loss factor for rS, a loss factor for rT (in Strain-charge form), or dispersion.
• Coupling Loss: Specify the domains of application, then choose if you want to define
a loss factor for e, a loss factor for d (in Strain-charge form), or Rayleigh damping.
• Conduction Loss (Time-Harmonic): Specify the domains of application, then choose
how you want to define the Electrical conductivity.
MECHANICAL PROPERTIES
3
• Density rho (SI unit: kg/m )
• Elasticity matrix cE (SI unit: Pa) in Stress-charge form.
• Compliance Matrix sE (SI unit: 1/Pa) in Strain-charge form.
ELECTROSTATIC PROPERTIES
• Relative Permittivity rS (dimensionless) in Stress-charge form.
• Relative Permittivity rT (dimensionless) in Strain-charge form.
COUPLING PROPERTIES
• Coupling matrix eES (SI unit: C/m2) in Stress-Charge form.
• Coupling matrix dET (SI unit: C/N) in Strain-charge form.
Ferroelectricity
The ferroelectroelasticity and ferroelectricity phenomena are related to phase
transitions in materials. In its ferroelectric phase, the material exhibits spontaneous
polarization, so that it is constituted of domains with nonzero polarization even at zero
applied field. This is similar to permanent magnetism in ferromagnetics, which explains
the name used for such materials. Electrostriction in ferroelectroelastic materials can
be related to the domain rotation. Thus, the applied electric field can both rearrange
the domains resulting into the net polarization and rotate the domains mechanically.
Thus, the material extends in the direction of the electric field and contracts in the
direction perpendicular to the field. The domain rotation can be affected by an applied
mechanical stress, which also results into the effective polarization. At very large
electric fields, the electrostrictive effect saturates, as all ferroelectric domains in the
material are aligned along the direction of the applied field. Domain wall interactions
can also lead to a significant hysteresis in the polarization and strain.
Multiphysics Interfaces
In COMSOL Multiphysics, the electrostrictive effect can be modeled using two
multiphysics interfaces: Electrostriction and Ferroelectroelasticity. Both interfaces are
constituted of one Solid Mechanics and one Electrostatics interface, which are coupled
The magnetostrictive strain has a nonlinear dependence on the magnetic field and the
mechanical stress in the material. However, the effect can be modeled using linear
coupled constitutive equations if the response of the material consists of small
deviations around an operating point (bias point). This type of coupling is reffed to as
Piezomagnetic Effect.
In COMSOL Multiphysics, there are two multiphysics interfaces for modeling either
linear or nonlinear magnetostriction, respectively:
The nonlinear model of magnetostrictive strain can be used for the whole range from
full demagnetization to saturation magnetization. In case of nonlinear
magnetostriction, the magnetization model can be selected. The following options are
For more details, see the corresponding theory section in Magnetostriction and
Piezomagnetism.
In the first two cases, by default all the domains in the model are assumed to be
magnetostrictive materials.
PIEZOMAGNETISM
When setting up the problem manually (that is, by adding single physics interfaces one
at a time) both Solid Mechanics and Magnetic Fields interfaces should be added. Then,
you have to specify which domains are in each physics, and which domains are to be
modeled as piezomagnetic materials.
1 On the Solid Mechanics interface Settings window, locate the Domain Selection
section. Select the domains which undergo structural deformation, including the
piezomagnetic material domains.
2 Go to the Solid Mechanics>Piezomagnetic Material node. If it is not yet available, add
such node. Select the domains where the piezomagnetic effect applies. You find all
the necessary material data inputs within the Piezomagnetic Material. This includes
the elasticity, magnetic permeability and coupling matrices. Non piezomagnetic
domains can be modeled using any other available structural material model.
3 Go to the Magnetic Fields node and under Domain Selection select the domains the
magnetics equations must be solved. These domains include all the piezomagnetic
domains, and any other magnetic domains.
Only domains that have both Ampere’s Law, Piezomagnetic selected in the
Magnetic Fields interface and Piezomagnetic Material selected in the Solid
Mechanics interface are applicable in the selection for Piezomagnetic Effect
coupling. The selection of this multiphysics coupling node cannot be
edited directly.
• Piezomagnetic Material
• Ampère’s Law, Piezomagnetic
NONLINEAR MAGNETOSTRICTION
When setting up the problem manually (that is, by adding single physics interfaces one
at a time) both Solid Mechanics and Magnetic Fields interfaces should be added. Then,
you have to specify which domains are in each physics, and which domains are to be
modeled as magnetostrictive materials.
Instead of Magnetic Fields interface, you can also use Rotating Machinery,
Magnetic interface.
1 On the Solid Mechanics interface Settings window, locate the Domain Selection
section. Select the domains which undergo structural deformation, including the
magnetostrictive material domains. The following structural material nodes are
supported to represent the magnetostrictive domains: Linear Elastic Material and
Hyperelastic Material (available with the Nonlinear Structural Materials Module).
2 Go to the Magnetic Fields node and under Domain Selection select the domains the
magnetics equations must be solved. These domains include all the magnetostrictive
domains, and any other magnetic domains.
3 Go to the Magnetic Fields>Ampere’s Law, Nonlinear Magnetostrictive node. If it is not
yet available, add such node. Select the domains where the magnetostriction effect
needs to be modeled. You can specify the magnetization model and enter the related
material data on the node.
4 A Multiphysics>Nonlinear Magnetostriction node is already present if the coupling was
added using either the Model Wizard or Add Physics window. If the model is set up
manually (that is, single physics interfaces are added), right-click the Multiphysics
Couplings node to add a Nonlinear Magnetostriction coupling. If several Solid
Mechanics or Magnetic Fields interfaces are present, select the correct ones. On the
coupling node, you can also select the magnetostriction model depending on the
material symmetry and enter the corresponding coupling data. By default, the
model will be solved as fully coupled. Alternatively, on the Nonlinear Magnetostriction
In this section:
• About Damping
• Rayleigh Damping
• Loss Factor Damping
• Viscoelastic Damping
• Explicit Complex-Valued Damping
• Equivalence Between Loss Factor, Rayleigh, and Viscous Damping
• Piezoelectric Damping
• Adding Damping in the Modal Solver
Damping Sources
There are many sources of damping in a system. Some of them are:
It is often difficult to separate and quantify these effects, so damping modeling is one
of the biggest challenges in structural dynamics.
2
d u du
m ---------2- + c ------- + ku = f t (2-21)
dt dt
In this equation u is the displacement of the degree of freedom, m is its mass, c is the
damping parameter, and k is the stiffness of the system. The time (t) dependent forcing
term is f(t). This equation is often written in the form:
2
d u- du 2 f t
--------- + 2 0 ------- + 0 u = --------- (2-22)
dt
2 dt m
where c2m0 and 02 km. In this case is the damping ratio ( 1 for
critical damping) and 0 is the undamped resonant frequency of the system. In the
literature it is more common to give values of than c. The damping ratio can also
be readily related to many of the various measures of damping employed in different
disciplines. These are summarized in Table 2-10.
TABLE 2-10: RELATIONSHIPS BETWEEN MEASURES OF DAMPING
Logarithmic u t0 d 2
decrement d = ln ----------------------
u t0 +
« 1
where t0 is a reference time and is the period
of vibration for a decaying, unforced degree of
freedom.
Quality factor Q = Q 1 2
In the frequency domain, the time dependence of the force and the displacement can
be represented by introducing a complex force term and assuming a similar time
dependence for the displacement. The equations
jt jt
f t = Re Fe and u t = Re Ue
are written where is the angular frequency and the amplitude terms U and F can in
general be complex (the arguments provide information on the relative phase of
signals). Usually the real part is taken as implicit and is subsequently dropped.
Equation 2-21 takes the following form in the frequency domain:
2
– mU + jcU + kU = F (2-23)
where the time dependence has canceled out on both sides. Alternatively, this equation
can be written as:
2 2 F
– U + 2j 0 U + 0 U = ----- (2-24)
m
There are three basic damping models available in the structural mechanics interfaces
for explicit modeling of material damping — Rayleigh damping, viscous damping,
and loss factor models based on introducing complex quantities into the equation
system. There are also other phenomena which contribute to the damping. Some
material models, such as viscoelasticity and plasticity are inherently dissipative. It is also
possible to model damping in spring conditions.
Rayleigh damping introduces damping in a form based on Equation 2-21. This means
that the method can be applied generally in either the time or frequency domain. The
parameter c in Equation 2-21 is defined as a fraction of the mass and the stiffness using
two parameters, dM and dK, such that
c = dM m + dK k (2-25)
Substituting this relationship into Equation 2-21 and rearranging into the form of
Equation 2-22 gives:
2
d u- 2 du 2 ft
---------
2
+ dM + dK 0 ------- + 0 u = ---------
dt dt m
When there are many degrees of freedom m, k, and c become matrices and the
technique can be generalized.
1 dM
= --- ----------- + dK 0 (2-26)
2 0
Note that Equation 2-26 holds separately for each vibrational mode in the system at
its resonant frequency. In the frequency domain it is possible to use frequency
dependent values of dM and dK. For example, setting dM 0 and dK 2/0
produces an equivalent viscous damping model at the resonant frequency 0.
While Rayleigh damping is numerically convenient, the model does not agree with
experimental results for the frequency dependence of material damping over an
extended range of frequencies. This is because the material damping forces behave
more like frictional forces (which are frequency independent) than viscous damping
forces (which increase linearly with frequency as implied by Equation 2-23). In the
frequency domain it is possible to introduce loss factor damping, which has the desired
property of frequency independence.
Using Equation 2-26, this relationship at two frequencies, f1 and f2, with different
damping ratio, 1 and 2, results in an equation system that can be solved for dM and
dK:
1-
----------- f
4f 1 1 dM 1
=
1- dK 2
----------- f
4f 2 2
1 f2 – 2 f1
dM = 4f 1 f 2 ---------------------------
2 2
f2 – f1
2 f2 – 1 f1
dK = ---------------------------
2 2
f2 – f1
Using the same damping ratio, 1 = 2 = 0, does not result in a constant damping
factor inside the interval f1 f f2. It can be shown that the damping factor is lower
inside the interval, as Figure 2-13 shows.
Damping ratio
Rayleigh damping
(f)
0
Specified damping
f1 f2 f
f1 2 f2
---- ----- ----
f2 1 f1
For many applications it is sufficient to leave dM as zero and to define damping only
using the dK coefficient. Then, according to Equation 2-26, a damping which
increases linearly with frequency is obtained. If the damping ratio f0 or loss factor
f0 is known at a given frequency f0, the appropriate value for dK is:
dK = f 0 = 2f 0
c
D = 1 + j s D
The use of loss factor damping traditionally refers to a scalar-valued loss factor s. But
there is no reason that s must be scalar. Because the loss factor is a value deduced from
true complex-valued material data, it can be represented by a matrix of the same
dimensions as the anisotropic stiffness matrix. Especially for orthotropic materials,
there should be a set of loss factors of all normal and shear elasticity modulus
components, and COMSOL allows all these options, so a more general expression is.
c
D mn = 1 + j s mn D mn
W s
S q = j s ----------
E
Loss factor damping is available for frequency response analysis and damped
eigenfrequency analysis in all interfaces.
Viscous Damping
You can add an explicit viscous damping to several material models. Viscous damping
can be used both in time domain and frequency domain. In the viscous damping
model, and extra contribution, proportional to the strain rate, is added to the stress
tensor, as described in Viscous Damping.
You can specify the viscous damping for volumetric strains and shear strains
independently.
dK = --------- = ----
2f
If, on the other hand, you would want to use a viscous damping, corresponding to a
certain Rayleigh stiffness damping, the conversion to bulk and shear viscosity can be
made using the expressions
b = K dK
v = G dK
where K and G are the bulk and shear moduli, respectively. Equivalently, you can
transform between loss factor damping and viscous damping,
Viscoelastic Damping
In some cases, damping is included implicitly in the material model. This is the case for
a viscoelastic material, where damping operates on the shear components of stress and
strain.
When viscoelasticity is modeled in the frequency domain, it will act as a loss factor
damping. The complex modulus G*() is the frequency-domain representation of the
stress relaxation function of viscoelastic material. It is defined as
G = G + jG = 1 + j s G
where G' is the storage modulus, G'' is the loss modulus, and their ratio sG''G' is
the loss factor. The term G' defines the amount of stored energy for the applied strain,
whereas G'' defines the amount of energy dissipated as heat; G', G'', and s can all be
frequency dependent.
Piezoelectric Damping
Piezoelectric losses are more complex and include coupling and electrical losses in
addition to the material terms. For damping in piezoelectric materials, see Piezoelectric
Losses.
For more details, see the section Modal Solver in the COMSOL
Multiphysics Reference Manual.
• Thin structures, where the deflection is of the same order of magnitude as the
thickness.
• Where the structure exhibits large rotations. A rigid body rotation of only a few
degrees causes significant strains and stresses in a material where a linear strain
representation is used.
• Where the strains are larger than a few percent.
• Contact problems.
• Where a prestress must be taken into account for computing the dynamic response
of a structure.
• Buckling problems.
• Where a deformed mesh is used.
• Fluid-structure interaction problems.
• Contact Modeling
• Fluid-Structure Interaction
• The more formal theory is described in Analysis of Deformation
Figure 2-14 shows that as the beam deforms, the shape, orientation, and position of
the element at its tip changes significantly. Each edge of the infinitesimal cube
undergoes both a change in length and a rotation that depends on position.
Additionally, the three edges of the cube are no longer orthogonal in the deformed
See Frames and Coordinate Systems in the theory chapter for more
details.
An even more fundamental quantity is the deformation gradient, which contains the
derivatives of the deformed coordinates with respect to the original coordinates:
x
F = ------- = u + I
X
The deformation gradient contains all information about the local deformation in the
solid, and can be used to form many other strain quantities. As an example, the Green–
Lagrange strain is
1 T
= --- F F – I
2
An element at a point (X, Y, Z) specified in the material frame moves with a single
piece of material throughout a solid mechanics simulation. It is often convenient to
define material properties in the material frame as these properties move and rotate
naturally together with the volume element at the point at which they are defined as
the simulation progresses. In Figure 2-14 this point is illustrated by the small cube
highlighted at the end of the beam, which is stretched, translated, and rotated as the
beam deforms. The three mutually perpendicular faces of the cube in the Lagrange
Figure 2-14: An example of the deformation of a beam showing the undeformed state
(left) and the deformed state (right) of the beam itself with an element near its tip
highlighted (top), of the element (center), and of lines parallel to the x-axis in the
undeformed state (bottom).
• A Solid Mechanics interface and a moving mesh feature (for example, Deforming
Domain or Rotating Domain) have a common selection. In this case, the selection in
When two physics interfaces are competing for frame control, you will get the error
message “Multiple moving frame specifications on the same selection” when trying to
run the study. To identify the problem, go to the settings for the study step, and select
Modify model configuration for study step. There, you will get an overview of the frame
control in the model. Note that it is quite possible that several features control the
spatial frame, as long as it is on different geometric selections.
When you select a physics interface in the tree view, you can click the Control Frame
Deformation button ( ) to toggle whether that interface should control the spatial
frame or not.
For the Cauchy stress tensor, both the force components and the normal to the area
on which the force is acting have fixed directions in space. This means that if a stressed
body is subjected to a pure rotation, the actual values of the stress components will
change. What was originally a uniaxial stress state might be transformed into a full
tensor with both normal and shear stress components. In many cases, this is neither
what you want to use nor what you would expect.
Consider for example an orthotropic material with fibers having a certain orientation.
It is much more plausible that you want to see the stress in the fiber direction, even if
the component is rotated. The Second Piola–Kirchhoff stress has this property as it is
defined along the material directions. In the figure below, an originally straight
cantilever beam has been subjected to bending by a pure moment at the tip. The
xx-component of the Cauchy stress and Second Piola–Kirchhoff stress are shown.
Since the stress is physically directed along the beam, the xx-component of the Cauchy
stress (which is related to the global x direction) decreases with the deflection. The
Second Piola–Kirchhoff stress, however, has the same through-thickness distribution
all along the beam, even in the deformed configuration.
Unfortunately, even without a rotation, the actual values of all these stress
representations are not the same. All of them scale differently with respect to local
volume changes and stretches. This is illustrated in the graph below. The
xx-component of four different stress measures are plotted at the fixed end of the beam
from the example above. At this point, the beam axis coincides with the x-axis, so the
directions of all stress tensor components coincide. In the center of the beam, where
strains, and thereby volume changes are small, all values approach each other. For a
case with large rotation but small strains, the different stress representations can be
seen as pure rotations of the same stress tensor.
If you want to compute the resulting force or moment on a certain boundary based on
the stresses, there are in practice only two possible choices: Either integrate the Cauchy
stress over the deformed boundary, or integrate the First Piola–Kirchhoff stress over
the same boundary in the undeformed configuration. In COMSOL Multiphysics this
corresponds to selecting either Spatial frame or Material frame in the settings for the
integration operator.
ALE METHOD
In the case of solid mechanics, the material and spatial frames are associated directly
with the Lagrangian and Eulerian frames, respectively. In a more general case (for
example, when tracking the deformation of a fluid, such as a volume of air surrounding
a moving structure) tying the Lagrangian frame to the material frame becomes less
desirable. Fluid must be able to flow both into and out of the computational domain,
without taking the mesh with it. The arbitrary Lagrangian-Eulerian (ALE) method
allows the material frame to be defined with a more general mapping to the spatial or
Eulerian frame. In COMSOL Multiphysics, a separate equation is solved to produce
this mapping — defined by the mesh smoothing method (Laplacian, Winslow,
hyperelastic, or Yeoh) with boundary conditions that determine the limits of
If any active feature in the model requires the analysis to be geometrically nonlinear,
the Include geometric nonlinearity check box is selected automatically, and it cannot be
cleared. The following physics features force this behavior:
• Contact, because the deformation between the contacting boundaries must be part
of the solution.
• Moving mesh (when at least one deforming domain is active).
• Large strain plasticity.
• Hyperelastic materials, which are always formulated for large strains.
Usually you would also want to use geometric nonlinearity when a Moving Mesh
interface is present, but this setting is not forced by the program.
When you select a geometrically nonlinear study step, the behavior of several features
differs from that in a geometrically linear case:
The geometric nonlinearity in the Beam interface is limited to large rotations and
displacements, but small strains.
The effect of using geometric nonlinearity in these interfaces is limited to the stress and
strain representation as the other effects described in Geometric Nonlinearity for the
Solid Mechanics Interface are not relevant.
Prestressed Structures
You can analyze eigenfrequency, frequency domain, or time-dependent problems
where the dynamic properties of the structure are affected by a preload, such as a
tensioned string.
Usually, a study of a prestressed problem includes using two study steps. The first step
is a Stationary step in which the static preload is applied. The preload step can be
computed with or without taking geometric nonlinearity into account. In the second
study step, where you compute the eigenfrequency or the frequency response, it is
necessary to take geometric nonlinearity into account. Even if the displacements and
strains are small, this is what gives the prestress contribution to the equations.
There are three predefined study sequences for prestressed dynamic analysis:
• Eigenfrequency, Prestressed
• Frequency Domain, Prestressed
• Frequency Domain, Prestressed, Modal
The prestressed study types assume that the loading causes small perturbations around
the prestressed state.
The same principles apply also to a linear buckling analysis, except that both study steps
should be geometrically linear. The nonlinear contribution is included in the
formulation of the buckling eigenvalue itself.
FOLLOWER LOADS
Loads that change orientation with deformation, such as a pressure, actually contribute
not only to the load but also to the stiffness. This is a physical effect and not just a
numerical artifact. Whether such loads are included or not in an Eigenfrequency study
step will affect the computed eigenfrequencies. If you for some reason do not want this
effect, you must suppress the load in the Physics and Variables section of the
Eigenfrequency node.
If you use a local coordinate system for describing a load, you must, in
case of geometric nonlinearity, pay attention to whether that coordinate
system has constant axis orientations or not. As an example, the default
boundary system has Frame set to Deformed Configuration, so that a load
represented in that system will behave as a follower load. Change to
Reference Configuration if the load should act in fixed directions.
There are three Preset study types which can be used to set up these two
study steps: Eigenfrequency, Prestressed; Frequency Domain, Prestressed;
and Linear Buckling.
If you want to explicitly prescribe the stress field for a prestressed analysis
rather than solving for it, you should not use the two study step
procedure. In such a case, prescribe the stress field using an Initial Stress
and Strain, External Stress, or External Strain node. Then add a separate
Eigenfrequency or Frequency Domain study and select Include Geometric
Nonlinearity in the settings for the study step.
1 u i u j u k u k
ij = --- -------- + -------- +
2 X j X i --------
- ---------
X i X j
(2-27)
k
The first two items above result in another set of constitutive equations for large
deformation piezoelectricity:
T
S = cE – e Em
P m = e + 0 rS – I E m
= rS – I
–1
D m = P m + 0 JC E m
T
C = F F
Fields in the global orientation result from the following transformation rules:
–T
E = F Em
–1
P = J FP m
(2-28)
–1
D = J FD m
–1
v = V J
where F is the deformation gradient; J is the determinant of F; and v and V are the
volume charge density in spatial and material coordinates, respectively. The
deformation gradient is defined as the gradient of the present position of a material
point x X + u:
x
F = -------
X
T
S = cE – e Em
–1
D m = e + 0 rS E m + 0 JC – I E m
The electrical part separates into two different cases: For solid domains, the electric
energy consists of polarization energy within the solid and the electric energy of free
space occupied by the deformed solid — the same as for the piezoelectric materials.
For nonsolid domains this separation does not occur, and the electric displacement in
D = 0 r E
The use of ALE has impacts on the formulation of the electrical large deformation
equations:
• The first impact is that with ALE, the gradient of the electric potential directly
results in the electric field in the global orientation, and the material electric field
results after transformation.
• The most visible impact is on the boundary conditions. With ALE, any surface
charge density or electric displacement is defined per the present deformed
boundary area, whereas for the case without ALE, they are defined per the
undeformed reference area.
References
1. J. Yang, An Introduction to the Theory of Piezoelectricity, Springer Science and
Business Media, N.Y., 2005.
2. A.F. Bower, Applied Mechanics of Solids, CRC Press, Boca Raton, FL (http://
www.solidmechanics.org), 2010.
In a contact analysis, you solve for the contact state and the contact forces. The most
fundamental quantity in a contact analysis is the distance between the objects which
may come into contact, the gap. If the gap is positive then there is no interaction. The
task of the contact algorithm is to ensure that the gap never becomes negative; that is,
to avoid overclosure. In order to accomplish this task, contact forces must be
introduced.
For normal contact, the state only consists of being in contact or not, and the force
variable is the contact pressure in the normal direction. With friction included, there
are two possible states for the relative tangential displacement when in contact:
sticking or sliding. The tangential force vector is added as force variable.
In this section:
• In the finalization step of the geometry sequence, you should normally have Action
set to Form an assembly. If Form a union is used, then the contacting boundaries must
be geometrically separated in the initial configuration.
• Add Contact Pair nodes under Definitions. A contact pair consists of two sets of
boundaries, which are called the source and destination boundaries. Contact pairs
can also be added automatically, based on boundary adjacency when the Form an
assembly action is used. The geometric gap distance is a variable set up by the
contact pair, which also define operators for mapping variables or expressions
between the selected boundaries.
• Use the default Contact node or add new Contact nodes in the physics interface. In
the Contact node, you select the contact pairs to be used, and provide the settings
for the physical and numerical properties of the contact model.
• If relevant, add Friction, Slip Velocity, Adhesion, Decohesion, or Wear subnodes to
Contact.
In a multiphysics analysis, a contact problem can also incorporate for example changes
in the heat flux or electric current through the contacting boundaries. You will then
also need to add corresponding features in the other participating interfaces, like a
Thermal Contact node in the Heat Transfer in Solids interface. The contact state and
contact pressure used by other physics interfaces are always supplied by the structural
mechanics interface.
The fact that you add a Contact node to your model will automatically make all study
steps geometrically nonlinear. For the default Contact node, this requirement can be
removed by selecting Disconnect pair.
• In many cases, only the normal forces are significant for the general force
distribution in the structure, while the frictional forces only cause minor local
effects.
• The values of the friction coefficients are difficult to obtain, and unless the structure
is assembled under well controlled conditions, the magnitude of the friction can
vary a lot. If there is a large degree in uncertainty in the input data, it can be argued
that the inclusion of friction does not add much to the quality of the results.
• Adding friction to a contact problem will often increase the computation time
significantly, or even cause convergence problems.
There are a number of situations when friction modeling cannot be avoided. Some of
them are:
• When a significant portion of the load is carried by shear tractions acting on the
contact boundaries.
• When a tangential force is necessary to maintain stability and to avoid rigid body
motions. In many cases, it is however possible to replace the friction by a suitable
boundary condition instead, as long as there are no resultant forces being resisted
by such a constraint.
• When modeling wear.
• When the frictional dissipation is an important component of a dynamics problem,
or when it is needed as a heat source in a thermal analysis.
In some cases, such as when a brake pad slides along a brake disc, the size and
orientation of the slip velocity is known. You can then employ a simplified form of
friction modeling by assuming the tangential contact to always be in a slip state, which
simplifies the computation of the friction forces. This is done using the Slip Velocity
node. This is particularly useful for wear modeling.
If adhesion is active between the contact boundaries, it is possible to break the bond
by adding a decohesion model. You can choose between several different decohesion
models.
Adhesion and friction can be combined, but during the time that two boundaries are
bonded to each other through adhesion, any settings for friction are ignored.
INCLUDING WEAR
It is possible to model adhesive or abrasive wear of the material when the contacting
boundaries are sliding along each other. The removal of material from the domain
adjacent to the contact boundary can be modeled using three, fundamentally different,
approaches.
For the Solid Mechanics and the Multibody Dynamics interfaces, the most general
approach relies on the deformed geometry concept, where material is actually removed
by using an adaptive mesh technique. In the second, simplified, approach wear is
incorporated as an offset in the contact condition. This approach is computationally
less expensive, and is suitable as long as only small amounts of material are removed,
and wear does not change the orientation of the normal to the boundary significantly.
In the Shell and Membrane interfaces the structural domain is represented by only a
meshed surface. Therefore, a different, more suitable approach is used, in which the
thickness variable and the midsurface offset to the meshed boundary are modified.
Computing the amount of wear involves solving a rate equation, hence, it is only
possible to compute wear in time-dependent studies. The wear rate is typically a
function of the contact pressure and the relative slip velocity between the contacting
boundaries.
Penalty Method
The default penalty method is rather simple and robust method to introduce the
contact conditions. Roughly speaking, it is based on inserting a stiff distributed spring,
active only in compression, between the contacting boundaries. In addition to the
robustness, it has the advantage that no special solver is required, which makes it easier
to set up multiphysics problems and time-dependent studies. However, the penalty
method only enforces the contact conditions approximately. The contact forces
computed are thus less accurate than when using, for example, the augmented
Lagrangian method, and there is always some overclosure between the contacting
surfaces. For stick friction, there is also an ‘elastic’ deformation due to the penalization
of the constraint.
When using the penalty method, there is always a tradeoff between accuracy and
stability. While a large penalty factor will reduce nonphysical overclosures, the problem
may become ill-conditioned and unstable if it is too large. It might therefore be
beneficial to accept some penetration between the contacting objects. Note, however,
that if the penalty factor is too small, the contact condition may be violated.
Nitsche Method
The Nitsche method for contact is an extension of the general method suggested by
J. Nitsche in 1971 to weakly impose Dirichlet conditions. Conceptually, it can be
viewed as an enhancement of the penalty method where the surface traction of the
adjacent domains is used to improve the accuracy of the contact condition. Also, as for
the penalty method, it has the advantage that no special solver is required, which makes
it easier to set up multiphysics problems and time-dependent studies.
For the first two problems in particular, the accuracy and stability can be further
improved by increasing the quadrature order for the contact equations, see also
Quadrature Settings. Hence, the Nitsche method uses a higher quadrature order by
default than the other contact methods.
• Symmetric
• Skew-symmetric
• Nonsymmetric
The default nonsymmetric formulation is recommended for the majority of cases and
provides a good tradeoff between performance and stability. In cases where the default
formulation fails, the skew-symmetric formulation can be an alternative, but it includes
additional equations that are expensive to evaluate that can increase the solution time.
CONTACT DETECTION
The contact search is made only towards one side of the source and destination
boundaries, determined by the positive direction of the contact normal of the selected
boundary. For contact to be detected, this means that the source and destination
contact normals must point towards each other. Generally, the contact normal
coincides with the geometry normal that is determined by the direction of the
boundary. However, there are situations where this is not the case, or where the
Contact node can control the sign of the contact normal. Some common situations are
summarized in the following:
• If the contact boundary is on the exterior of a domain, then the geometry normal
always point outwards from the domain. This is the most common case, and then
no other considerations are needed.
• In some cases, like fluid-structure interaction, there may be a domain with a moving
mesh between the two contacting objects. The source and destination boundaries
will then, in the geometrical sense, be interior boundaries. In this case, the physics
interface (Solid Mechanics or Multibody Dynamics) defines a normal vector which
is pointing outward from the boundaries that are external to the physics interface.
Thus, this case is also handled automatically.
• When a boundary without an adjacent domain is selected, you need to be careful so
that the normal is pointing in the intended direction.
• When using the Shell, Layered Shell or Membrane interfaces, contact can potentially
occur on both sides of the boundary. In a single Contact node, you can only model
contact on one side. In the Contact Surface section, you can select whether the
contact should occur on the top side or bottom side. The ‘top’ and ‘bottom’ sides
are defined by the orientation of the physics interface normal, which may differ from
the geometry normal. In these interfaces, the orientation of the physics normal is
controlled by the Boundary System that is attached to each boundary through the
The actual normal vector used in the contact search algorithm can be
visualized by plotting <pair_tag>.n<coord_label>. For example, the
variable p1.nx gives the x-component of the spatial normal used by
contact pair p1. Plotting the contact normal vector can be useful to verify
that the source and destination normals point towards each other.
In 3D, the rules for the orientation of a boundary are more complicated. In general,
you have to visualize the normals, in order to check its orientation. This can be done
in different ways:
• For a mesh without physics, you will need to use the general result presentation
tools. Follow these steps:
- Run Get Initial Values for an arbitrary study in order to create data for result
presentation.
- Add a 3D Plot Group under Results.
- Add an Arrow Surface plot to the new plot group.
- In the Replace Expression ( ) dialog, select the geometry normal.
If the source boundary is not part of the a physics interface with a Contact node, the
gap is computed using only the current location of its mesh, ignoring any physical
properties that may exist there. In this case, the Contact node has no control of the
contact normal used by the source boundary. Care must therefore be taken to ensure
There are two main scenarios where you may want to use a source selection that is not
in the current physics interface:
Contact Pairs
• If there is a significant difference in the stiffness in the normal direction between the
source and destination boundaries, select the stiffer side as source. This is especially
important if the difference in stiffness is quite large (for example, more than a factor
of ten).
• When the contacting parts have approximately the same stiffness, consider the
geometry of the boundaries instead. Make a concave boundary the source and a
convex boundary the destination rather than the opposite.
• If one of the boundaries belongs to a part that is rigid it should be selected as the
source boundary. Rigid boundaries can be created in several ways, for example
- The underlying domain has the Rigid Material material model.
- The boundary or its underlying domain is constrained by, for example, Fixed
Constraint or Prescribed Displacement.
- The source boundary is meshed, but has no physics attached.
- The source boundary position is controlled by Prescribed Deformation under
Moving Mesh.
• If only one side of the contact pair is within a physics interface that has a Contact
node, that side must be the destination side.
If it is difficult to follow the above guidelines, the same boundaries can be selected as
both source and destination. Doing so results in an unbiased (or symmetric)
formulation that is less sensitive to, for example, difference in stiffness or mesh density
between the contacting boundaries. However, such a formulation involves evaluating
equation and contact mapping at additional points and can thus be more expensive to
use.
A contact pair is just a definition, and does not have to be referenced by any Contact
node.
Neither the source, nor the destination, in a contact pair has to be geometrically
contiguous. In practice, this means that you often only need a few contact pairs in a
model. The number of pairs actually needed will be determined by how many different
settings that are required in the Contact Pair and Contact nodes.
If you have many contact pairs in your model, it is a good idea to manually set the Label
of each pair in order to simplify the identification during subsequent selections in the
Contact nodes.
The automatic pair generation will not know which side to use as source
or destination. Based on the suggestions in Selecting Source and
Destination above, you may need to switch selections using the Swap
Source and Destination ( ) button in the Source Boundaries section of the
contact pair settings.
In such cases, it is possible to simplify the problem by selecting the Mapping Method to
be Initial Configuration in the Contact Pair node. With this setting, a material point on
the destination boundary will see the same material point on the source boundary
throughout the entire simulation; that is, the mapping is constant. This setting will
make the contact analysis run faster and convergence to be more stable.
The analysis is geometrically nonlinear also when using this option, and the contact
region can still have arbitrarily large displacements and rotations.
Friction can be modeled. Even though there is no sliding in a geometrical sense, the
difference in tangential displacement is computed.
You cannot mix contact pairs with the two different types of Mapping
Method within the same Contact node.
It is always important that the geometry is well resolved, so that a curved contact
boundary actually is seen as curved rather than “faceted”. The density of the mesh
often needs to be finer than what would be needed to resolve stresses on a similar
boundary without the contact conditions. If the normal to the contact boundary
changes much from one element to the next, there is a risk that the contact analysis
does not converge.
If the source boundary is rigid, there are no requirements on the mesh size of the
destination boundary. In this case, you may use a significantly finer mesh on the source
boundary than on the destination boundary. This is sometimes necessary in order to
resolve the geometry well. On the other hand, if you have a flat rigid boundary, you
only need to mesh it with a single element.
Penalty Factor
An important parameter in the settings for the Contact node is the penalty factor.
When running into convergence problems, check the penalty factor settings and
consider changing the current value. It is used by all contact methods, but its
interpretation differs:
• In the penalty method, the penalty factor has a straightforward interpretation as the
stiffness of a distributed spring connecting the two contacting boundaries. A higher
value will decrease the unphysical penetration, but can be detrimental to the
numerical conditioning of the stiffness matrix. A too small value can, however, result
in violation of the contact condition.
• In the augmented Lagrangian method, the penalty factor is a numerical parameter,
which affects the convergence properties of the algorithm. Its interpretation is
different depending on the chosen solution method. For a segregated method, the
penalty factor mainly affects the rate of convergence of the outer iteration loop; a
higher value typically leads to faster convergence, but can decrease the stability of
the algorithm for certain problems. For the coupled method, the penalty factor
mainly affects the structure of the underlying equations, and the solution is typically
In the augmented Lagrangian method and the Nitsche method, the value
of the penalty factor does not affect the accuracy of the final solution, like
it does in the penalty method.
High values of the penalty factor can lead to an ill-conditioned stiffness matrix and
convergence problems in the Newton iterations. This is often identified by that the
damping factor reported by the solver becomes less than 1 for many Newton iterations,
or by that the structure “jumps” into an unphysical state. Large errors returned from
the linear solver are also an indication that the stiffness matrix is ill-conditioned. If this
occurs, decrease the penalty factor.
The default value for the penalty factor is based on a characteristic stiffness. The
default is an “equivalent” Young’s modulus (Eequ) of the material on the destination
side. For linear elastic isotropic materials, Eequ is the actual Young’s modulus. For
other types of materials, Eequ is defined by an estimate of a similar stiffness at zero
strain. For materials that are user defined or in other ways nonstandard (for example,
anisotropic with large differences in stiffness in different directions), Eequ might need
to be replaced with another estimate.
In a situation where the contact is well established, using relaxation will however cost
extra iterations, and it may even lead to a loss of convergence.
For this method, the penalty factor can be tuned in several ways. You have three basic
choices, ranging from simple to advanced:
• With a Preset penalty factor, you can choose having it tuned for Stability, Speed, or
Bending. With Stability, relaxation is used in every step. With Speed, a constant
penalty factor is used, and the value that is used is equal to the final value obtained
when using Stability. For Bending, a constant and low penalty factor is used in all
iterations. The value corresponds to the initial value when using Stability. As the
name implies, this option is intended for bending dominated problems where the
structural stiffness is much lower than the bulk stiffness of the material.
• With Manual tuning, you can adjust the magnitude of the penalty factor, and also the
relaxation strategy.
• With User defined, you can enter any expression for the penalty factor.
For details about these settings, see the documentation of the settings for
the Contact node.
Some hints for selecting the penalty factor for the segregated augmented Lagrangian
method:
• Use relaxation only when large changes in the contact state is expected.
• If an analysis takes a long time, check the convergence graphs. If the contact
variables show a steady, but slow, convergence it may help to increase the penalty
factor. Increase by a factor of 2 – 5.
• If a model fails to converge, and the contacting parts appear to have moved far from
each other, it is probable that the penalty factor is too high. You can then either
decrease the total stiffness or add more relaxation.
Trigger Cutback
If, during the iterations, a contact problem comes into a state where it is far from the
converged solution, it is unlikely that the solution will ever converge. In such a case,
much computing time can be spent before the maximum number of iterations is
reached, and the solver makes an attempt with a smaller time or parameter step. To
• When analyzing problems with for example shrink fit, nominal dimensions can be
used for the geometry, and the overclosure included in the gap by using a positive
offset.
• When there is a small clearance between two boundaries, a negative offset can be
used instead of changing the geometry.
• If geometries are such that a large overclosure exists in the initial configuration, the
contact algorithms may not converge. You can then add a negative offset, which is
slowly removed by letting it depend on time or on the parameter in the parametric
continuation solver.
• A positive offset can be used to avoid a complete collapse of a mesh that exists
between the source and destination boundaries. This is discussed in more detail in
Multiphysics Contact Analysis.
When the source and destination boundaries are curved, the discretization introduced
by the meshing may lead to small variations in the computed distance between the
source and destination boundaries, even though the geometrical shapes of the two
objects are ideal. When modeling for example a shrinkage fit, this effect can cause
significant fluctuations in the contact pressure. If you select Force zero initial gap, the
computed distance from destination to source will be adjusted by the initial gap
distance detected by the contact search. Positive gap distances smaller than the
tolerance gap are adjusted to be zero. By default, gap is set to Inf, which means that
all gaps and overclosures detected are adjusted to be zero. This adjustment can be
combined with an offset. The offset is applied to the adjusted gap value.
You can only apply an offset to the source boundaries if they belong to the same physics
interface as the destination boundaries.
Initial Value
In the augmented Lagrangian method, where the contact pressure is a dependent
variable, it can be given an initial value. In force-controlled contact problems where no
other stiffness than the contact prohibits the deformation, the initial contact pressure
is crucial for convergence. If it is too low, the parts might pass through each other in
the first iteration. If it is too high, they will never come into contact.
Discretization
When using the augmented Lagrangian method it is possible to change the type and
order of the shape functions used for the contact pressure and friction force degrees of
freedom. The default is linear shape functions, which ensures that there are no
over-constraints in the contact interface. It is allowed to use a shape function order that
is equal to or lower than the order of the displacement field. Increasing the order of
the contact variables from the default can increase the accuracy of how well the contact
conditions are enforced, but can impair convergence and increase the computational
cost.
For a discretization other than Linear, the lumped solver is no longer optimal for the
contact pressure update when using a segregated solution method. In such cases, a
standard segregated step should be used. The default solver generation takes this into
account, but if you later modify the discretization, you should update the solver
sequence.
Quadrature Settings
The weak equations set up by the Contact node and its subnodes typically involve
discontinuous functions. These originate from the contact mapping, where the source
and destination meshes, in most cases, are nonconforming. The default quadrature
used in the numerical integration of these integrals is equal to the order used by the
displacement field. For a quadratic displacement field, this means integration order
equal to four.
In most situations, the default quadrature provides sufficient accuracy of the numerical
integration. However, it is sometimes necessary to increase the integration order for
the contact weak equations. This can improve the stability of the contact model since
Jacobian Contribution
In the Advanced section of the Contact node, there is an option to specify the type of
Jacobian contribution from the contact equations. The default Automatic option will
choose a suitable setting depending on the mapping method used by the contact pair.
However, if controlled manually the Nonsymmetric option is the preferable choice
especially when the source boundaries undergoes large deformations, since it is more
robust. The Symmetric option can be attractive for large models since it preserves the
symmetry of the global stiffness matrix, as long as no other features cause it to be
nonsymmetric. This can decrease the solution time and memory requirements when
solving the model.
Friction Parameters
The friction model specifies the threshold for the friction force in the contact pair. If
the computed (trial) friction force is above this threshold value, the contact is in a stick
state; if the (trial) friction force is above the threshold, the contact is in a slip (or
sliding) state.
Two predefined friction models are available based on the classical Coulomb law,
where the friction force is proportional to the contact pressure through the friction
coefficient. Both Coulomb laws are additionally generalized by allowing specification
of minimum (cohesion) sliding resistance and a cap that sets the maximum tangential
traction.
It is also possible to define the threshold for the friction force as an arbitrary expression
that may depend on any quantity in the model, for example temperature or position.
The only limitation is that the expression may not implicitly depend on the current
value of the friction force.
This section is not available when the Nitsche method is used. The same penalty factor
as set in the parent Contact node will be used for friction as well.
Knowing the slip velocity greatly simplifies the computation of friction forces. There is
no need to determine the transition between stick and slip contact, which can be
difficult.
The adhesive layer is conceptually without thickness, but by specifying an offset in the
Contact node, you can to some extent include the dimensions of the adhesive layer.
The adhesive layer always has a finite stiffness. For a glue layer, this represents the true
stiffness. For a more conceptual joining of two boundaries, this stiffness should be
considered in the same way as the penalty stiffness in the contact formulation. The
stiffness can differ between tension and compression: In compression the stiffness is
always taken as the penalty stiffness, whereas you can change the tensile stiffness.
Two alternative CZM are available. The Displacement-based damage models defines
damage growth as a function of a mixed mode displacement quantity. It comes with
several traction separation laws that associate the onset of damage with the peak
strength of the interface. For some of them, it is possible to choose between different
mixed mode failure criteria. The Energy-based damage models define damage growth
as a function of the stored undamaged elastic energy density of the interface. It also
comes with several different traction separation laws. However, these are more general
and define the onset of damage at an arbitrary elastic energy density. In principle, you
can define the model so that damage initiates immediately during loading of the
adhesive layer, that is for zero energy density. The strength of the interface is then
determined by the critical energy release rate and the shape of the damage evolution
function. In this way, the energy-based damage models can be viewed as a
regularization of linear elastic fracture mechanics.
The most general technique to model the removal of material during the wear process
relies on the deformed geometry concept. When selecting the Deformed geometry
formulation, the wear feature adds a (hidden) Deforming Domain feature that controls
the material frame through an adaptive mesh smoothing. The wear, as computed in the
Wear node is fed as a (hidden) Prescribed Normal Mesh Displacement boundary
condition to the deforming domain, and thus describes the actual removal of material
from the geometry. When using this formulation, you must be aware that the adaptive
mesh means that state variables stored in Gauss points, for example plastic strains or
creep strains, will not represent the same material points all the time. Whether or not
this effect is acceptable must be judged on a case-by-case basis. Unless the amount of
material that is removed is large, or gradients are strong, this is mainly an issue close
to the boundary where material is removed by the wear process.
The slip velocity used for the wear computation can be obtained from either a Friction
node or a Slip Velocity node, so one of these two nodes should be present and active
under the same Contact parent node. For the Generalized Archard wear model, this is a
In general, modeling wear on the destination side is slightly more accurate, since it is
there that the contact pressure and slip velocity are originally computed. When
modeling wear on the source side, these quantities are mapped from the destination
boundary. Multiple Wear nodes under the same Contact node contributes with each
other, which means that it is possible to model wear on both source and destination
simultaneously. However, adding multiple wear contributions to either source or
destination may give unphysical results. On the source side it is possible to also use the
Rigid Material material model; this is not permitted on the destination due to general
restrictions of the Contact node.
• Create the geometry or set initial values for the displacement variables so that there
is a small penetration in the initial configuration.
• Use boundary conditions giving a prescribed displacement rather than a prescribed
force. When possible, this is usually the best way to stabilize problem. Note that you
can always obtain the force actually applied from the reaction forces.
and similarly in any other direction that needs to be restrained. It is not important
whether the spring is applied to domains, boundaries, or edges, but it should not
create significant local forces. The value for the stiffness k should be chosen so that
the force generated by the spring balances the external load at a sufficiently small
displacement. A too weak spring will give a too large initial overclosure of the
contact boundaries. A too stiff spring might influence the solution too much.
To model dynamic contact events, two specialized contact methods are provided, the
Penalty, dynamic and the Augmented Lagrangian, dynamic methods. Both are based on
a viscous formulation that constrains the gap rate to be zero, ensuring that the normal
contact is dissipative and does not introduce any spurious energy contribution to the
system. Since the methods are dissipative, they are mainly intended for short duration
events, such as soft impact between two bodies. For prolonged interaction between
two bodies, energy dissipation can become significant, and overclosures can become
large, since the gap rate is only approximately zero. Both the dissipation and the
accuracy are controlled by a penalty factor that for these two methods conceptually
represents a dashpot, rather than a spring. It therefore has a characteristic time user
input that sets its magnitude. As a rule-of-thumb, it should be of the order of the
contact event duration, but the best choice must be decided on a case-by-case basis.
The Penalty, dynamic method also provides the possibility to combine the stiffness and
viscous based penalization of the normal contact. For impact analysis, it is often best
to use only the viscous formulation by setting the stiffness Penalty factor control to
Viscous only.
When modeling dynamic contact, the main interest is often the kinematics between the
contacting bodies. If you rely on the (default) adaptive time-stepping algorithm, the
Regardless of the method that is used, and how the solver is set up, it is good practice
to do an a posteriori check of conservation of momentum and energy, to ensure that
the solution is acceptable.
This class of problems often exhibit a high degree of nonlinearity, which may lead to
convergence problems in the nonlinear solver. As an example, consider heat transfer
through the contact area, where initially only a small spot is in contact. The solution
for the temperature field is then extremely sensitive to the size of the contact area. If,
at the same time, the solid deforms due to thermal expansion, there may be large
changes in the contact area between iterations,
It is important to resolve the size of the contact area accurately, that is, to
use a very fine mesh in the contact area when modeling fully coupled
multiphysics problems.
If the contact area is larger, a fine mesh is not required because then the temperature
solution is not that sensitive to the size. If possible, start with an initial configuration
where the contact area is not very small.
You can use the contact variables (gap and contact pressure) in expressions for
quantities in other physics interfaces. As an example, a thermal resistance in the contact
region can depend on the contact pressure.
This is the case in, for example, fluid-structure interaction (FSI) problems. Here, the
equations in the fluid are solved on a domain with a moving mesh, so that the shape
of the fluid domain is controlled by the displacements of the solid. Another case of the
same type is when there is an electric field in an air gap.
If contact is established, the mesh in the original gap between the source and
destination boundaries will collapse. This must be avoided. The remedy is to add an
offset in the contact settings to either the source boundary, the destination boundary,
or both. If you do this, contact forces will be transmitted without the geometrical gap
being fully closed.
INTERFERENCE FIT
Interference fits can be analyzed using contact modeling. This is necessary if you are
interested in checking the connection with respect to, for example, slipping or local
stresses.
There are two possible approaches for modeling interference fits, both equally valid:
• The actual geometry of the two parts before assembly is modeled. In this case, there
will be some overlap between the two domains.
• A nominal geometry in which the contacting boundaries have the same location is
created. In this case, the overlap is described as part of the contact modeling.
An imported CAD geometry can use either of these approaches, depending of the
strategy used during the geometry creation. Often, the geometrical parts are modeled
as nominal, and instead equipped by tolerance information that describe the amount
of interference.
True Geometry
With a true geometry, you can often immediately solve the contact problem.
Sometimes convergence problems may, however, appear, in particular if the material
model is nonlinear. The cause is often that the initial overlap is too far from the final
solution.
To deal with such a problem, add an offset in the in the settings for the Contact node.
The offset should be defined by a parameter, so that the boundaries of the two domains
are barely in contact in the initial state. Now, the offset can be reduced to zero
step-by-step, using an auxiliary sweep in the solver.
SELF-CONTACT
To model self-contact, include the same boundaries in both the source and destination
selections of the Contact Pair definition. This will cause the boundaries to act as both
source and destination in the contact search and mapping. For mechanical contact, this
results in a unbiased (or symmetric) contact formulation, as the contact conditions are
formulated on both sides of the contact pair. Note that a source is not allowed to
partially intersect the destination when used for mechanical contact.
This technique to model self-contact means that some of the considerations discussed
in this chapter regarding contact modeling do not apply. For example, instead of the
recommendations in Meshing for Contact Analysis, it is recommended to use a
uniform mesh element size along the contacting boundaries. Self-contact is a case
where it might be necessary to increase the quadrature order used in the weak
equations, see Quadrature Settings.
In the Solid Mechanics interface, the following boundary conditions are by default
recognized as fallback conditions if their selection intersects any source and destination
boundary:
• Boundary Load
• Spring Foundation
• Added Mass
• Free
In the Shell and Membrane interfaces, the following boundary conditions can be
recognized as fallback conditions:
• Face Load
• Spring Foundation
• Added Mass
To enable the fallback condition in these nodes, change the Allowed region to Fallback
and nonpair regions in the Applicable Pair Regions section of settings window of the
node. To display the Applicable Pair Regions section, click the Show More Options
button ( ) and select Advanced Physics Options in the Show More Options dialog box.
The most common case when a fallback condition is used is when there is a pressure
load acting on the part of the boundary that is currently not in contact. In this case
you would add a Boundary Load or Face Load with a selection that intersects with the
source and destination boundaries.
GENERAL SETTINGS
The following solver settings can help to successfully perform contact simulations in
general:
• The convergence check relies on the scaling of the degrees of freedom, but since
contact pressures and friction forces often are zero over parts of the simulation, you
should not rely on automatic scaling. When the solver sequence is first created, both
contact pressure and friction forces are given a manual scaling which is relevant for
typical metal-to-metal contact. You should in most cases change this to values
appropriate for your application. The variable scaling is accessed under Dependent
Variables in the solver sequence. Set the scale for each variable to a value that is
representative for the expected result. Too large values may give a too early
convergence, while too small values may lead to an excessive number of iterations.
• The default solver sequence generates one lumped step in the segregated solver for
each Contact node. This split of variables into different lumped steps does not
influence the solution as such; you can equally well group the contact variables in a
single lumped step. Each lumped step will however generate an individual curve in
the convergence plot, making it easier to pinpoint the source of possible
convergence problems. You can also increase the granularity even more by changing
Solver log to Detailed in the Advanced node in the solver sequence. This will give a
separate convergence curve for each dependent variable.
The default solver generates a solver sequence that is stable and gives and accurate
solution for the majority of contact problems. However, if convergence is difficult
other settings can be tested:
• The convergence check relies on the scaling of the degrees of freedom, but since
contact pressures and friction forces often are zero over parts of the simulation, you
should not rely on automatic scaling. When the solver sequence is first created, both
contact pressure and friction forces are given a manual scaling typical for
metal-to-metal contact. You should in most cases change this to values appropriate
for your application. The variable scaling is accessed under Dependent Variables in the
solver sequence.
• Since the solution to the augmented Lagrangian can be non smooth, the default
double dogleg nonlinear solver in stationary studies is sometimes too conservative.
The convergence can in such cases often be improved by using a Newton solver, for
example, the Constant (Newton) with a full Jacobian update.
• Although the penalty factor does not affect the accuracy of the solution, it can have
significant influence on the convergence properties of the model.
DYNAMIC CONTACT
The dynamic contact methods in general inherit the properties and suggested solver
settings from the corresponding standard method. With regards to solver settings, the
main difference for the dynamic methods is how to set up the time-dependent solver,
where it is recommended to use manual time stepping, see Dynamic Contact Analysis.
Using the Results while solving functionality in the study step is a good practice. You
can either use a stress plot, or a plot of the contact pressure. In most cases, the scale of
a deformed plot should be set to 1 when monitoring contact problems. If you select
Results while solving in the Segregated node, the plot is updated after each iteration,
thus allowing you to monitor the convergence in detail.
For each contact pair, two global variables that can be used in probe plots are available.
These are the maximum contact pressure (<phys>.Tnmax_<pair>) and the minimum
gap distance (<phys>.gapmin_<pair>).
Looking at the convergence plot will give valuable information about the convergence
properties. There will, as a default, be one graph per Contact node in the Model Tree,
which will help you pinpoint the source of a convergence problem. You can also
increase the granularity even more by changing Solver log to Detailed in the Advanced
node in the solver sequence. This will give a separate convergence curve for each
dependent variable.
You can also select to include information about the contact state in the solver log. To
do that, select the Add contact status to solver log check box in the Advanced section of
the settings for the Contact node. For each contact pair, messages like
will be generated for each time or parameter step. Only changes are reported.
• Variables changed until convergence is reached during the iterations. These variables
appear in the Lumped Step nodes in the Segregated solver or in the Fully Coupled
node.
• Variables used to store the state, once the iterations have converged for a certain
time or parameter step, so called internal degrees of freedom. Such state variables
are not immediately visible in the solver sequence, but it you select the Displacement
Field node under Dependent Variables in the solver sequence, you will see them listed
as ‘Internal variables’.
If you change settings in the Contact or Friction nodes after the solver
sequence has been generated, dependent variables may be added or
removed. The second case is never a problem, but when new dependent
variables are created, they are not automatically added to the groups in the
segregated solver. You may then encounter the error message
“Segregated solver steps do not involve all components.” You
can then either regenerate the solver sequence, or manually insert the
missing variables into the Lumped Step node.
In Table 2-11 the dependent variables that can be created by the Contact and its
subnodes are summarized. To shorten the variable names, the full scope has been
removed. As an example, the contact pressure variable for pair p1 in component comp1,
generated in the Solid Mechanics interface solid, will have the full name similar to
comp1.solid.Tn_p1. In the table, it is shown as Tn.
The Activation subnode can be used for this purpose. You enter an activation
expression to determine whether material is active or not. When this Boolean
expression is satisfied, the material is activated. Rather than truly adding or removing
material, the Activation subnode alters the stiffness and density of the material to
emulate this.
ACTIVATION
Add an Activation subnode when you want to activate or deactivate one or several
domains selected in a Linear Elastic Material node. The Activation expression field is used
to define when material should be activated, and the Activation scale factor is used to
reduce the elastic stiffness and density of the material which is not active.
The activation condition can be any type of expression or function. A common case is
that it is a function of the temperature. The activation expression is evaluated in each
Gauss point. This means that an element can be partially activated. If you want to force
whole elements to be activated, you can for example put the activation expression
inside the centroid() operator.
VARIABLES
Two useful variables are created when you add the Activation subnode. The variable
isactive is set to ‘1’ when the activation condition is satisfied, and it is ‘0’ otherwise.
The variable wasactive is used to record if the material has been active at any previous
step in the analysis. This variable can be used to “lock” the state of activation, once it
has been reached. Suppose that you want the material in the previous example to
remain active even if para later becomes less than 1.5. The activation expression for an
interface with the name solid could then be expressed as:
Similarly, the variable wasinactive is used to record if the material has been inactive
at any previous step in the analysis.
VARIABLE DESCRIPTION
RESULTS
If you have performed an analysis in which only part of the material is active, it is useful
to apply a Filter and only display the regions that actually are active; see Figure 2-18.
When an Activation subnode is added in a Solid Mechanics or Membrane interface,
such a Filter node is automatically added to the default stress plots.
The features are completely analogous, with the difference that a Spring Foundation
connects the structural part on which it is acting to a fixed “ground,” while the Thin
Elastic Layer acts between two parts, either on an interior boundary or on a pair.
A Spring Foundation is most commonly used for simulating boundary conditions with
a certain flexibility, such as the soil surrounding a construction. Another important use
is for stabilizing parts that would otherwise have a rigid-body singularity. This is a
common problem in contact modeling before an assembly has actually settled. In this
case a Spring Foundation acting on the entire domain is useful because it avoids the
introduction of local forces.
A Thin Elastic Layer used as a pair condition can simulate thin layers with material
properties that differ significantly from the surrounding domains. Common
applications are gaskets and adhesives.
When a Thin Elastic Layer is applied on an interior boundary, it usually models a local
flexibility, such as a fracture zone in a geological model.
• Spring Data
• Loss Factor Damping
• Viscous Damping
SPRING DATA
The elastic properties can be defined either by a spring stiffness or by a force as function
of displacement. The force as a function of displacement can be more convenient for
nonlinear springs. Each spring node has three displacement variables defined, which
can be used to describe the deformation dependency. These variables are named
<interface_name>.uspring1_<tag>, <interface_name>.uspring2_<tag>, and
<interface_name>.uspring3_<tag> for the three directions given by the local
coordinate system. In the variable names, <tag> represents the tag of the feature
defining the variable. The tag could, for example, be spf1 or tel1 for a Spring
Foundation or a Thin Elastic Layer, respectively. These variables measure the relative
extension of the spring after subtraction of any predeformation.
f sl = 1 + i f s
VISCOUS DAMPING
It is also possible to add viscous damping to the Spring Foundation and Thin Elastic
Layer features. The viscous damping adds a force proportional to the velocity (or in
the case of Thin Elastic Layer: the relative velocity between the two boundaries). The
viscosity constant of the feature can be made dependent on the velocity by using the
variables named <interface_name>.vdamper1_<tag>, <interface_name>
.vdamper2_<tag>, and <interface_name>.vdamper3_<tag>, which contain the
velocities in the three local directions.
In this section:
• Thermal-Structure Interaction
• Acoustic-Structure Interaction
• Thermal-Electric-Structural Interaction
Thermal-Structure Interaction
The Thermal Stress, Solid Interface included with this module has a predefined
one-way approach for thermal-structure interaction (thermal stress), which combines
a Solid Mechanics interface with a Heat Transfer interface from the Heat Transfer
Module or COMSOL Multiphysics.
There are also similar multiphysics interfaces available for thin structures, as described
in The Thermal Stress, Shell Interface, The Thermal Stress, Membrane Interface, and
The Thermal Stress, Layered Shell Interface. The latter requires the Composite
Materials Module.
There are several physics interfaces available that are documented and
described in the Acoustic-Structure Interaction Interfaces chapter in the
Acoustic Module User’s Guide
Thermal-Electric-Structural Interaction
The Joule Heating and Thermal Expansion Interface enables
thermal-electric-structural interaction. This is a combination of three physics
interfaces: Solid Mechanics, Heat Transfer in Solids, and Electric Currents.
Using a single iteration does not produce a correct result if there are
thermal properties or electrical that depend on the displacements, making
the thermal-structure part into a bidirectional coupling.
Temperatures can either be computed using another physics interface, usually Heat
Transfer in Solids, or directly be prescribed in the input for the various physics nodes.
In this section:
If a material property under the Materials branch has a temperature dependence, you
have to input the temperature to be used in the Model Input section in the settings
window for the node in the physics interface that references the property. It is possible
that not all aspects of a material are defined in the same node in the Model Builder tree.
For example, if a problem is run with thermal expansion and plasticity, then:
• Young’s modulus, Poisson’s ratio, and mass density are given in the Linear Elastic
Material node.
For each of these nodes, there is a Model Input section in the Settings window. Some
of these sections may be empty if none of the properties given in that node have a
temperature dependence. In general, you have to supply the temperature in all the
Model Inputs sections.
As a default, the value of the temperature T is obtained from a Common model input.
You can also select User defined to enter a value or expression for the reference
temperature locally. This can be done either by explicitly giving a temperature or by
selecting a temperature variable from another physics interface.
If you want to create a model input value which is local to your current selection, click
the Create Model Input button . This will create a new Model Input node under
Definitions->Shared Properties in the current component, having the same selection as
in the current node.
See also Default Model Inputs and Model Input in the COMSOL
Multiphysics Reference Manual.
If any material in the model has a temperature-dependent mass density, the Volume
reference temperature list will appear in the Model Input section of the material settings.
As a default, the value of Tref is obtained from a Common model input. You can also
select User defined to enter a value or expression for the reference temperature locally.
All effects of volume change with temperature are incorporated through the thermal
expansion effects.
See also
Thermal Expansion
As the temperature changes, most materials react by a change of volume. For a
constrained structure, the stresses that evolve even with moderate temperature changes
can be considerable. The volume change can be represented a thermal strain th,
which produces stress-free deformations. For a linear elastic material, the constitutive
law is
= C – th
As long as you are using materials from the COMSOL Material Library, everything is
handled internally. When you want to enter data from your own measurements or from
the literature, you do, however, need to be aware of some details in the definitions
used.
dL
-------- = t T dT (2-29)
L
where t is the tangential thermal expansion coefficient. This form, which is the
thermodynamic definition, is conceptually simple, because t is uniquely defined at
each temperature. It is, however, less convenient to use in practice because an
integration is required for determining the actual change in length for a finite
temperature difference.
The secant formulation, which is the default in COMSOL Multiphysics, is often used
in engineering:
L
-------- = T T
L0
In the secant formulation, the actual values of will however depend on the choice of
reference temperature, Tref, at which the material has the reference length L0:
ln ------ =
L
L0 t d (2-30)
T ref
Define
Thus,
L L I T,T ref
-------- = ------ –1 = e –1
L0 L0
I T,T ref
e –1
T,T ref = ---------------------------------
T – T ref
For most materials and temperature ranges I T,T ref « 1 , which makes it possible to
approximate with the simpler expression
I T,T ref
T,T ref = ------------------------- (2-31)
T – T ref
If you have access to tangent data, you can choose between two different methods for
using them in COMSOL Multiphysics:
• In most of the physics interfaces, you can enter tangent data directly by selecting
Tangent coefficient of thermal expansion in the settings for Thermal Expansion. When
using this option, a numerical integration of Equation 2-30 will be performed each
L T = 1 + m T T – T m L T m (2-32)
When using the measured data, it is possible that the strain-free state occurs at a
temperature Tref which differs from Tm. The dilation at any temperature T would then
be defined as LTL(Tref, where L(Tref) can be written as.
Here r(T) is the redefined thermal expansion coefficient, based on Tref. It can be
derived from the relations above. Using Equation 2-32 and Equation 2-34 there are
two ways of writing the current length L(T), so that
m T – m T ref
m T + T ref – T m -----------------------------------------------
T – T ref
r T = ------------------------------------------------------------------------------------------------------ (2-36)
1 + m T ref T ref – T m
The Material Contents section in Figure 2-19 shows the material property alpha,
which is the redefined thermal expansion coefficient r(T). The complete expression
for alpha is as follows:
(alpha_solid_1(T[1/K])[1/K]+(Tempref-293[K])*
if(abs(T-Tempref)>1e-3,(alpha_solid_1(T[1/K])[1/K]
-alpha_solid_1(Tempref[1/K])[1/K])/(T-Tempref),
d(alpha_solid_1(T[1/K]),T)[1/K]))/
(1+alpha_solid_1(Tempref[1/K])[1/K]*(Tempref-293[K]))
This is essentially Equation 2-36, but with a small modification to avoid problems if
T=Tref.
r(T)
Figure 2-19: An example in COMSOL Multiphysics showing the Materials branch and
where to find the temperature-dependent thermal expansion coefficient.
Take care when describing the units. Temperature unit conversions can be
the cause of subtle errors because of the shift in zero-point value. Use
kelvin (K) as the temperature unit to the largest possible extent. As an
alternative, you can use the other absolute temperature scale, Rankine
(R). Avoid using Fahrenheit and Celsius unless you are completely
familiar with how the temperature unit conversion works.
In many cases, not only the structure which actually is modeled deforms due to the
changes in temperature, but also the surroundings (which are approximated by
constraints) will deform. You can take this effect into account by adding a Thermal
Expansion subnode to the constraints. The constraints will then provide an extra
displacement based on a given temperature field. For thermal strains, which have a
simple variation in space (for example, linear temperature variations), it is possible to
completely offset the constraint stresses using this method. For more general cases, the
stresses caused by the constraint can be significantly reduced.
The thermal expansions of the constraints are independent of that of the material in
the adjacent domain, so that the surrounding structure can be made from another
material, or have a different temperature distribution.
Thermoelastic Damping
In most engineering problems, the coupling between temperatures and structural
problems can be considered as unidirectional. Only the thermal expansion is
considered.
The opposite effect, where changes in stress cause heat generation, may be important
in small structures vibrating at high frequencies. The Thermoelasticity interface,
available with the MEMS Module, is designed for analyzing such problems.
It is also possible to take this effect into account by adding the Thermoelastic Damping
node to the Heat Transfer in Solids interface. When you add a Thermal Expansion node
to a material in the Solid Mechanics interface, the heat source term is computed and
made available to the Heat Transfer in Solids interface.
When you add a Thermal Expansion node under the Multiphysics Couplings branch, it is
possible to select whether the thermoelastic damping effect should be taken into
account or not. The heat source contribution is then included automatically without
adding any data in the heat transfer interface.
See also
• The pressure and viscous forces in the fluid provides a load on the boundary of the
solid. Usually, the pressure is dominant.
• The deformation of the structure changes the geometry of the fluid domain.
• The fluid sees the structure as a moving wall, which imposes a velocity at the
interface.
You can model FSI with four different structural mechanics interfaces: Solid
Mechanics, Multibody Dynamics, Shell, and Membrane. The fluid-flow can be
modeled with any domain-level physics interface from the Single-Phase Flow group and
the Two-Phase Flow groups under Multiphase Flow.
• Fluid-Solid Interaction
• Fluid-Shell Interaction
• Fluid-Membrane Interaction
• Fluid-Multibody Interaction
• Fluid-Multibody Interaction, Assembly
The Deforming Domain node is, however, not added for multiphysics interfaces denoted
‘Fixed Geometry’, which are intended for cases where the deformation of the fluid
domains is small everywhere.
A deforming domain represents domains and boundaries where the mesh can deform.
By default, the Deforming Domain node has an empty selection. You can then select any
fluid domain. However, this is only needed if the geometry of such a domain
experience significant changes due to the deformation or rotation of the adjacent solid
domains. Otherwise, the moving mesh computations could introduce unnecessary
Under the Deforming Domain, you can also choose the Mesh smoothing type, by default
set to Hyperelastic. More information of the smoothing type can be found in
Deforming Domain chapter in the COMSOL Multiphysics Reference Manual
By default, the mesh is free at all external boundaries of the geometry, and it follows
the solid boundaries at the solid-fluid interfaces. You can also add other types of
boundary conditions for the mesh motion; for details, see Deformed Geometry and
Moving Mesh in the COMSOL Multiphysics Reference Manual.
Union or Assembly
In most cases you model FSI problems so that the geometry sequence is set up to form
a union, and the same multiphysics coupling, Fluid-Structure Interaction, is used
irrespective of the type of structural mechanics interface. This coupling will
automatically find all boundaries that are shared between the structure and the fluid.
There are, however, cases where the assembly mode must be chosen, particularly when
having mechanisms, as is common in the Multibody Dynamics interface. In that case,
the interface between the solid and the fluid is no longer formed by a common
boundary. Rather, it consists of two boundaries, located at the same place in space.
These boundaries will in general slide with respect to each other. To model this, you
use the Fluid-Structure Interaction, Pair multiphysics coupling. You must create
appropriate pairs containing the boundaries from both types of physics under
Definitions, and manually select them in the Fluid-Structure Interaction, Pair node.
For the moving mesh, you must specify the deformation of the mesh manually. Add a
Prescribed Mesh Displacement node under the Moving Mesh node, in which you give the
structural displacement as the mesh displacement expression. The variables for the
displacement in the structure is provided by the multiphysics coupling.
By adding an offset in the settings for the Contact node, you can force the two sides of
the solid to experience contact at some distance before they meet in the geometrical
sense. This approach only will, however, leave a thin channel through which the fluid
can pass. The reduction in flow may be sufficient, but you can block it even further by
increasing the viscosity in the channel when the gap is closed. To do that, you can, for
example, compute the minimum gap anywhere in the contact pair, and then make the
viscosity a function of it. Another option is to compute the wall distance in the fluid
from both sides of the contact pair and use that information to modify the viscosity.
Do not increase the viscosity more than a couple of orders of magnitude, to avoid
numerical problems.
In configurations where you more or less completely cut off the whole flow, you must
pay particular attention to your boundary conditions. A prescribed flux will cause an
extreme pressure build-up upstream of the valve and thus unrealistically large forces on
the structure.
Deformed geometry FSI, or fixed geometry fully coupled FSI, should not
be solved using a one-way approach.
Below are the steps to follow to compute a one-way FSI problem sequentially:
1 In the study step settings windows, under the Physics and Variables Selection section,
clear the physics that is solved in the second step, so that only the governing physics
is selected.
2 Add a second study step to the study and in the settings windows make sure you
have the governing physics cleared.
If the selected study steps are of stationary type, you can generate the default solver
configuration, edit it if necessary, and compute the solution. The mapping of the
solution from the first to the second study step is done automatically.
In case of a transient problem, continue with the steps below:
3 In the second study step settings window, expand the section Physics and Variables
Selection. Under Initial values of variables solved for make sure the settings are
defined as in the table below:
In the case of a fluid loading to structure coupling type, the structural mechanics
problem can be treated as quasistatic. This can be handled by running the structural
analysis as a parametric sweep over a number of static load cases, where time is used as
the parameter.
The reduced component accounts for the constraints applied on the component. It can
also contain, for example, information about loads applied to the component prior to
model reduction.
Since the reduced components only have a small number of degrees of freedom (often
of the order of 10-100), they are computationally more efficient than the original full
FEM components. To create a reduced component, it is necessary to both perform an
eigenfrequency analysis and to solve for a number of static load cases. These studies
are, however, computed at the component level, and are thus usually computationally
much cheaper than analyzing the full model.
• Reduced-Order Modeling
• Model Reduction
COMSOL Implementation
The Component Mode Synthesis (CMS) technique is currently implemented in the
Solid Mechanics, Shell and Multibody Dynamics interfaces. You use the Reduced
Reduced components are connected to each other and to nonreduced parts of the
model through Attachment features. It means that all physics features that can access
an attachment can be used to connect the reduced components. The most common
option is to use joints, but there are other alternatives, such as Spring-Damper, or
bearing foundations in the rotordynamics interfaces.
The global model can consist of either a set of reduced components, or a mix of
ordinary FE discretized domains and reduced components. The reduced components
can be used as linear parts in an otherwise nonlinear model.
You can compute all type of results and visualize reduced components just as if they
are ordinary domains in your model. One important property of the reduced
component is that only the results for the reduced set of degrees of freedom are stored.
It means that for time-dependent studies with many time steps, the file sizes can be
reduced by orders of magnitude.
Applications
Some examples of applications of the CMS reduction technique include:
By default, disconnected geometries connect by the Continuity and Thin Elastic Layer
pair features are merged when automatically generating Component Definition
subnodes. This also applies to Boundary to Boundary, Edge to Boundary, and Edge to
Edge connection features in the Shell interface. The automatic handling of such
connections can be disabled by clearing Include connections and pairs in component
definition in the Reduced Flexible Components node.
When working with CMS, the full static and dynamic behavior of a component is
represented by a number constraint modes and constrained eigenmodes that are
computed in the training phase. These are independent of each other and are used to
construct one ROM for each component. Each Attachment node connected to a
component will add a number of static load cases that describe the constraint modes:
six in 3D, and three in 2D. The number of eigenmodes to be used is controlled
manually. You can either define it for all components in the Reduced Flexible
Components node, or individually in each Component Definition subnode. To get an
accurate representation of the reduced component, always make sure to use a sufficient
number of eigenmodes to describe the dynamics of each component. During the
eigenfrequency training step, all attachments are treated as fixed constraints, hence, the
eigenmodes are always constrained.
To create the ROMs, a special study sequence needs to be set up for each Reduced
Flexible Components node. It should sweep over all Component Definition subnodes, and
for each component, compute the static load cases and the requested eigenmodes in
training study steps as outlined above. The results from these study steps are then used
in a Model reduction step to generate a ROM. By using the Configure CMS Study ( )
button in the Reduced Flexible Components node, the set-up of this special study
Figure 2-21: Nodes in the model tree that are added automatically when working with
CMS.