Каждый раз в цикле вы можете добавлять в свой список, а затем обращаться ко всем значениям.
result_list = []
stiffness_matrix_total = np.zeros((2*nnode,2*nnode), dtype=np.float)
for i in range(nele):
ele_stiff_mat = np.zeros((2*nnode,2*nnode), dtype=np.float)
a = np.array([[cos[i,0]**2,sin[i,0]*cos[i,0],-(cos[i,0]**2),-cos[i,0]*sin[i,0]],[cos[i,0]*sin[i,0],sin[i,0]**2,-cos[i,0]*sin[i,0],-(sin[i,0]**2)],[-cos[i,0]**2,-sin[i,0]*cos[i,0],cos[i,0]**2,cos[i,0]*sin[i,0]],[-cos[i,0]*sin[i,0],-(sin[i,0]**2),cos[i,0]*sin[i,0],sin[i,0]**2]])
result_list.append(a)
print (result_list)