晶体塑性有限元弥勒指数转欧拉角

Question

最近刚入门晶体塑性有限元，用的基于黄永刚的UMAT改的VUMAT子程序，得到了一些6列SDV数据，SDV37-39和SDV73-75，我手动把这些数据提取出来了,我知道这些是变形后的弥勒指数，现在想把这些数据转换成欧拉角来画极图，有什么脚本或者代码来实现吗，有偿

Answer 1

引自免费微信小程序：皆我百晓生

为了将弥勒指数（Muller indices）转换为欧拉角，我们可以采用以下步骤：

1. 将SDV37-39和SDV73-75转换为旋转矩阵。
2. 使用旋转向量表示法从旋转矩阵提取欧拉角。

首先，我们需要导入必要的Python库，然后定义一个函数来执行这个过程。

import numpy as np

def muller_to_euler(sdv37_39, sdv73_75):
    # Step 1: Convert SDV to rotation matrix
    R = np.zeros((3, 3))
    
    R[0, 0] = 1 - 2 * (sdv37_39 ** 2 + sdv73_75 ** 2)
    R[0, 1] = 2 * (sdv37_39 * sdv38 + sdv73_75 * sdv74)
    R[0, 2] = 2 * (sdv37_39 * sdv39 - sdv73_75 * sdv73)

    R[1, 0] = 2 * (sdv37_39 * sdv38 - sdv73_75 * sdv74)
    R[1, 1] = 1 - 2 * (sdv38 ** 2 + sdv74 ** 2)
    R[1, 2] = 2 * (sdv38 * sdv39 + sdv74 * sdv75)

    R[2, 0] = 2 * (sdv37_39 * sdv39 + sdv73_75 * sdv73)
    R[2, 1] = 2 * (sdv38 * sdv39 - sdv74 * sdv75)
    R[2, 2] = 1 - 2 * (sdv39 ** 2 + sdv73 ** 2)

    # Step 2: Extract Euler angles using the ZYX convention
    # Note: Depending on your specific application, you may need to use a different convention
    theta_z = np.arctan2(2 * (R[0][2] * R[1][1] - R[0][1] * R[1][2]), R[0][0] - R[1][1])
    theta_y = np.arcsin(2 * (R[0][2] * R[1][0] + R[0][0] * R[1][2]))
    theta_x = np.arctan2(2 * (R[2][1] * R[0][2] - R[2][2] * R[0][1]), R[2][2] * R[1][1] - R[2][1] * R[1][2])

    return theta_z, theta_y, theta_x

# Assuming data is in a pandas DataFrame called df with columns 'SDV37', 'SDV38', 'SDV39', 'SDV73', 'SDV74', 'SDV75'
df['theta_z'], df['theta_y'], df['theta_x'] = zip(*df.apply(muller_to_euler, axis=1, result_type='expand'))

# To plot the results in a polar plot, you can use matplotlib:
import matplotlib.pyplot as plt

fig, ax = plt.subplots(subplot_kw={'projection': 'polar'})
ax.plot(df['theta_y'], df['theta_x'])
ax.set_rmax(np.max(np.abs(df['theta_x'])) * 1.1)  # Adjust the radius limit if needed
ax.set_rticks([0, np.pi / 2, np.pi, 3 * np.pi / 2])  # Set radial tick marks
ax.set_rlabel_position(-22.5)  # Move radial labels away from plotted line
ax.grid(True)
ax.set_title('Polar Plot of Euler Angles')
plt.show()

这段代码首先定义了一个名为muller_to_euler的函数，用于将SDV值转换为旋转矩阵，并进一步计算出ZYX欧拉角。接着，我们假设数据存储在名为df的pandas DataFrame中，每一行代表一组SDV值。我们应用这个函数到DataFrame上，得到对应的θz、θy和θx，并将其添加为新的DataFrame列。

最后，我们使用matplotlib创建一个极坐标图，展示θy对θx的关系。根据你的实际需求，可能需要调整极坐标图的设置，例如轴限、标签等。

请注意，这里的代码假设你的数据已经按照上述格式组织好并存在DataFrame df 中。如果你的数据结构不同，你需要适当地修改代码以适应你的数据。另外，对于大角度的欧拉角，可能会遇到角度折叠问题，这取决于你的旋转顺序和选择的欧拉角表示方法。如果你在实际应用中有特定的需求或疑问，请随时提问。