Shrink-fitting with code_aster: two numerical approaches for the same stress state

Introduction

A shrink-fit assembly is usually obtained by heating the hub, inserting the shaft into the expanded bore, then letting the whole thing cool down. Simulating this entire sequence is possible, but it isn't always necessary when the goal is simply to find the residual stress state. Here we compare two more direct ways of introducing the interference in a code_aster model: an interface that is initially flush with a contact offset, and a geometry that is initially interpenetrated. How far can these two shortcuts be pushed before they stop agreeing with each other, and what do they actually represent of the real industrial process?

In short

Two axisymmetric models of the same shaft-hub shrink-fit were built with code_aster. They share the same final dimensions, the same materials, the same frictional contact, and the same elastoplastic law. Only the initial geometric configuration and the way the interference is introduced differ.

• In the flush approach, the shaft and hub surfaces initially coincide. The interference is imposed virtually through the contact definition.
• In the interpenetrated approach, the interference is carried by the initial geometry: the bore radius is smaller than the shaft radius.

Both models produce very close contact pressures: 52.25 and 52.55 MPa in the elastic regime, then 256.39 and 257.39 MPa once plasticity sets in. The gap on the average pressure stays at 0.57% and 0.39% respectively.

This agreement shows that, for this case, both numerical representations of the interference lead to essentially the same final mechanical state.

Real shrink-fitting is first and foremost a thermal process

The principle of shrink-fitting is to assemble two parts whose dimensions at room temperature are incompatible. The shaft radius is slightly larger than the bore radius. A direct assembly would therefore require a large force and risk damaging the surfaces.

In a hot shrink-fitting process, the sequence usually goes as follows:

1. 1/ the hub is heated;
2. 2/ its thermal expansion increases the bore diameter;
3. 3/ the shaft is inserted while the surfaces are separated or lightly loaded;
4. 4/ the hub cools down and contracts;
5. 5/ the clearance closes, contact is established, and the pressure rises;
6. 6/ the assembly reaches its residual state at service temperature.

The most faithful simulation would therefore reproduce this whole thermal and mechanical history: initial hub temperature, heat exchange with the shaft and the environment, expansion, progressive onset of contact, cooling, and possible plastification.

This strategy, however, requires additional data: expansion coefficients, temperature-dependent properties, exchange conditions, assembly temperature and duration, and possibly a friction coefficient that evolves with the surface state. It also increases the cost of preparing and solving the computation.

When the goal is only to know the pressure and stresses after returning to room temperature, a natural question arises: can the thermal process be replaced by a purely mechanical interference?

Replacing the process by its final effect

Both approaches studied here rely on the same simplification: neither simulates the heating nor the cooling. They directly impose the geometric incompatibility that remains at the final temperature.

This reduction is particularly attractive in an industrial context:

• the computation doesn't need a thermal model;
• the only assembly data required are the final dimensions and the interference;
• the contact pressure and residual stresses are obtained in a single mechanical computation;
• the model is faster to build, to converge, and to use in a parametric study.

The difficulty is shifted onto the initialization of the contact. In the undeformed state, two solids can't occupy the same space, while the final tightening precisely imposes a dimensional incompatibility. Solvers therefore offer several tricks to introduce this interference.

We compare two of them here.

Approach 1 — Flush mesh with the interference carried by contact

In the first model, the shaft and the hub are built with the same radius at the interface. Their surfaces are geometrically flush, but their nodes remain distinct so that contact can be created.

The tightening is introduced using a slave distance imposed in the contact definition. This distance virtually thickens the hub surface towards the shaft. The solver then detects a virtual penetration and generates the displacements and pressure needed to remove it.

The interference is applied progressively:

• it is zero at the start of the computation;
• it reaches 0.015 mm at a first state meant to remain elastic;
• it then increases up to 0.100 mm, a level at which plastification appears.

In the code_aster input, this evolution is carried by a positive function assigned to DIST_ESCL. The following excerpt is deliberately limited to the essential contact keywords:

# Both surfaces are geometrically flush.
# DIST_ESCL progressively creates the virtual interference.
DIST_ES = DEFI_FONCTION(
    NOM_PARA="INST",
    VALE=(
        0.0,  0.000,
        0.15, 0.015,
        1.0,  0.100,
    ),
)

CTAC = DEFI_CONTACT(
    MODELE=MODE,
    FORMULATION="CONTINUE",
    FROTTEMENT="COULOMB",
    ZONE=_F(
        GROUP_MA_MAIT="CONT_ARBRE",
        GROUP_MA_ESCL="CONT_MOYEU",
        ALGO_CONT="STANDARD",
        ALGO_FROT="PENALISATION",
        CONTACT_INIT="NON",
        DIST_ESCL=DIST_ES,
        COULOMB=0.15,
    ),
)

The positive sign of DIST_ESCL here virtually thickens the slave surface towards the shaft. The final value of the function is therefore directly equal to the target radial interference. CONTACT_INIT="NON" is consistent with there being no geometric overlap at the start of the computation.

This approach has a real practical advantage: the mesh can be built from a common nominal geometry, with no visible initial overlap. The interference becomes a parameter of the computation and can be changed without regenerating the mesh.

In return, it requires a good grasp of sign conventions, normal orientation, and the role of the master and slave surfaces. A mistake on any of these choices can open the contact instead of closing it, or double the intended interference.

Approach 2 — Interference carried by the geometry

In the second model, the bore is directly built with a radius smaller than the shaft's. The two domains therefore overlap, in the initial configuration, by exactly the amount of the final tightening.

A contact distance opposite to this overlap is applied at the start of the computation to create an initial state with no mechanical load. It is then progressively brought back to zero. In the final state, the compensation has disappeared and the full geometric interference acts on the contact.

The contact definition this time uses a negative initial value. For a geometric overlap of 0.100 mm, the ramp is written as:

# The mesh already contains a radial overlap of 0.100 mm.
# DIST_ESCL compensates it at first, then releases the tightening.
DIST_ES = DEFI_FONCTION(
    NOM_PARA="INST",
    VALE=(
        0.0,  -0.100,  # effective interference: 0.000 mm
        0.15, -0.085,  # effective interference: 0.015 mm
        1.0,   0.000,  # effective interference: 0.100 mm
    ),
)

CTAC = DEFI_CONTACT(
    MODELE=MODE,
    FORMULATION="CONTINUE",
    FROTTEMENT="COULOMB",
    ZONE=_F(
        GROUP_MA_MAIT="CONT_ARBRE",
        GROUP_MA_ESCL="CONT_MOYEU",
        ALGO_CONT="STANDARD",
        ALGO_FROT="PENALISATION",
        CONTACT_INIT="INTERPENETRE",
        DIST_ESCL=DIST_ES,
        COULOMB=0.15,
    ),
)

The interference actually handled by the solver is here 0.100 mm + DIST_ESCL. The keyword CONTACT_INIT="INTERPENETRE" tells code_aster to take charge of the initially overlapping configuration. Once DIST_ESCL reaches zero, no virtual compensation remains: contact resolves the full interference carried by the mesh.

This representation is intuitive when the manufactured dimensions of the shaft and bore are known: the interference is visible in the geometry and directly corresponds to the difference in radii.

It comes, however, with several constraints:

• the tightening value must stay consistent between the geometry and the contact definition;
• changing the interference may require regenerating the mesh;
• the initial overlap doesn't represent an admissible physical configuration, but a numerical state to be corrected;
• a poor initialization can immediately generate very large forces and jeopardize convergence.

A common protocol to isolate the effect of the method

For the comparison to be meaningful, both models use the same mechanical framework:

• 2D axisymmetric geometry;
• solid shaft of radius 25 mm;
• hub of outer radius 60 mm and height 100 mm;
• small strains;
• Coulomb frictional contact, with a coefficient of 0.15;
• von Mises plasticity with linear isotropic hardening;
• Young's modulus of 210,000 MPa;
• Poisson's ratio of 0.30;
• yield stress of 300 MPa;
• hardening modulus of 2,000 MPa.

The shaft is defined as the master surface and the hub as the slave surface in both computations. The same effective interference history is applied. Results are compared on the pressure along the interface and on the stresses extracted along a radial line at $z=20$ mm, away from the fixed base.

This protocol doesn't yet aim to reproduce a particular industrial assembly. It aims to answer a more basic question: for identical physics and loading, do the two numerical encodings of the tightening lead to the same stress state?

In the elastic regime, the two approaches overlap

At the first state studied, the effective radial interference is 15 µm. The behavior remains elastic throughout the shaft and hub.

15 µm interference	Flush model	Interpenetrated model	Gap
Average contact pressure	52.25 MPa	52.55 MPa	0.57%
Maximum local pressure gap	—	—	0.30 MPa

The pressure is nearly uniform over most of the interface. The largest variations appear near the ends, where boundary conditions and the termination of the contact surfaces break the idealization of an infinitely long assembly.

The radial and hoop stress profiles along this line are also very close. In the shaft, the radial and hoop stresses are practically uniform and compressive. In the hub, the radial compression decreases towards the outer surface while the hoop stress is maximal near the bore.

At this loading level, the choice between a flush geometry and an interpenetrated geometry therefore has very little influence on the resulting mechanical state.

After plastification, the agreement remains very good

The interference is then raised to 100 µm. The yield stress is exceeded in the hub near the interface and a plastic zone develops.

100 µm interference	Flush model	Interpenetrated model	Gap
Average contact pressure	256.39 MPa	257.39 MPa	0.39%
Maximum local pressure gap	—	—	1.05 MPa

Despite the nonlinearity of both the material and the contact, the average pressures differ by less than 1%. The stress fields show the same overall organization, and plasticity develops in the same region.

The V1 component of the internal variables of the VMIS_ISOT_LINE behavior represents the cumulated equivalent plastic strain. It is zero at 15 µm. At 100 µm, its maximum value is around $2.7\times10^{-3}$ and the non-zero zone is concentrated near the interface in the hub.

This second comparison is more demanding than the first: in an elastoplastic problem, the result can depend on the loading path. The very good agreement obtained indicates that the two interference ramps produce mechanical histories close enough here to lead to almost the same final state.

What the comparison shows (and what it doesn't)

The main result is the low sensitivity of the computation to the way the interference is introduced. For this model:

• the average pressure varies by less than 0.6% between the two approaches;
• the conclusion holds both before and after the onset of plasticity;
• both methods can be used as mechanical representations of the final tightening.

This equivalence is valuable for the engineer. It makes it possible to pick whichever strategy best fits the modeling workflow: a reusable nominal mesh with the flush approach, or a geometry that explicitly carries the manufactured dimensions with the interpenetrated approach.

But both models share the same shortcut: neither simulates the heating or the cooling. Their agreement is a verification of the interference implementation, not a validation of the actual shrink-fitting process.

This distinction matters as soon as the assembly history influences the final state.

When is the mechanical shortcut representative?

With purely elastic behavior, uniform temperatures, and a contact with no complex evolution, the final state mainly depends on the interference reached. Replacing the cooling by a mechanical interference is then generally an efficient approximation.

Representativeness becomes more uncertain when path-dependent phenomena come into play:

• plasticity during cooling: the yield stress and hardening modulus may depend on temperature;
• thermal gradient in the hub: the outer surface and the bore don't necessarily cool at the same rate;
• heating of the shaft after insertion: heat transfer temporarily reduces the clearance and shifts the moment contact is established;
• sliding with friction: the progressive closing of the contact can produce a sliding history different from that of a mechanically imposed interference;
• initial residual stresses: manufacturing, machining, or heat treatment can alter the state before assembly;
• temperature-dependent properties: expansion, Young's modulus, yield stress, and friction aren't necessarily constant.

In these situations, two mechanical models can agree perfectly with each other while still departing from the real process. The right question is no longer just "do the two methods give the same result?", but "can the phenomena left out change the design decision?".

Lamé as a complementary check

The Lamé analytical solution remains useful, but it plays a secondary role here. It provides an independent check on the pressure level and on the shape of the stress profiles as long as the behavior stays elastic.

The Lamé solution comes from the radial equilibrium of an axisymmetric cylinder in linear elasticity:

$$\frac{d\sigma_{rr}}{dr}+\frac{\sigma_{rr}-\sigma_{\theta\theta}}{r}=0,$$

which leads to the profiles $\sigma_{rr}=A-B/r^2$ and $\sigma_{\theta\theta}=A+B/r^2$. For a solid shaft under external pressure $p$, this gives $\sigma_{rr}=\sigma_{\theta\theta}=-p$. In the hub, with $\sigma_{rr}(R_a)=-p$ and $\sigma_{rr}(R_e)=0$:

$$\sigma_{rr}=p\frac{R_a^2}{R_e^2-R_a^2}\left(1-\frac{R_e^2}{r^2}\right), \qquad \sigma_{\theta\theta}=p\frac{R_a^2}{R_e^2-R_a^2}\left(1+\frac{R_e^2}{r^2}\right).$$

The pressure $p$ is then found by imposing that the expansion of the hub and the contraction of the shaft together absorb the radial interference $\delta$. These expressions assume purely elastic behavior, axisymmetric geometry, a pressure uniform along the axis, and the absence of end effects.

For a 15 µm interference, the Lamé variant consistent with a plane-stress-like state predicts a pressure of 52.06 MPa. The flush model gives 52.25 MPa, a 0.36% gap. The axial stress extracted at $z=20$ mm is small compared to the radial and hoop stresses, which confirms that this plane-stress assumption is appropriate for the case studied.

At 100 µm, the fully elastic Lamé solution is no longer a valid quantitative reference overall: it predicts a pressure of 347.08 MPa, while code_aster gives 256.39 MPa. The plasticity developed in the hub therefore reduces the contact pressure and also affects the stress state of the shaft, even though the shaft itself stays elastic.

The internal variable V1 makes it possible to locate this loss of validity. It becomes non-zero between the interface and a plastic front located around $r=39.9$ mm. In this zone, the numerical profiles clearly depart from the elastic solution. Beyond the plastic front, in the outer part of the hub that remains elastic, the radial and hoop stresses progressively rejoin the Lamé profiles.

This reading should stay nuanced: plasticity is local, but its effect on the contact pressure is global. The absence of plasticity at a given point doesn't by itself guarantee agreement with the fully elastic solution. This is especially visible in the shaft, where the stresses remain uniform but follow the contact pressure reduced by plastification in the hub.

This analytical comparison is therefore an additional check, not the central subject of the benchmark. It verifies that the elastic pressure level is consistent and helps pinpoint the moment an elastic formula stops being appropriate.

Conclusion

This shrink-fitting case shows that two very different numerical representations of the interference can lead to the same stress state.

With a nominally flush interface, the tightening is carried by the contact definition. With an interpenetrated geometry, it is carried by the initial dimensions and then progressively released. The gaps on the average pressure stay below 0.6%, including after a plastic zone appears.

For a study focused on the residual state at room temperature, both methods are therefore consistent and efficient mechanical shortcuts. The choice can be guided by ease of parametrization, mesh management, and the robustness of the contact initialization.

Their agreement, however, doesn't replace an analysis of how representative the model is of the real process. A shrink-fit is obtained through a thermal history, and that history can become decisive once plasticity, temperature gradients, or friction influence the final state.

---

Computation carried out with code_aster. Pre- and post-processing can be done in SALOME and ParaVis. The values shown correspond to a demonstration case and don't constitute generic design data.