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:
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}$$
And...?
- 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?
- 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.
