I am an orthopaedic surgeon working on a finite element comparison of fibula fracture fixation devices, and I would be grateful for advice on a contact convergence problem. I am running CalculiX 2.21 (built with PARDISO, Linux, 16 threads). The decks are generated by my own scripts (gmsh mesh) and solved directly with ccx rather than built in the PrePoMax GUI, but I hope the contact question is relevant here as many of you know CalculiX contact very well.
Model
CT-based fibula (distal 150 mm), C3D10 mesh (~285 k nodes), linear elastic materials (cortical / cancellous / marrow, E 7.3 GPa / 1.1 GPa / 2 MPa; implant Ti 113.8 GPa).
A 3.6 mm titanium intramedullary nail inside the medullary canal. The canal was cut by an exact boolean of the nail, so the canal surface and the nail surface are coincident conforming meshes (nail-side nodes duplicated from the bone nodes; curved interface ≈ 1600 mm²).
Bone elements along the canal are mixed: by area roughly 60 % cancellous, 31 % marrow (E 2 MPa) and 9 % cortical.
The nail is bonded to the bone at both ends (distal locking screws and a proximal locking element share nodes with the bone), so there are no rigid body modes.
Two contact pairs, both face-to-face penalty (TYPE=SURFACE TO SURFACE):
nail–canal: same normal stiffness, *FRICTION 0.3, 200
Load: *STEP, INC=1000, *STATIC 0.125, 1.0, 1e-6, 0.25; 70 N axial force at a reference node with a *DISTRIBUTING coupling on the distal end; proximal 5 mm fixed. (Step 2 adds 1 N·m torque.)
What happens
The same bone/fracture model with an intramedullary screw instead of the nail (fracture pair only) converges in 15–50 iterations per increment.
With the nail–canal pair, the first increment never converges: 60 iterations, cutback, repeated three times. The residual force drops to single digits and then diverges again (e.g. 3.4 → 1800). The number of contact elements in the .cvg file oscillates strongly between iterations (≈ 55 k ↔ 186 k).
Without the proximal locking element the residual stayed around 1e5, so the rigid-body issue is solved; what remains looks like contact-state chatter on the coincident curved surface.
While auditing the mesh I found one quadratic mid-edge node on the interface that was still shared between nail and canal although all corner nodes of that face were duplicated. I have fixed this (mid-edge nodes are now duplicated whenever all three corners are), and I am currently running a diagnostic series on a coarser mesh: nail–canal contact off, frictionless variants, and a frictionless normal-stiffness sweep (2e3 / 2e4 / 2e5 N/mm³).
Questions
Is open/close chattering expected for face-to-face penalty contact when the two surfaces are exactly coincident and curved? Is there a recommended way to stabilise it (e.g. *CONTROLS, PARAMETERS=CONTACT — delcon / alea / kscalemax / itf2f)? I could not find much guidance on how these are used in practice.
Could the very soft marrow elements (2 MPa) behind part of the canal surface be the main driver, given a penalty of 2e4 N/mm³? Is there a rule of thumb for choosing the face-to-face penalty relative to the underlying element stiffness in such a mixed-stiffness case?
Would node-to-face penalty contact (with SMALL SLIDING) be the better formulation for this conforming case?
Would you rather recommend a small initial radial clearance between nail and canal (e.g. 0.05–0.1 mm, as a reamed canal would have) instead of exactly coincident surfaces?
Any experience with the mortar options (TYPE=MORTAR / LINMORTAR) for this kind of problem with friction?
I can share .cvg / .sta excerpts and the full contact definitions. The geometry comes from patient CT data, so I would prefer not to post the complete mesh, but I can build a simplified synthetic test case if that helps.
Yes, chattering may happen. Unfortunately, CalculiX (unlike Abaqus) doesn’t provide contact stabilization/damping to reduce it. Usually the best recommendation is to reduce the contact stiffness. Also, whenerver possible, use first-order hex elements in contact analyses. And make sure your mesh is dense enough.
Small sliding may indeed improve convergence provided that there’s no significant tangential motion between the contact surfaces. Node to face contact is usually worse for pressure distribution, more prone to master-slave assignment and second-order element issues, though.
You can add a small gap or even interference and close them with ADJUST in contact to make sure the slave surface lies exactly on the master surface. Otherwise, there might be initial rigid body motions.
Mortar contact can provide more accurate contact pressure results, but it’s unlikely to converge better. In fact, it often diverges, especially in CalculiX 2.22 (you should switch to 2.23 if you want to try this type of contact).
On the CalculiX forum you can find some further discussions about contact types and constants (also for sliding). Perhaps you should ask there, but if you need more detailed help, it would be best if you could share the .pmx file (can be via private messages).
Thank you very much @FEAnalyst, this is really helpful, and it is reassuring to know that help is available here.
A short status update, in case it is useful to others:
I ran a diagnostic series on a coarser mesh. With both interfaces frictionless, the nail model now converges without a single cutback (13 increments, ~10–25 iterations each). With the nail–canal pair removed it also converges. The difficulty starts when the frictional fracture contact (μ 0.46) is combined with the nail–canal contact.
Using the K = 0.5·(E1+E2)/L_char rule from the thread you linked, our fracture-plane penalty (2e4 N/mm³) is roughly 2× too high for cortical and >10× too high for cancellous bone, while the nail–canal value is, if anything, on the low side. So a lower fracture-plane K is the next thing we will try.
Following your suggestion in Q4, a run with a small uniform nail–canal clearance (*CLEARANCE 0.02 / 0.05 mm) is running now.
I don’t have a .pmx (the decks are generated by scripts), but the CT data come from an open dataset, so if we still cannot get the frictional case to converge, I would be glad to share the input deck with you. Thank you again for the offer, it gives us a lot of confidence.
Friction can help eliminate RBMs (provided that contact is established properly), but it’s also really bad for convergence. Especially with higher coefficients. If possible, it’s good to use displacement control instead of force control. You can also try adjusting the constants for the frictional behavior, but I’d rather simplify/reduce the mode for debugging first (avoid other nonlinearities for now) and make sure that a good mesh is used (possibly hexa).
Yeah, even if you use characteristic length based on the average element edge length in the contact region, the stiffness may still be too high (the default one is usually way too large). So, if it still diverges, you can try reducing it further, but watch out for excessive penetrations.
You can also try adjusting *CONTROLS if everything else fails. It has to be done carefully, though.
*CLEARANCE can be a good way to avoid discretization errors and model interference fit or even bolt preload, but in CalculiX it used to have some issues. With ccx 2.22 or 2.23 it should be fine, though. Especially with the latter release if you want to use Mortar formulation too.
So you work with pure .inp files ? Ok, I see. Then it would be even more recommended to post here where more CalculiX users (and also developers) are present and can have a look at the input decks: https://calculix.discourse.group/
Thank you again, this is very useful, especially the point about displacement control and the warning about *CLEARANCE in older releases (we are on 2.21, so I will keep that in mind when interpreting the clearance runs).
For now we are testing a lower fracture-plane K, a lower nail–canal friction (μ 0.15, which is also used in the literature) and, as a fallback, a tied nail–canal interface (consistent with how the screw–bone interfaces of the other implants are modelled). I will post the outcome here, and once my CalculiX forum account is approved I will share the decks there as you suggested.
Maybe you could also consider using FEBio. This open-source software is meant for biomechanics and offers more material models for that. However, CalculiX should be fine too as long as regular elasticity (nonlinear needed for larger strains) and plasticity is sufficient.
Contact interfaces for such implants are indeed often modeled as bonded, especially if there’s some degree of osseointegration. Tie constraints are the most robust, but tied contact allows you to check the contact pressure at the interface. Without smoothing or at least hex elements, it can be quite irregular, though.
Thank you again for your detailed suggestions; they shaped how we debugged this. A short summary in case it helps someone with a similar problem:
Solver switch. We moved the production runs from CalculiX 2.21 to FEBio 4.13 (sliding-elastic contact, full Newton, non-symmetric stiffness). The deck is converted from the same Abaqus-style mesh, and loads/BCs are identical. On the intact bone model, CalculiX (NLGEOM) and FEBio agree to 0.00 %.
Two bugs in our own converter. (a) The surface face winding was inverted, so contact normals pointed inward, and FEBio treated penetration as separation (huge first-iteration displacement, negative Jacobians). (b) At exactly zero gap, sliding-elastic contributes no stiffness, so we use a tiny offset.
The real culprit: marrow in the contact surfaces. Our fracture and canal surfaces included the faces of the marrow elements (E ≈ 2 MPa). Under a penalty of 2×10⁴ N/mm³ these elements were crushed and inverted. In every failed run, 100 % of the inverted elements were marrow, all within 0.6 mm of the interface. Removing the marrow faces from the contact pairs fixed convergence: cortical and cancellous contact are unchanged, and marrow carries no load anyway. On the runs that had already converged, stiffness changed by 0.0–2.4 %, and solve time dropped 3–8×.
Your point about friction harming convergence and about choosing the penalty relative to element size and the softer material was spot on. The soft material was simply one we hadn’t thought of.
The nail–canal contact is still being tuned; I’ll report back. Thanks again for taking the time.
Glad to hear that your systematic approach worked and you’ve managed to eliminate the root cause of the convergence issue.
Indeed, contact between parts with such large stiffness differences can be problematic, even if you ensure good meshing and master-slave surface assignment (the stiffer part should be the master).
In some cases, contact interface can be modeled using element layers with compression-only material or with compression-only constraints.
FEBio has some features that could help model the tissues such as marrow in a more accurate way, including biphasic materials (and their contact).