I wrote a UEL (User Element in Abaqus) for one element and compared to a reference UEL which used standard FEM, where the results agreed satisfactorily, except the reaction force. The stress, strain, strain energy all were pretty accurate.

I am unable to understand this, how the reaction force is calculated?

I am simply restraining the element in the bottom (two nodes) and applying displacement (incrementally) at the top (two nodes).

As far as I know, the stiffness matrix is used to calculate reaction forces in these cases, the stiffness values are quite similar in my UEL and the FEM UEL. Also if this was the case, the other parameters wouldn't have had the same values.

Would be grateful if anybody could help me figure this out.

  • 1
    $\begingroup$ Welcome to Scicomp! This question is better-suited for an abaqus-specific forum or for the abaqus support team. They will know the specific internals of the abaqus solver, and for this reason, application-specific questions are generally regarded as off-topic here. If, however, your question is about the mathematics behind FEM, then rephrasing the question in this way is more likely to be on-topic here. $\endgroup$ Apr 21, 2018 at 2:32
  • $\begingroup$ Thanks for the comment! I will try to rephrase the question since it is mostly about the mathematics/mechanics behind the FEA procedure. $\endgroup$
    – Schneider
    Apr 21, 2018 at 15:43

2 Answers 2


To calculate the reaction forces at a node, Abaqus (or any structural FE code) simply sums the internal forces for all elements attached to that node. The reaction forces are the negative of that sum.

For an Abaqus user element, the internal forces for the element are returned from subroutine UEL in the RHS array. The returned stiffness matrix (Jacobian), AMATRX, is not used in the reaction force calculations.

Of course, for a simple linear element, the internal force vector will equal the negative of the product of the element stiffness matrix and the nodal displacement vector.

  • $\begingroup$ note that with one element and all displacement bc's you could return a completely wrong rhs and it will not cause any issue except giving the bad value as a reaction force. $\endgroup$
    – george
    Apr 21, 2018 at 13:13
  • $\begingroup$ @george: How would that be possible? That all other possible parameters are correctly output but only the Reaction force is wrong? I mean the RHS is used in each iteration for calculations right? Then are how are the other values correct? $\endgroup$
    – Schneider
    Apr 21, 2018 at 15:40
  • $\begingroup$ @Bill Greene: This problem involves a non-linear element, there are basically two elements defined in the UEL and I calculate the RHS as the product of Transpose of B matrix with the Stress vector. $\endgroup$
    – Schneider
    Apr 21, 2018 at 15:51
  • $\begingroup$ with displacement boundary conditions on every node in the model the solver doesn't do anything with the rhs vector. You really should be using some force conditions to test your uel. (I would recomend testing in a model mixed with some standard elements as well) $\endgroup$
    – george
    Apr 21, 2018 at 17:33
  • $\begingroup$ Thanks for your comments! I did run a test with more than one element, there the UEL failed to show convergence! It didn't even converge for the first increment! The strange thing is RHS values provided by my UEL and the FEM one is almost similar. $\endgroup$
    – Schneider
    Apr 21, 2018 at 17:39

Once the solution of the problem is known, i.e. you know displacement vector, to calculate reaction/internal forces an integral is evaluated

\begin{equation} \mathbf{f}^\textrm{int} = \sum_{e=1}^{n_e} \int_\Omega \mathbf{B}^\textrm{T}({\boldsymbol\sigma}(\boldsymbol\varepsilon)) \textrm{d}\Omega = \sum_{i=1}^{n_e} \sum_{i=1}^{n_g} j_i w_i \left( \mathbf{B}_i^\textrm{T}{\boldsymbol\sigma}(\boldsymbol\varepsilon_i) \right) \end{equation} where $j_i$ and $w_i$ is jacobian and integration weight respectively, $\mathbf{B}$ is differential operator evaluated at integration point $i$, such that \begin{equation} \boldsymbol\varepsilon_i = \mathbf{B}_i\mathbf{u}^e \end{equation} where $\mathbf{u}^e$ is vector of nodal degrees of freedom on element $e$.

Elements of the vector $\mathbf{f}^\textrm{int}$ which are at constrained degrees of freedom are reaction forces.

Equations above are general, apply to the linear and nonlinear problem. To make it work you need to provide a physical equation for example in UMAT, i.e. \begin{equation} {\boldsymbol\sigma} = {\boldsymbol\sigma}( \boldsymbol\varepsilon ) \end{equation}

  • $\begingroup$ Just wanted to verify something I discussed in the previous answer as well. Suppose I am using displacement boundary conditions only. In my model, One Element, two nodes are fixed and other two nodes have a fixed displacement(applied incrementally), so the only output variable to be compared(with another reference code with same boundary conditions) is the internal force(or reaction force),since all the displacements are known, every other parameter such as Strain-Stress-Strain Energy would anyways be the same? $\endgroup$
    – Schneider
    Apr 23, 2018 at 0:32
  • $\begingroup$ Yes. You can check if $K u - f^\textrm{int} = 0$, and energy which should be $E=\frac{1}{2} u f^\textrm{int}$. $\endgroup$
    – likask
    Apr 23, 2018 at 7:23
  • $\begingroup$ Would you kindly point me to a derivation of $f^{int}$? $\endgroup$
    – Olumide
    Jan 23, 2022 at 4:49
  • $\begingroup$ To be specific, the part that I don't understand is $\int_{\Omega} B^T \sigma d\Omega$. Also, does this formula cope with/ignore rotations? $\endgroup$
    – Olumide
    Jan 23, 2022 at 5:48

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.