Question:

# -*- coding: utf-8 -*- """2023MEM1062_Fem_Code_Project.ipynb Automatically generated by Colab. Original file is located at https://colab.research.google.com/drive/1yz_EdoPO4NYHS9E5NxXGJcSA1LhNnYar Name :- Manas Kishor Thakur Branch :- Mechanical Engineering Mtech Mechanics and Design Roll.No. :- 2023MEM1062 Subject :- Finite Element Methods In Engineering-Term Project ### Importing The Liabraries ### """ import numpy as np import pandas as pd import matplotlib.pyplot as plt import sympy as sp from matplotlib.colors import Normalize from matplotlib.tri import Triangulation from matplotlib.patches import Polygon from matplotlib.collections import PatchCollection import pandas as pd """Defining the function to compute to compute jacobian matric and further which is used to cumpute the shape functions wrt x and y derivatives""" def gauss_points(option): # Define a function gauss_points that takes an argument option. if option == 1: loc = np.array([[-1 / np.sqrt(3), -1 / np.sqrt(3)], # Define the coordinates of Gauss points [1 / np.sqrt(3), -1 / np.sqrt(3)], [1 / np.sqrt(3), 1 / np.sqrt(3)], [-1 / np.sqrt(3),1 / np.sqrt(3)]]) wt = np.ones((4,1)) # Initialize the weights for full integration, s etting all to 1. else: loc =np.zeros((1,2)) # Define the coordinates of Gauss points for partial integration, setting them to zero. wt = 4 # Set the weight for partial integration to 4. return wt,loc # Return the weights and locations. def Jacobian(nodeCoord,der_n): # Define a function Jacobian that takes node coordinates and derivatives as inputs. J = (nodeCoord.T) @ der_n dn_zhi_eta = (der_n).T J_Transpose = J.T dn_x_y = np.linalg.solve(J_Transpose,dn_zhi_eta) # Solve for the derivative of x and y with respect to zhi and eta der_x_y = (dn_x_y).T return J, der_x_y # Return the Jacobian matrix and the derivative of x and y. """Defining the function to compute elemental stiffness matrix and then converting the matrix obtained from each element to the global stiffness matrix""" def Stiffness_matrix(x_y_deof,num_Elements,elements,nP,Node_List,E_matrix,h): K = np.zeros((x_y_deof,x_y_deof)) # Initialize the stiffness matrix with zeros. weightsonpoints, gauss_point = gauss_points(1) # Get the Gauss points and weights. for e in range(num_Elements): id = elements[e,:] # Get the node IDs of the current element. eDof = np.zeros((8,1)) # Initialize element degrees of freedom. eDof[0:4,0] = id # Assign the node IDs to the first half of element degrees of freedom. eDof[4:8,0] = id + nP # Assign the node IDs with added node number to the second half. eDof = eDof.flatten() x_y_deof = id.size # loop for Gauss point for q in range(weightsonpoints.size): GaussPoint = gauss_point[q,:] zhi = GaussPoint[0] eta = GaussPoint[1] # shape functions and derivatives shape_mat,der_n = N_Shape(zhi,eta) # Jacobian matrix, inverse of Jacobian J,der_x_y = Jacobian(Node_List[id-1,:],der_n) # B_mat matrix (Linear strain - displacement matrix) B_mat = np.zeros((3,2*x_y_deof)) B_mat[0,0:x_y_deof] = np.transpose(der_x_y[:,0]) B_mat[1,x_y_deof:(2*x_y_deof)] = np.transpose(der_x_y[:,1]) B_mat[2,0:x_y_deof] = np.transpose(der_x_y[:,1]) B_mat[2,x_y_deof:(2*x_y_deof)] = np.transpose(der_x_y[:,0]) # stiffness matrix B_mat_Transpose = np.transpose(B_mat) J_det = np.linalg.det(J) K_elemental = np.matmul(np.matmul(B_mat_Transpose,E_matrix),B_mat)*h*J_det*weightsonpoints[q] for ii in range(np.size(K_elemental,0)): row = int(eDof[ii])-1 for jj in range(np.size(K_elemental,1)): col = int(eDof[jj])-1 K[row,col] = K[row,col] + K_elemental[ii,jj] return K def N_Shape(zhi,eta): shape_mat = np.zeros((4,1)) # Initialize the shape matrix. shape_mat[:,0] = 1/4*np.array([(1-zhi)*(1-eta),(1+zhi)*(1-eta),(1+zhi)*(1+eta),(1-zhi)*(1+eta)]) der_n = np.zeros((4,2)) der_n[0,:] = 1/4*np.array([-(1-eta),-(1-zhi)]) # Compute derivatives for the shape function. der_n[1,:] = 1/4*np.array([1-eta,-(1+zhi)]) der_n[2,:] = 1/4*np.array([1+eta,1+zhi]) der_n[3,:] = 1/4*np.array([-(1+eta),1-zhi]) return shape_mat, der_n # Return the shape matrix and the derivatives of shape functions. E = 2e11 h = 0.01 nue = 0.3 S =E/(1-nue**2) E_matrix = S*np.array([[1, nue, 0], [nue, 1, 0], [0, 0, (1-nue)/2]]) """Plotting Mesh""" def Plotting_Mesh(Node_List, elements, nel): nnel = np.size(elements, 1) X = np.zeros((nnel, nel)) Y = np.zeros((nnel, nel)) for iel in range(nel): for i in range(nnel): ndi = elements[iel, i] X[i, iel] = Node_List[ndi-1, 0] Y[i, iel] = Node_List[ndi-1, 1] plt.figure(figsize=(10, 10)) # Increased figure size plt.axis('equal') plt.fill(X, Y, facecolor='grey', edgecolor='black', linewidth=2) # Changed edge color and increased line width plt.xticks(fontsize=16) # Increased font size for x-axis ticks plt.yticks(fontsize=16) # Increased font size for y-axis ticks plt.show() N_list = pd.read_csv('NL_K.csv', header=None) N_list = N_list.dropna(axis=0) N_list1 = N_list[N_list[0].isin(["NODE"])] N_list = pd.concat([N_list, N_list1, N_list1]).drop_duplicates(keep=False) nP = N_list.count(axis=0) nP = nP[0] x_y_deof = 2*nP Node_List = N_list.iloc[:,1:3] Node_List = Node_List.reset_index() Node_List = Node_List.iloc[:,1:3] Node_List = Node_List.to_numpy(dtype=np.float32) Nodes = pd.read_csv('EL_K.csv', header=None) Nodes = Nodes.dropna(axis=0) Nodes = Nodes[[6,7,8,9]] elements = Nodes[[6,7,8,9]] num_Elements = elements.count(axis=0) num_Elements = num_Elements[6] elements.head() elements = elements.astype(int) elements = elements.to_numpy(dtype=np.int32) elements.shape Plotting_Mesh(Node_List,elements,num_Elements) """Calculating Stiffness Matrix""" K = Stiffness_matrix(x_y_deof,num_Elements,elements,nP,Node_List,E_matrix,h) K.shape print(K) """calculating eigen values and eigen vectors to verify stiffness matrix""" eigv, eigvec=np.linalg.eig(K) eigv NL=pd.read_csv('NL_U.csv',header=None).values # Importing the Node List EL=pd.read_csv('EL_U.csv',header=None).values # Import the Element List L= np.array([56,59,62,65,68,71,74]) # Defining array of nodes to apply force """Calculating Force Vector""" def forcing(g2,g3,force): zhi, eta =sp.symbols('zhi eta') N1=(1/4)*(zhi-1)*(eta-1) N2=(-1/4)*(zhi+1)*(eta-1) N3=(1/4)*(zhi+1)*(eta+1) N4=(-1/4)*(zhi-1)*(eta+1) dN2e1=sp.diff(N2,eta) dN3e1=sp.diff(N3,eta) JL11=(dN2e1*(NL[g2,0])) JL12=(dN3e1*(NL[g3,0])) JL21=(dN2e1*(NL[g2,1])) JL22=(dN3e1*(NL[g3,1])) JL1a=(JL11+JL12)**2 JL2a=(JL21+JL22)**2 JL1q=sp.sqrt(JL1a+JL2a) JL=JL1q.subs(zhi,1) force_vec=np.zeros([2*len(NL),1]) force_vec[2*g2]=JL*force force_vec[2*g3]=JL*force return force_vec global_force=np.zeros([150,1]) force=5*(10**7) for i in range(len(L)-1): f1=forcing(L[i],L[i+1],force) global_force=f1+global_force global_force fix=np.array([0,2,4,6,8,10,84,86,88,73,75,77,79,81,83,109,111,113]) #Fix nodes k_reduced_1=np.delete(K,fix,axis=0) #Reduced Stiffness Matrix k_reduced=np.delete(k_reduced_1,fix,axis=1) k_reduced.shape force_red=np.delete(global_force,fix,axis=0) force_red.shape k_inv=np.linalg.inv(k_reduced) U = np.dot(k_inv,force_red) #calculating displacement U U1=U.flatten() U2=U1.tolist() print(U1) def insert_zeros(U, fix): # adding zeros to the displacement vectors where reduced the matrix new_list = U.copy() for index in fix: new_list.insert(index, 0) # Insert zero at the specified index return new_list new_list = insert_zeros(U2, fix) print(new_list) len(new_list) u_deform=[] #Storing array of deformed in u direction v_deform=[] #Storing array of deformed in v direction for i in range(len(new_list)): if i% 2 == 0: u_deform.append(new_list[i]) else: v_deform.append(new_list[i]) print(u_deform) print(v_deform) u_deform1=np.array([u_deform]).T v_deform1=np.array([v_deform]).T deformed=np.column_stack([u_deform1,v_deform1]) print(deformed) deformed.shape max_u_deform = max(u_deform) min_u_deform = min(u_deform) max_v_deform = max(v_deform) min_v_deform = min(v_deform) print("Maximum deformation in u direction:", max_u_deform) print("Minimum deformation in u direction:", min_u_deform) print("Maximum deformation in v direction:", max_v_deform) print("Minimum deformation in v direction:", min_v_deform) """Plotting the u and v displacements""" deformed1=NL+deformed deformed1 plt.plot(deformed1[:,0],deformed1[:,1],'*',color='red') plt.show() plt.plot(NL[:,0],NL[:,1],'*',color='red') plt.show() def Plotting_Mesh(Node_List, elements, nel): nnel = np.size(elements, 1) X = np.zeros((nnel, nel)) Y = np.zeros((nnel, nel)) for iel in range(nel): for i in range(nnel): ndi = elements[iel, i] X[i, iel] = Node_List[ndi-1, 0] Y[i, iel] = Node_List[ndi-1, 1] plt.figure(figsize=(10, 10)) # Increased figure size plt.axis('equal') plt.fill(X, Y, facecolor='grey', edgecolor='black', linewidth=2) # Changed edge color and increased line width plt.xticks(fontsize=16) # Increased font siz e for x-axis ticks plt.yticks(fontsize=16) # Increased font size for y-axis ticks plt.show() Plotting_Mesh(deformed1,elements,num_Elements) def strs_strn(L): u1,u2,u3,u4=deformed1[L[0],0],deformed1[L[1],0],deformed1[L[2],0],deformed1[L[3],0] v1,v2,v3,v4=deformed1[L[0],1],deformed1[L[1],1],deformed1[L[2],1],deformed1[L[3],1] x1,x2,x3,x4=NL[L[0],0],NL[L[1],0],NL[L[2],0],NL[L[3],0] y1,y2,y3,y4=NL[L[0],1],NL[L[1],1],NL[L[2],1],NL[L[3],1] zhi, eta =sp.symbols('zhi eta') N1=(1/4)*(zhi-1)*(eta-1) N2=(-1/4)*(zhi+1)*(eta-1) N3=(1/4)*(zhi+1)*(eta+1) N4=(-1/4)*(zhi-1)*(eta+1) x_z = (N1*x1)+(N2*x2)+(N3*x3)+(N4*x4) y_z = (N1*y1)+(N2*y2)+(N3*y3)+(N4*y4) dxz = sp.diff(x_z,zhi) dxe=sp.diff(x_z,eta) dyz=sp.diff(y_z,zhi) dye=sp.diff(y_z,eta) Ja=sp.Matrix([[dxz,dyz],[dxe,dye]]) J=Ja.det() Jinv =Ja.inv() dN1x=Jinv[0]*(sp.diff(N1,zhi)) + Jinv[1]*(sp.diff(N1,eta)) dN2x=Jinv[0]*(sp.diff(N2,zhi)) + Jinv[1]*(sp.diff(N2,eta)) dN3x=Jinv[0]*(sp.diff(N3,zhi)) + Jinv[1]*(sp.diff(N3,eta)) dN4x=Jinv[0]*(sp.diff(N4,zhi)) + Jinv[1]*(sp.diff(N4,eta)) dN1y=Jinv[2]*(sp.diff(N1,zhi)) + Jinv[3]*(sp.diff(N1,eta)) dN2y=Jinv[2]*(sp.diff(N2,zhi)) + Jinv[3]*(sp.diff(N2,eta)) dN3y=Jinv[2]*(sp.diff(N3,zhi)) + Jinv[3]*(sp.diff(N3,eta)) dN4y=Jinv[2]*(sp.diff(N4,zhi)) + Jinv[3]*(sp.diff(N4,eta)) E=((200*(10**9))/(1-(0.09)))*np.array([[1,0.3,0],[0.3,1,0],[0,0,0.35]]) B=sp.Matrix([[dN1x,0,dN2x,0,dN3x,0,dN4x,0],[0,dN1y,0,dN2y,0,dN3y,0,dN4y],[dN1y,dN1x,dN2y,dN2x,dN3y,dN3x,dN4y,dN4x]]) disp_vec = np.array([[u1],[v1],[u2],[v2],[u3],[v3],[u4],[v4]]).reshape(8,1) #disp_vec = sp.Matrix([u1, v1, u2, v2, u3, v3, u4, v4]) strain1=B*disp_vec gp=sp.Matrix([[-1,-1],[1,-1],[1,1],[-1,1]]) strain=np.zeros([3,4]) for i in range(4): strain2=strain1.subs([(zhi,gp[i,0]),(eta,gp[i,1])]) strain3=np.array(strain2).flatten() strain[:,i]=strain3 strain_x=np.sum(strain[0,:])/4 strain_y=np.sum(strain[1,:])/4 strain_xy=np.sum(strain[2,:])/4 stress=np.zeros([3,4]) for i in range(4): strain_c=strain[:,i].reshape(3,1) stress1=np.dot(E,strain_c) stress2=stress1.flatten() stress[:,i]=stress2 stress_x=np.sum(stress[0,:])/4 stress_y=np.sum(stress[1,:])/4 stress_xy=np.sum(stress[2,:])/4 return strain_x, strain_y,strain_xy,stress_x,stress_y,stress_xy Strain_XX=[] Strain_YY=[] Strain_XY=[] Stress_XX=[] Stress_YY=[] Stress_XY=[] for i in range(len(EL)-1): S1,S2,S3,S4,S5,S6=strs_strn(EL[i]) Strain_XX.append(S1) Strain_YY.append(S2) Strain_XY.append(S3) Stress_XX.append(S4) Stress_YY.append(S5) Stress_XY.append(S6) print(Strain_XX) print(Strain_YY) print(Strain_XY) print(Stress_XX) print(Stress_YY) print(Stress_XY) import matplotlib.pyplot as plt from matplotlib.patches import Polygon from matplotlib.collections import PatchCollection fig, ax = plt.subplots() patches = [] stress_values = [] # Iterate over valid indices in EL2 for i in range(len(EL2)): element = EL2[i] coords = deformed1[element] polygon = Polygon(coords, closed=True) patches.append(polygon) p = PatchCollection(patches, cmap='jet', edgecolors='black', linewidths=1) p.set_array(Stress_XX) ax.add_collection(p) # Set colorbar cb = fig.colorbar(p, ax=ax) cb.set_label('Stress_XX') # Set plot limits ax.autoscale() plt.axis('equal') plt.xlabel('X') plt.ylabel('Y') plt.title('Stress_XX') plt.show() import matplotlib.pyplot as plt from matplotlib.patches import Polygon from matplotlib.collections import PatchCollection # Make sure `deformed1` is properly defined and contains the necessary data # deformed1 = ... fig, ax = plt.subplots() patches = [] stress_values = [] for i in range(len(deformed1)-1): # Make sure `element` is within the bounds of `deformed1` # Make sure `deformed1` is properly defined and contains the necessary data coords = deformed1[EL2[i]] polygon = Polygon(coords, closed=True) patches.append(polygon) print(i) p = PatchCollection(patches, cmap='jet', edgecolors='black', linewidths=1) p.set_array(Stress_XX) ax.add_collection(p) # Set colorbar cb = fig.colorbar(p, ax=ax) cb.set_label('Stress_XX') # Set plot limits ax.autoscale() plt.axis('equal') plt.xlabel('X') plt.ylabel('Y') plt.title('Stress_XX') plt.show() import matplotlib.pyplot as plt from matplotlib.patches import Polygon from matplotlib.collections import PatchCollection fig, ax = plt.subplots() patches = [] for element in EL: if isinstance(element, int): # Check if element is an integer coords = deformed1[element] polygon = Polygon(coords, closed=True) patches.append(polygon) p = PatchCollection(patches, cmap='jet', edgecolors='black', linewidths=1) p.set_array(Stress_XX) ax.add_collection(p) # Set colorbar cb = fig.colorbar(p, ax=ax) cb.set_label('Stress_XX') # Set plot limits ax.autoscale() plt.axis('equal') plt.xlabel('X') plt.ylabel('Y') plt.title('Stress_XX') plt.show() Use the code to rewrite the code full for a rectangular plate with length L and width 'a' with hole in circle with radius r, the relation is a/r>=5, I am not having the csv files, create own nodes and element connection, the craeting mesh should be right half of the plate, plot the stress etc. like in code.