[Feature Request] Add option to calculate and visualize Von Mises plane stress for shell elements

[GJoe2]

Hello @yexiang92, I hope you're doing well.

I noticed that plane and solid elements in opstool have the option to visualize Von Mises stress, but this feature does not appear to be available for shell elements. While I understand that Von Mises stress is traditionally associated with quad/plane elements, I believe it can still be very useful for shell elements to assess how the principal in-plane forces interact.

Here’s the standard 2D Von Mises formula that could be applied to membrane forces in shell elements:

$\sigma_{\text{VM}} = \sqrt{ \sigma_{11}^2 - \sigma_{11} \sigma_{22} + \sigma_{22}^2 + 3\sigma_{12}^2 }$

$F_{\text{VM}} = \sqrt{ F_{11}^2 - F_{11} F_{22} + F_{22}^2 + 3F_{12}^2 }$

Thanks again for this great tool, and have a great day!

opstool/opstool/post/_get_response/_get_solid_resp.py

Lines 348 to 392 in 8adef7b

	def _calculate_stresses_measures_4D(stress_array, dtype):
	"""
	Calculate various stresses from the stress values at Gaussian points.
	Parameters:
	stress_array (numpy.ndarray): A 4D array with shape (num_elements, num_gauss_points, num_stresses).
	Returns:
	dict: A dictionary containing the calculated stresses for each element.
	"""
	# Extract individual stress components
	sig11 = stress_array[..., 0] # Normal stress in x-direction
	sig22 = stress_array[..., 1] # Normal stress in y-direction
	sig33 = stress_array[..., 2] # Normal stress in z-direction
	sig12 = stress_array[..., 3] # Shear stress in xy-plane
	sig23 = stress_array[..., 4] # Shear stress in yz-plane
	sig13 = stress_array[..., 5] # Shear stress in xz-plane
	# Calculate von Mises stress for each Gauss point
	sig_vm = np.sqrt(
	0.5 * ((sig11 - sig22) ** 2 + (sig22 - sig33) ** 2 + (sig33 - sig11) ** 2)
	+ 3 * (sig12**2 + sig13**2 + sig23**2)
	)
	# Calculate principal stresses
	# Using the stress tensor to calculate eigenvalues
	stress_tensor = _assemble_stress_tensor_4D(stress_array)
	# Calculate principal stresses (eigenvalues)
	principal_stresses = np.linalg.eigvalsh(stress_tensor) # Returns sorted eigenvalues
	p1 = principal_stresses[..., 2] # Maximum principal stress
	p2 = principal_stresses[..., 1] # Intermediate principal stress
	p3 = principal_stresses[..., 0] # Minimum principal stress
	# Calculate maximum shear stress
	tau_max = np.abs(p1 - p3) / 2
	# Calculate octahedral normal stress
	sig_oct = (p1 + p2 + p3) / 3
	# Calculate octahedral shear stress
	tau_oct = 1 / 3 * np.sqrt((p1 - p2) ** 2 + (p2 - p3) ** 2 + (p3 - p1) ** 2)
	data = np.stack([p1, p2, p3, sig_vm, tau_max, sig_oct, tau_oct], axis=-1)
	return data.astype(dtype["float"])

opstool/opstool/post/_get_response/_get_plane_resp.py

Lines 381 to 408 in 8adef7b

	def _calculate_stresses_measures(stress_array, dtype):
	"""
	Calculate various stresses from the stress values at Gaussian points.
	Parameters:
	stress_array (numpy.ndarray): A 4D array with shape (num_elements, num_gauss_points, num_stresses).
	Returns:
	dict: A dictionary containing the calculated stresses for each element.
	"""
	# Extract individual stress components
	sig11 = stress_array[..., 0] # Normal stress in x-direction
	sig22 = stress_array[..., 1] # Normal stress in y-direction
	sig12 = stress_array[..., 2] # Normal stress in z-direction
	# Calculate von Mises stress for each Gauss point
	sig_vm = np.sqrt(sig11**2 - sig11 * sig22 + sig22**2 + 3 * sig12**2)
	# Calculate principal stresses (eigenvalues)
	p1 = (sig11 + sig22) / 2 + np.sqrt(((sig11 - sig22) / 2) ** 2 + sig12**2)
	p2 = (sig11 + sig22) / 2 - np.sqrt(((sig11 - sig22) / 2) ** 2 + sig12**2)
	# Calculate maximum shear stress
	tau_max = np.sqrt(((sig11 - sig22) / 2) ** 2 + sig12**2)
	data = np.stack([p1, p2, sig_vm, tau_max], axis=-1)
	return data.astype(dtype["float"])