Penalty Method-Based Contact Formulation
Theory and Physics
Fundamentals of Contact Problems
Professor, what makes FEM "contact problems" so difficult?
Contact problems are inherently nonlinear. There are three reasons:
1. State Change — The contact surface has a binary state of "Contact/No Contact". The state changes with load.
2. Inequality Constraint — Penetration must not occur: $g \geq 0$
3. Friction — Coulomb friction involves state changes between "stick/slip".
Inequality constraints... That's fundamentally different from the equality constraints (like SPC) in regular FEM.
Regular FEM solves the equality $[K]\{u\} = \{F\}$, but contact must satisfy the KKT conditions (Karush-Kuhn-Tucker conditions): "$g \geq 0$ and $g \cdot p = 0$". $g$ is the gap, $p$ is the contact pressure.
Principle of the Penalty Method
The Penalty Method is the simplest contact handling technique. It adds a "reaction force proportional to penetration":
$k_p$ is the penalty stiffness (contact stiffness). $g_n$ is the normal penetration amount.
It's like putting a spring at the contact surface.
Exactly. The penalty method places a very stiff spring at the contact surface. When penetration occurs, a reaction force is generated, pushing back the penetration.
Setting Penalty Stiffness
The setting of penalty stiffness $k_p$ greatly influences the results:
- $k_p$ too large → Condition number worsens. Convergence becomes difficult.
- $k_p$ too small → Penetration amount becomes excessive. Inaccurate.
Guideline: $k_p \approx 10 \sim 100 \times E \cdot A / L$ (on the order of the contact surface stiffness). Many solvers calculate this automatically.
Is it safe to rely on automatic calculation?
For most cases, yes. The automatic penalty stiffness in Abaqus or Ansys is a good starting point. Only adjust manually when problems arise.
Advantages and Disadvantages of the Penalty Method
| Advantages | Disadvantages |
|---|---|
| Simple implementation | Penetration does not become exactly zero |
| No additional DOFs required | Penalty stiffness setting affects results |
| Compatible with explicit methods | Large $k_p$ worsens condition number |
| Default for many solvers | — |
Is "penetration not becoming zero" the biggest weakness?
As long as $k_p$ is finite, slight penetration is unavoidable. If the penetration amount is less than 1% of the plate thickness, it's practically acceptable. If more, consider the Lagrange multiplier method.
Summary
Key Points:
- Contact is a nonlinear problem with inequality constraints — State change + Friction
- Penalty Method = Spring at contact surface — $F = k_p \cdot g_n$
- $k_p$ setting is key — Too large causes convergence difficulty, too small causes excessive penetration
- Many solvers calculate automatically — Manual adjustment only when problems occur
- Penetration does not become exactly zero — Aim for less than 1% of plate thickness
Courant's 1943 Penalty Method
The mathematical origin of the penalty method lies in Richard Courant's 1943 paper "Variational methods for the solution of problems of equilibrium and vibrations". The idea of transforming a constrained variational problem into an unconstrained one by adding a term multiplied by a large coefficient (penalty) was initially used for handling boundary conditions of elliptic partial differential equations. In the 1970s, Oden & Kim (1977) and others systematically applied it to contact FEM.
Physical Meaning of Each Term
- Inertia Term (Mass Term): $\rho \ddot{u}$, i.e., "mass × acceleration". Have you ever experienced being thrown forward when slamming on the brakes? That "feeling of being carried away" is precisely the inertial force. Heavier objects are harder to start moving and harder to stop once moving. Buildings shake during earthquakes because the ground moves suddenly while the building's mass "gets left behind". In static analysis, this term is set to zero, assuming "forces are applied slowly so acceleration can be ignored". It cannot be omitted in impact loads or vibration problems.
- Stiffness Term (Elastic Restoring Force): $Ku$ or $\nabla \cdot \sigma$. When you pull a spring, you feel a "force trying to return it", right? That's Hooke's law $F=kx$, the essence of the stiffness term. Now a question — an iron rod and a rubber band, which stretches more when pulled with the same force? Obviously the rubber. This "resistance to stretching" is the Young's modulus $E$, which determines stiffness. A common misconception: "High stiffness = strong" is incorrect. Stiffness is "resistance to deformation", strength is "resistance to failure" — different concepts.
- External Force Term (Load Term): Body force $f_b$ (gravity, etc.) and surface force $f_s$ (pressure, contact force, etc.). Think of it this way — the weight of a truck on a bridge is a "force acting on the entire volume" (body force), the force of the tires pushing on the road is a "force acting only on the surface" (surface force). Wind pressure, water pressure, bolt tightening force... all are external forces. A common mistake here: getting the load direction wrong. Intending "tension" but it becomes "compression" — sounds like a joke, but it actually happens when coordinate systems are rotated in 3D space.
- Damping Term: Rayleigh damping $C\dot{u} = (\alpha M + \beta K)\dot{u}$. Try plucking a guitar string. Does the sound continue forever? No, it gradually fades. That's because vibration energy is converted to heat by air resistance and internal friction in the string. Car shock absorbers work on the same principle — intentionally absorbing vibration energy for a smoother ride. What if damping were zero? Buildings would keep shaking forever after an earthquake. Since that doesn't happen in reality, setting appropriate damping is important.
Assumptions and Applicability Limits
- Continuum assumption: Treats material as a continuous medium, ignoring microscopic heterogeneity.
- Small deformation assumption (for linear analysis): Deformation is sufficiently small compared to initial dimensions, stress-strain relationship is linear.
- Isotropic material (unless specified otherwise): Material properties are independent of direction (anisotropic materials require separate tensor definition).
- Quasi-static assumption (for static analysis): Ignores inertial and damping forces, considers only balance between external and internal forces.
- Non-applicable cases: Large deformation/large rotation problems require geometric nonlinearity. Nonlinear material behavior like plasticity or creep requires constitutive law extension.
Dimensional Analysis and Unit Systems
| Variable | SI Unit | Notes / Conversion Memo |
|---|---|---|
| Displacement $u$ | m (meter) | When inputting in mm, unify loads and elastic modulus to MPa/N system. |
| Stress $\sigma$ | Pa (Pascal) = N/m² | MPa = 10⁶ Pa. Be careful of unit inconsistency when comparing with yield stress. |
| Strain $\varepsilon$ | Dimensionless (m/m) | Note the distinction between engineering strain and logarithmic strain (for large deformation). |
| Elastic Modulus $E$ | Pa | Steel: ~210 GPa, Aluminum: ~70 GPa. Note temperature dependence. |
| Density $\rho$ | kg/m³ | In mm system: tonne/mm³ (= 10⁻⁹ tonne/mm³ for steel). |
| Force $F$ | N (Newton) | Unify as N in mm system, N in m system. |
Numerical Methods and Implementation
Implementation of the Penalty Method
How is the penalty method set up in each solver?
Abaqus
```
*CONTACT PAIR, INTERACTION=contact_prop
slave_surface, master_surface
*SURFACE INTERACTION, NAME=contact_prop
*SURFACE BEHAVIOR, PENALTY
```
Abaqus default is the penalty method. Specify explicitly with *SURFACE BEHAVIOR, PENALTY.
Nastran
```
BCTABLE, ...
BCPARA, CTYPE, UGLYPEN $ Penalty Method
```
Contact definition in Nastran SOL 101/106/400.
Ansys
```
MP, MU, cid, 0.3 ! Friction coefficient
KEYOPT, cid, 2, 1 ! Penalty Method
```
Or in Workbench Contact settings, Formulation=Penalty.
LS-DYNA
```
*CONTACT_AUTOMATIC_SURFACE_TO_SURFACE
$ Penalty method is default
```
All LS-DYNA contacts are penalty method based. The most natural method for explicit analysis.
Contact Detection
Penalty method procedure (each iteration):
1. Contact Detection — Search for the closest point of slave nodes on the master surface.
2. Gap Calculation — $g_n$ = distance between slave node and master surface.
3. Penetration Judgment — If $g_n < 0$, contact. Reaction force $F_n = k_p |g_n|$.
4. Friction Force Calculation — $F_t = \mu F_n$ (slip) or stick.
5. Reflect in Global Equation — Add penalty force to the right-hand side.
Does contact detection account for most of the computational cost?
For large-scale models (hundreds of thousands of contact element pairs), contact detection can account for 30-50% of computation time. Fast search algorithms like Bucket sort and KD-tree are used.
Summary
Selection Rule for Penalty Stiffness
Choosing the contact penalty coefficient ε is the most important practical issue in the penalty method. If ε is too small, penetration exceeds allowable values; if too large, the condition number of the stiffness matrix worsens and convergence slows. LS-DYNA's default setting uses the "SOFT=0" algorithm which automatically sets ε to 0.1 times the minimum element stiffness k_e (=E×t×A) of the contact surface. For dissimilar material contact (e.g., steel-resin), manual tuning is recommended.
Linear Elements (1st Order Elements)
Linear interpolation between nodes. Computational cost is low, but stress accuracy is low. Beware of shear locking (mitigated by reduced integration or B-bar method).
Quadratic Elements (with Mid-side Nodes)
Can represent curved deformation. Stress accuracy improves significantly, but degrees of freedom increase about 2-3 times. Recommended: when stress evaluation is important.
Full Integration vs Reduced Integration
Full Integration: Risk of over-constraint (Locking). Reduced Integration: Risk of hourglass mode (zero-energy mode). Choose appropriately.
Adaptive Mesh
Automatic refinement based on error indicators (e.g., ZZ estimator). Efficiently improves accuracy in stress concentration areas. There are h-method (element subdivision) and p-method (order increase).
Newton-Raphson Method
Standard method for nonlinear analysis. Updates tangent stiffness matrix each iteration. Achieves quadratic convergence within convergence radius, but computational cost is high.
Modified Newton-Raphson Method
Updates tangent stiffness matrix using initial value or every few iterations. Cost per iteration is low, but convergence speed is linear.
Convergence Criteria
Force residual norm: $||R|| / ||F_{ext}|| < \epsilon$ (generally $\epsilon = 10^{-3}$〜$10^{-6}$). Displacement increment norm: $||\Delta u|| / ||u|| < \epsilon$. Energy norm: $\Delta u \cdot R < \epsilon$
Load Increment Method
Instead of applying full load at once, apply in small increments. The Arc-length method (Riks method) can track beyond extremum points of the load-displacement curve.
Analogy: Direct Method vs Iterative Method
The direct method is like "solving simultaneous equations accurately with pen and paper" — reliable but takes too long for large problems. The iterative method is like "repeatedly guessing to approach the correct answer" — starts with a rough answer but accuracy improves with each iteration. It's the same principle as looking up a word in a dictionary: opening to an estimated page and adjusting forward/backward (iterative method) is more efficient than searching sequentially from the first page (direct method).
Related Topics
なった
詳しく
報告