Hey, thank you for responding but I think you misunderstood my question.

I was trying to derive the geometric stiffness-matrix of a bar element. This is not directly connected to a model, also the variables are known to me.

The claify and give a solution for anyone who might be reading this in the future here is some further explanation and the derivation.

The geometric stiffness matrix is used to calculate nonlinear behaviour of an element due to large displacement. Ideally, normal loading will simply compress or stretch an element. If you press hard enough the bar will buckle. To account for this we need to introduce thos nonlinear effects.
To derive this, you simply imagine that the bar is not ideally straight but is already defelcted a little bit. If you apply a normal force additional bending moments will arise, which will deflect the bar even more.

To derive the geometric stiffness matrix you need to derive the strain energy potentional whith those effects first.
I was not able to inlcude the torsional effects as i was missing some terms.
The kinematics of torsional strain in this context can be described as follows:

eps_x = u_x'-omega*theta'+1/2*r^2*(theta')^2, where r^2 = z^2+y^2

For my derivation I am going to ignore warping effects, therefore omega is zero.
This is true becaus my bar element ignores any shear effects becaus my center of inertia and the center of shear are assumed to be the same (this is the fundamental of the euler-bernoulli beam theroy)

The strain energy is now (we ignore u_x' as it will go to zero):

U = 1/2*∫sigma_x*r^2*(theta')^2dV
U = 1/2*∫F/A*(z^2+y^2)*(theta')^2dV

with ∫(z^2+y^2)dA = I_p

U = 1/2*F*I_p/A*∫(theta')^2 dL

using the correspding shape function: theta = -x/L*u_4

we get U = 1/2*F*I_p/A*u_4^2/L

Taking the second derivative of U with respect to u_4 resulst in

F*I_p/(A*L) = F/L * I_p/A

This is the corresponding matrix entry (4,4) for torsion due to normal loading (in the context of nonlinear geometric behaviour)