Attached below is some code I wrote to solve a basic problem: finding the node displacements of a cube with two vertices constrained (vertices 6 and 7 with coordinates [1,1,1], [1,-1,1]), and the rest of the vertices loaded with point loads in the y direction. What I expect is for the displacement at vertex 0, with coordinates [-1,-1,-1], to be strictly inside the plane spanned by [0,1,0], [1,0,1] because of the symmetry of the problem. What I get is a displacement that is not on that plane, with components [-0.00289863, 0.07050896, -0.02001042].
The way I solve the problem is by calculating the stiffness matrix with Gaussian integration. Then I remove the rows and columns in the matrix that correspond to the constrained vertices (actually I set them to 0, and the element on the diagonal to 1, but it should give the same results, right?). And then I solve the equation kd = f where k is the stiffness matrix, d is the displacement vector (24 dimensions) and f is the nodal forces vector.
I tested the shape function derivatives (test code included) by numerically calculating the derivatives and comparing with the ones returned by the function. I also tested the strain matrix by calculating the strain numerically. These are not to blame.
Please, help me find my mistake.
Edit: The problem is not well-formed, because there is nothing preventing the cube from rotating around the vertices 6-7 axis. However, why is the determinant of the stiffness matrix not 0? Shouldn't the matrix be singular?
import numpy as np
from math import sqrt
_samplePoints = [[0.0],
[-sqrt(1./3),sqrt(1./3)],
[-sqrt(0.6),0.0,sqrt(0.6)],
[-0.861136, -0.339981, 0.339981, 0.861136],
[-0.906180, -0.538469, 0, 0.538469, 0.906180],
[-0.932470, -0.661209, -0.238619, 0.238619, 0.661209, 0.932470]]
_coeffs = [[2.0],
[1.0,1.0],
[0.555556,0.888889,0.555556],
[0.347855, 0.652145, 0.652145, 0.347855],
[0.236927, 0.478629, 0.568889,0.478629, 0.236927],
[0.171324, 0.360762, 0.467914, 0.467914, 0.360762, 0.171324]]
def GetGaussQuadrature(numPoints):
return _samplePoints[numPoints-1],_coeffs[numPoints-1]
_naturalCoordinatesVertices = np.array([[-1,-1,-1],[-1,1,-1],[-1,1,1],[-1,-1,1],[1,-1,-1],[1,1,-1],[1,1,1],[1,-1,1]])
def ShapeFuncVec(ksiEtaZeta):
return 0.125 * np.prod(1 + ksiEtaZeta.reshape(1,3) * _naturalCoordinatesVertices,axis=1)
def DShapeFuncVec_DNaturalCoordinates(ksiEtaZeta):
return 0.125 * _naturalCoordinatesVertices[:,[0,1,2]]*(1+_naturalCoordinatesVertices[:,[1,2,0]]*ksiEtaZeta[[1,2,0]])*(1+_naturalCoordinatesVertices[:,[2,0,1]]*ksiEtaZeta[[2,0,1]])
def CalcJacobian(ksiEtaZeta,vertices):
return DShapeFuncVec_DNaturalCoordinates(ksiEtaZeta).transpose().dot(vertices)
def DShapeFuncVec_DXYZ(ksiEtaZeta,vertices):
dShapeFuncVec_DKsiEtaZeta = DShapeFuncVec_DNaturalCoordinates(ksiEtaZeta)
J = dShapeFuncVec_DKsiEtaZeta.transpose().dot(vertices)
Jinv = np.linalg.inv(J)
dShapeFuncVec_dXYZ = dShapeFuncVec_DKsiEtaZeta.dot(Jinv.transpose())
return dShapeFuncVec_dXYZ, J
def CalcStrainMatrix(ksiEtaZeta, vertices):
dN_dxyz,J = DShapeFuncVec_DXYZ(ksiEtaZeta,vertices)
assert(np.linalg.det(J)>0)
B = np.zeros((6,8*3))
for i in range(8):
Bi = np.zeros((6,3))
dNi_dxyz = dN_dxyz[i,:]
Bi[[0,1,2],[0,1,2]] = dNi_dxyz
Bi[[3,3,4,4,5,5],[1,2,0,2,0,1]] = dNi_dxyz[[2,1,2,0,1,0]]
B[:,(i*3):(i*3+3)] = Bi
return B,J
def GetMaterialTensor():
youngsModulus = 50.0
poissonRatio = 0.4
denom = ((1-2*poissonRatio)*(1+poissonRatio))
c11 = youngsModulus * (1-poissonRatio) / denom
c12 = youngsModulus * poissonRatio / denom
materialTensor = np.zeros((6,6))
materialTensor[:3,:3] = c12
materialTensor[[0,1,2],[0,1,2]] = c11
materialTensor[[3,4,5],[3,4,5]] = (c11-c12)/2
return materialTensor
materialTensor = GetMaterialTensor()
import itertools
def CalcStiffnessMatrix(vertices, materialTensor):
n=3
#perform gauss integration with n points
samplePoints,coeffs = GetGaussQuadrature(n)
Ke = np.zeros((24,24),dtype=float)
count=0
for i,j,k in itertools.product(range(n),repeat=3):
count+=1
ksiEtaZeta = np.array([samplePoints[i],samplePoints[j],samplePoints[k]])
coeff = coeffs[i]*coeffs[j]*coeffs[k]
B,J = CalcStrainMatrix(ksiEtaZeta,vertices)
Ke += np.linalg.det(J)*(B.transpose().dot(materialTensor).dot(B)) * coeff
assert(count == n**3)
return Ke
def TestDerivatives():
ksiEtaZeta_sample = np.array([0.4,0.3,0.2])
print(ShapeFuncVec(ksiEtaZeta_sample))
print(sum(ShapeFuncVec(ksiEtaZeta_sample)))
print(DShapeFuncVec_DNaturalCoordinates(ksiEtaZeta_sample))
print(CalcJacobian(ksiEtaZeta_sample,vertices))
dksiEtaZeta = np.array([0.000001,-0.00002,0.000003])
xyz_sample = ShapeFuncVec(ksiEtaZeta_sample).dot(vertices)
xyz_sample2 = ShapeFuncVec(ksiEtaZeta_sample + dksiEtaZeta).dot(vertices)
print(xyz_sample,xyz_sample2)
dxyz_sample = xyz_sample2 - xyz_sample
dxyz_sample_estimated = CalcJacobian(ksiEtaZeta_sample,vertices).transpose().dot(dksiEtaZeta)
print(dxyz_sample)
print(dxyz_sample_estimated)
assert(np.max(np.abs(dxyz_sample-dxyz_sample_estimated))/np.max(np.abs(dxyz_sample)) < 1e-5)
x1 = ShapeFuncVec(ksiEtaZeta_sample + dksiEtaZeta) - ShapeFuncVec(ksiEtaZeta_sample)
x2 = DShapeFuncVec_DXYZ(ksiEtaZeta_sample,vertices)[0].dot(dxyz_sample)
print(x1)
print(x2)
assert(np.max(np.abs(x1-x2))/np.max(np.abs(x1)) < 1e-5)
def TestStrain():
ksiEtaZeta = np.random.rand(3) *2-1
B,J = CalcStrainMatrix(ksiEtaZeta,vertices)
# print(J)
#displacements = np.array([[0.1,0.5,-0.1]] + [[0,0,0]]*7)
displacements = 0.2 * np.random.rand(8,3)-0.1
displacements_P = displacements.transpose().dot(ShapeFuncVec(ksiEtaZeta))
strain_fromB = B.dot(displacements.flatten())
strain_fromNumDiff = np.zeros((3,3))
dKezMag = 0.01
dKez = np.eye(3) * dKezMag
kezPdkez = np.vstack([ksiEtaZeta]*3) + dKez
for i in range(3):
thisDisplacement = displacements.transpose().dot(ShapeFuncVec(kezPdkez[i,:]))
dDisplacement = thisDisplacement - displacements_P
strain_fromNumDiff[:,i] = dDisplacement / dKezMag
strain_fromNumDiff_vec = np.array([strain_fromNumDiff[0,0],
strain_fromNumDiff[1,1],
strain_fromNumDiff[2,2],
strain_fromNumDiff[1,2]+strain_fromNumDiff[2,1],
strain_fromNumDiff[0,2]+strain_fromNumDiff[2,0],
strain_fromNumDiff[1,0]+strain_fromNumDiff[0,1]])
print(strain_fromB)
print(strain_fromNumDiff)
print(strain_fromNumDiff_vec)
assert(np.max(np.abs(strain_fromB-strain_fromNumDiff_vec)) < 1e-6)
print(np.max(np.abs(strain_fromB-strain_fromNumDiff_vec)))
vertices = _naturalCoordinatesVertices.copy()
Ke = CalcStiffnessMatrix(vertices,materialTensor)
constrainedVertices = np.vstack([[6,7]]*3)
constrainedVertices *= 3
constrainedVertices += np.arange(3).reshape(3,-1)
Ke[constrainedVertices,:]=0
Ke[:,constrainedVertices]=0
Ke[constrainedVertices,constrainedVertices]=1
print("det(Ke) = ",np.linalg.det(Ke))
print("eig(Ke) = ",np.linalg.eig(Ke)[0])
forces = np.zeros((8,3),dtype=float)
forces[:6,1]=0.1
forces=forces.flatten()
displacements = np.linalg.solve(Ke,forces).reshape(-1,3)
print()
print("result nodal displacements = ")
print(displacements)
print()
print("Tests:")
TestDerivatives()
TestStrain()
Output:
det(Ke) = 12856964.398889132
eig(Ke) = [1.98207971e+02 5.97329352e+01 5.36174927e+01 4.56626826e+01
3.28979711e+01 3.02076122e+01 2.88176063e+01 2.70455110e+01
2.10812642e+01 1.73934001e+01 1.69113782e+01 1.29241062e+01
9.90047362e+00 9.82465058e+00 4.45269904e+00 9.04339639e-15
1.37800889e+00 1.37375076e+00 1.00000000e+00 1.00000000e+00
1.00000000e+00 1.00000000e+00 1.00000000e+00 1.00000000e+00]
result nodal displacements =
[[-0.00289863 0.07050896 -0.02001042]
[ 0.02001042 0.07050896 0.00289863]
[ 0.01565451 0.04365447 0.01244406]
[-0.01565451 0.04365447 -0.02955586]
[-0.01244406 0.04365447 -0.01565451]
[ 0.02955586 0.04365447 0.01565451]
[ 0. 0. 0. ]
[ 0. 0. 0. ]]
Tests:
[0.042 0.078 0.117 0.063 0.098 0.182 0.273 0.147]
1.0
[[-0.07 -0.06 -0.0525]
[-0.13 0.06 -0.0975]
[-0.195 0.09 0.0975]
[-0.105 -0.09 0.0525]
[ 0.07 -0.14 -0.1225]
[ 0.13 0.14 -0.2275]
[ 0.195 0.21 0.2275]
[ 0.105 -0.21 0.1225]]
[[ 1.00000000e+00 -1.38777878e-17 1.38777878e-17]
[ 0.00000000e+00 1.00000000e+00 0.00000000e+00]
[-1.38777878e-17 0.00000000e+00 1.00000000e+00]]
[0.4 0.3 0.2] [0.400001 0.29998 0.200003]
[ 1.e-06 -2.e-05 3.e-06]
[ 1.e-06 -2.e-05 3.e-06]
[ 9.72493763e-07 -1.62249301e-06 -1.70250199e-06 1.85250124e-06
2.50249124e-06 -3.35249199e-06 -3.32251301e-06 4.67251376e-06]
[ 9.7250e-07 -1.6225e-06 -1.7025e-06 1.8525e-06 2.5025e-06 -3.3525e-06
-3.3225e-06 4.6725e-06]
[-0.01654885 -0.00454239 0.00107933 0.01714707 -0.00310323 -0.09393567]
[[-0.01654885 -0.03663831 -0.00499938]
[-0.05729736 -0.00454239 0.00809046]
[ 0.00189615 0.00905662 0.00107933]]
[-0.01654885 -0.00454239 0.00107933 0.01714707 -0.00310323 -0.09393567]
7.927686285214008e-16
