I need to add to an existing FEM solver some embedded reinforcement element. This would give me the possibility to model/solve concrete structure (reinforced with steel rebar) taking into account the bond slip.

What is allready implemented (no bond-slip)


In the current model, the steel rebar are taken into account with the assumption of a perfect bond with the concrete element. The local stiffness matrix of this element is the following: $$K_{local} = \frac{EA}{L} \begin{bmatrix} 1& 0& -1 & 0\\ 0& 0& 0&0 \\ -1&0 & 1 & 0\\ 0& 0& 0 & 0 \end{bmatrix}$$ It's possible to have the global stiffness matrix with a little rotating matrix: $$K_{global} = T^{^{T}}K_{local} T$$

What I plan to do


I will use the link element develloped by Ngo and Scordelis (1967).The link element may be conceptualized as two orthogonal springs, linking the reinforced concrete element and the discrete reinforcement element.

As we can see in the upper figure, the steel point j is connected to the concrete point i. We can define the following vector (corresponding dofs): $$\bar{r} = \begin{bmatrix} r_{i,x}\\ r_{i,y}\\ r_{j,x}\\ r_{j,y}\\ \end{bmatrix}$$

The tangeantial/radial displacement vector can be computed by: $$\begin{bmatrix} \Delta _{t}\\ \Delta _{r}\\ \end{bmatrix} = \begin{bmatrix} -cos(\phi) & -sin(\phi) & cos(\phi) &sin(\phi) \\ sin(\phi)&-cos(\phi) & -sin(\phi)& cos(\phi) \end{bmatrix} \bar{r}$$ The corresponding Force vector (in the two springs) can be computed by:

$$ \begin{bmatrix} F_{t}\\ F_{r} \end{bmatrix}=A\begin{bmatrix} k_{t} & 0 \\ 0 & k_{r} \end{bmatrix}\begin{bmatrix} \Delta _{t}\\ \Delta _{r}\\ \end{bmatrix}$$

With Kt corresponding to the stiffness of the interface, Kr set to a big value (typically 100 Kt) to neglect dowel action, and A being the circunferential area of the rebar.

Finally, we can rotate back the forces in the global (x-y) coordinate system:

$$\bar{F} = \begin{bmatrix} F_{i,x}\\ F_{i,y}\\ F_{j,x}\\ F_{j,y}\\ \end{bmatrix} = T^{T}\begin{bmatrix} F_{t}\\ F_{r} \end{bmatrix}$$

It is possible to remark that the final (assembled) equation is under this form:

$$\bar{F} = [K] \bar{r}$$

Or more generally: $$\begin{bmatrix} \bar{F_{conc}}\\ \bar{F_{steel}} \end{bmatrix} = \begin{bmatrix} [K_{cc}] & [K_{cs}]\\ [K_{cs}] & [K_{ss}] \end{bmatrix} \begin{bmatrix} \bar{d_{conc}}\\ \bar{d_{steel}} \end{bmatrix}$$

It's possible to apply a static condensation of the steel dofs: $$\bar{F_{c}^{*}} = \bar{F_{conc}} - [K_{cs}][K_{ss}]^{-1}\bar{F_{steel}}$$ $$[K_{c}^{*}]= [K_{cc}] - [K_{cs}][K_{ss}]^{-1}[K_{cs}]$$ With this equivalent formulation, the system can be seen under this form: $$\bar{F_{c}^{*}} = [K_{c}^{*}]\bar{d_{conc}}$$

As you can see, only the degrees of freedom of the concrete are involved. The big advantage is that we can now simply add the inside Force vector coming from the steel rebar, the concrete element and the bonding link. We can also sum the tangeant stiffness matrix of the concrete, steel rebar, and bonding link (and apply after a standart Newton/Rapson step).

So where is my problem?

I will take a very simple example to express my problem:

Link_element As you can see:

  • Horizontal bar $\phi = 0 $
  • Circunferential area of the rebar : $0.015 m^{2}$
  • $k_{t} = 10000 [Mpa]$ for the tangeantial imposed slip (0.5mm)
  • $k_{r} = 300000 [Mpa]$ to avoid radial displacement

We get quickly: $$\begin{bmatrix} \Delta _{t}\\ \Delta _{r}\\ \end{bmatrix} = \begin{bmatrix} 0.0005\\ 0\\ \end{bmatrix} = \begin{bmatrix} -1 & 0 &1 &0 \\ 0&-1 & 0& 1 \end{bmatrix} \begin{bmatrix} 0\\ 0\\ 0.0005\\ 0 \end{bmatrix}$$

Searching the force we get: $$ \begin{bmatrix} F_{t}\\ F_{r} \end{bmatrix}=\begin{bmatrix} 0.075\\ 0 \end{bmatrix}=0.015\begin{bmatrix} 10000 & 0 \\ 0 & 300000 \end{bmatrix}\begin{bmatrix} 0.0005\\ 0\\ \end{bmatrix}$$

Applying the rotating matrix: $$\bar{F} = \begin{bmatrix} F_{i,x}\\ F_{i,y}\\ F_{j,x}\\ F_{j,y}\\ \end{bmatrix} = \begin{bmatrix} -0.075\\ 0\\ 0.075\\ 0\\ \end{bmatrix} =T^{T}\begin{bmatrix} 0.075\\ 0 \end{bmatrix}$$

Which seems to be physically correct (opposite forces in concrete/steel elements, correct value...).

The corresponding stiffness matrix is the following: $$[K] = \begin{bmatrix} 150 & 0 & -150 & 0 \\ 0& 4500 & 0 & -4500 \\ -150 & 0 & 150 & 0\\ 0& -4500 & 0 & 4500 \end{bmatrix}$$

Applying the previous condensed equation, we can build: $$\bar{F_{c}^{*}} = \begin{bmatrix} 0\\ 0 \end{bmatrix}$$ $$[K_{c}^{*}]= \begin{bmatrix} 0 & 0\\ 0 & 0 \end{bmatrix}$$


  1. The equivalent bonding load vector is full of zeros! What is the meaning of this? Even with a consequent bond-slip, no force are applied in the condensed concrete model?
  2. The equivalent tangeant stiffness (bonding link) is full of zeros! This is a big problem, since it will raise a singular matrix error in the next step of the Newton/Rapson iteration!

Is there any physical explanation for this? For sure, I am misunderstanding a concept, but this is a real problem in my implementation purpose. Maybe it's a question of boundary condition... I am really not sure.


1 Answer 1


It looks like you statically condense dofs node by node. To make it work this should be done at least element by element. Your matrix $K_{cc}$ is 8 by 8, whereas matrix $K_{ss}$ is 4 by 4.

This is the assumption that you fibres are short or micro-cracking make the response of long fibre limited to short length. Static condensation, which you presenting has an advantage, that you can implement that on element level and in well-implemented cold be numerically efficient and be easily implemented in standard FE code, for example using UMAT matrix in ABAQUS.

However, you should consider applying condensation on the system level, to get realistic redistribution of stresses for long reinforced fibres. In that case, static condensation is applied for each fibre independently. Instead of that, you can consider using multiplicative/additive block solver or block solver with Schur complement. In fact, you don't have to implement this by yourself and use PETSc field split preconditioner, see http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/PC/PCFIELDSPLIT.html

Note bar does not have to be restricted to edge of solid/concrete element, you have to express displacements at the end of the bar by linear combination of dofs of concrete element. In principle method will be the same.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.