Jacobian

Jacobian#

Lets begin with some definitions. Jacobian matrix represents deformation of an element from its intrinsic coordinates

\[ J_{IK} = \frac{\partial x}{{\partial\eta_K}} = \frac{\partial N_i}{{\partial\eta_K}} X_{iI} \]

being \(X_{iI}\) the nodal coordinates in the \(I\) dimension, the current position and X the reference position, with \(i\) in sum

import numpy as np

import numpy as np
# Define material properties
E = 200e9  # Young's modulus in Pa
nu = 0.3   # Poisson's ratio
m_dim = 2
m_nodxelem = 4
# Define element properties
num_nodes_element = 4  # Number of nodes per element
element_length = 1.0   # Length of the element
m_gp_count = 4

dNdX = np.zeros((m_gp_count, m_dim, m_nodxelem)) 
dNdrs = np.zeros((m_gp_count, m_dim, m_nodxelem)) 
# Define shape functions and their derivatives for 2D quadrilateral element
def shape_functions(xi, eta):
    dNdX_ = np.zeros((m_dim, m_nodxelem))
    N = np.array([(1-xi)*(1-eta)/4,
                  (1+xi)*(1-eta)/4,
                  (1+xi)*(1+eta)/4,
                  (1-xi)*(1+eta)/4])
    dNdX_[0,:] = np.array([-(1-eta)/4, (1-eta)/4, (1+eta)/4, -(1+eta)/4])
    dNdX_[1,:] = np.array([-(1-xi)/4, -(1+xi)/4, (1+xi)/4, (1-xi)/4])
    return N, dNdX_
    print(dNdX)
# Gauss quadrature points and weights
gauss_points = np.array([[-0.577350269, -0.577350269],
                         [ 0.577350269, -0.577350269],
                         [ 0.577350269,  0.577350269],
                         [-0.577350269,  0.577350269]])

gauss_weights = np.array([1, 1, 1, 1])

gp_count = len(gauss_points)

In case of elasticity, we can define \( B^T C B\)

for gp in range(len(gauss_points)):
    xi, eta = gauss_points[gp]
    N, dNdrs[gp] = shape_functions(xi, eta)

# Finite element strain rate calculation
def calculate_jacobian(pos):
    J = np.zeros((gp_count, 2, 2))
    for gp in range(len(gauss_points)):
        xi, eta = gauss_points[gp]
        weight = gauss_weights[gp]
        N, dNdrs[gp] = shape_functions(xi, eta)
        J[gp] = np.dot(dNdrs[gp], pos)
        detJ = np.linalg.det(J[gp])
        print("det J\n", detJ)
        invJ = np.linalg.inv(J[gp])
        print ("invJ", invJ)
        dNdX[gp] = np.dot(invJ,dNdrs[gp])
    return J

Now use it

# Example usage
x    =  np.array([[0., 0.], [0.1, 0.], [0.1, 0.1], [0., 0.1]])
velocities = np.array([[0, 0], [0, 0], [0, -1], [0, -1]])  # Example velocities at nodes
#strain_rate = calculate_strain_rate(velocities)
J = calculate_jacobian(x)
print ("Jacobian\n", J[0])

det J
 0.0025000000000000005
invJ [[ 2.00000000e+01 -4.29380248e-17]
 [ 4.29380248e-17  2.00000000e+01]]
det J
 0.0025000000000000005
invJ [[ 2.00000000e+01 -4.29380248e-17]
 [ 2.34617731e-16  2.00000000e+01]]
det J
 0.0025000000000000005
invJ [[ 2.00000000e+01 -2.34617731e-16]
 [ 2.34617731e-16  2.00000000e+01]]
det J
 0.0025000000000000005
invJ [[ 2.00000000e+01 -2.34617731e-16]
 [ 4.29380248e-17  2.00000000e+01]]
Jacobian
 [[ 5.00000000e-02  1.07345062e-19]
 [-1.07345062e-19  5.00000000e-02]]

Velocity Gradient and Strain Rate#

Velocity gradient tensor is defined as:

\[ \nabla v = \frac{dv_I}{dx_J} \]

Due shape function gradients are calculated as gauss points, We can express this as
\( \nabla v_{IJ} = \frac{dN_k}{dX_J} V_{kI} \)
This means that, for each dimension,

# Define nodal velocities (dummy data for demonstration)
vel = np.full(m_dim * m_nodxelem, 0.1)
vel[5] = vel[7] = -1.0

def velocity_gradient_tensor(dNdX, vel):
    grad_v = np.zeros((m_gp_count,m_dim, m_dim))
    for gp in range (m_gp_count):
        for I in range(m_dim): 
            for J in range(m_dim):
                for k in range(m_nodxelem): 
                    #grad_v[gp,I, J] += dNdX[gp, J, k] * vel[k * m_dim + I]
                    grad_v[gp,I, J] += dNdX[gp, J, k] * vel[k, I]
    return grad_v

def calc_str_rate (dNdX,velocities):
    str_rate = np.zeros((m_gp_count,m_dim, m_dim))
    for gp in range (m_gp_count):
        grad_v = velocity_gradient_tensor(dNdX, velocities)
        print("Velocity gradients\n" ,grad_v[0])

        str_rate[gp] = 0.5*(grad_v[0]+grad_v[0].T)
    print("strain rate:\n" ,str_rate)
    return str_rate
str_rate = calc_str_rate (dNdX,velocities)
print("strain rate:\n" ,str_rate[0])

Velocity gradients
 [[  0.   0.]
 [  0. -10.]]
Velocity gradients
 [[  0.   0.]
 [  0. -10.]]
Velocity gradients
 [[  0.   0.]
 [  0. -10.]]
Velocity gradients
 [[  0.   0.]
 [  0. -10.]]
strain rate:
 [[[  0.   0.]
  [  0. -10.]]

 [[  0.   0.]
  [  0. -10.]]

 [[  0.   0.]
  [  0. -10.]]

 [[  0.   0.]
  [  0. -10.]]]
strain rate:
 [[  0.   0.]
 [  0. -10.]]

Now, how can we compute stresses? We have several ways to do this. One of them is to assume absoute elastic behavior. In any case, we have to compute the incremental displacements from strain rate

dt = 8.0e-5
stress = np.zeros((m_gp_count,m_dim, m_dim))
strain = dt * str_rate
print (strain)

# Define material matrix for plane stress
def material_matrix():
    C = E / (1 - nu**2) * np.array([[1, nu, 0],
                                     [nu, 1, 0],
                                     [0, 0, (1 - nu) / 2]])
    return C

def calculate_strain(str_rate,dt):
    strain = np.zeros((m_gp_count,m_dim, m_dim))
    strain = dt * str_rate
    return strain
    
def calculate_stress(eps,dNdX):
    stress = np.zeros((m_gp_count,m_dim, m_dim))
    # strain = np.zeros((m_gp_count,m_dim, m_dim))
    # eps[gp] +=  str_rate * dt
  # # PLAIN STRESS
    c = E / (1.0- nu*nu)
  
  # #!!!! PLAIN STRAIN
  # #c = dom%mat_E / ((1.0+dom%mat_nu)*(1.0-2.0*dom%mat_nu))
    for gp in range(len(gauss_points)):
        stress[gp,0,0] = c * ((1.0-nu)*eps[gp,0,0] + nu*eps[gp,1,1])
        stress[gp,1,1] = c * ((1.0-nu)*eps[gp,0,0] + nu*eps[gp,1,1])
        stress[gp,0,1] = stress[gp,1,0] = c * (1.0-2*nu)*eps[gp,0,1] 
    return stress

[[[ 0.      0.    ]
  [ 0.     -0.0008]]

 [[ 0.      0.    ]
  [ 0.     -0.0008]]

 [[ 0.      0.    ]
  [ 0.     -0.0008]]

 [[ 0.      0.    ]
  [ 0.     -0.0008]]]

strain =  calculate_strain(str_rate,dt)
stress =  calculate_stress(strain,dt)
print ("strain ",strain)
print ("stress", stress)

strain  [[[ 0.      0.    ]
  [ 0.     -0.0008]]

 [[ 0.      0.    ]
  [ 0.     -0.0008]]

 [[ 0.      0.    ]
  [ 0.     -0.0008]]

 [[ 0.      0.    ]
  [ 0.     -0.0008]]]
stress [[[-52747252.74725275         0.        ]
  [        0.         -52747252.74725275]]

 [[-52747252.74725275         0.        ]
  [        0.         -52747252.74725275]]

 [[-52747252.74725275         0.        ]
  [        0.         -52747252.74725275]]

 [[-52747252.74725275         0.        ]
  [        0.         -52747252.74725275]]]

# Finite element strain rate calculation
#We can calculate with B matrix
#F = BT x sigma = [dh1/dx dh1/dy ] x [ sxx sxy]
#               = [dh2/dx dh2/dy ]   [ syx syy]
def calc_forces(stress,dNdX,J):
    forces = np.zeros((m_nodxelem,m_dim))
    B = np.zeros((m_dim, m_nodxelem))
    
    for gp in range(len(gauss_points)):
        for i in range(m_nodxelem):
            B[0, i] = dNdX[gp,0,i]
            B[1, i] = dNdX[gp,1,i]
        forces +=  np.dot(B.T,stress[gp]) *  np.linalg.det(J[gp]) * gauss_weights[gp]
        #print(forces)
    return forces 

forces = calc_forces(stress,dNdX,J)
print (forces)
    

[[ 2637362.63736264  2637362.63736264]
 [-2637362.63736264  2637362.63736264]
 [-2637362.63736264 -2637362.63736264]
 [ 2637362.63736264 -2637362.63736264]]

Note that these are the internal forces

# Finite element strain rate calculation
#We can calculate with B matrix
def calculate_stress(str_rate,dt,dNdX):
    stress = np.zeros((m_gp_count,m_dim, m_dim))
    for gp in range(len(gauss_points)):
        B = np.zeros((3, 8))
        for i in range(num_nodes_element):
            B[0, 2*i] = dNdX[gp,0,i]
            B[1, 2*i+1] = dNdX[gp,1,i]
            B[2, 2*i] = dNdX[gp,1,i]
            B[2, 2*i+1] = dNdX[gp,0,i]
        C = material_matrix()
        stress[gp] = np.dot(C,B)
        #strain_rate += np.dot(B.T, np.dot(invJ.T, np.dot(B, C))) * detJ * weight
        #print ("Jacobian", J[gp])
    return stress

Jacobian

Contents

Jacobian#

Velocity Gradient and Strain Rate#