我正在使用streamplot来绘制围绕空心圆的应力轨迹。我不希望在圆半径内分析应力轨迹，这有两个原因：（1）应力不会像通过孔周围的介质那样通过空气传播，并且（2）数学不允许它。我一直在想着口罩的概念，但是我一直无法使它正常工作。可能有更好的方法。有谁知道我如何绘制这些轨迹而无需在孔的半径内绘制呢？我实际上需要某种命令来告诉流图在到达孔的外部半径时停止，但同时又知道在哪里重新取回。下面的代码的第一位只是用于导出应力轨迹方向的数学运算。我将其作为参考。在此之后，我绘制了轨迹。
import numpy as np
import matplotlib.pyplot as plt
from pylab import *
def stress_trajectory_cartesian(X,Y,chi,F,a):
# r is the radius out from the center of the hole at which we want to know the stress
# Theta is the angle from reference at which we want to know the stress
# a is the radius of the hole
r = np.sqrt(np.power(X,2)+np.power(Y,2))*1.0
c = (1.0*a)/(1.0*r)
theta = np.arctan2(Y,X)
A = 0.5*(1 - c**2. + (1 - 4*c**2. + 3*c**4.)*np.cos(2*theta))
B = 0.5*(1 - c**2. - (1 - 4*c**2. + 3*c**4.)*np.cos(2*theta))
C = 0.5*(1 + c**2. - (1 + 3*c**4.)*np.cos(2*theta))
D = 0.5*(1 + c**2. + (1+ 3*c**4.)*np.cos(2*theta))
E = 0.5*((1 + 2*c**2. - 3*c**4.)*np.sin(2*theta))
tau_r = 1.0*F*c**2. + (A-1.0*chi*B) # Radial stress
tau_theta = -1.*F*c**2. + (C - 1.0*chi*D) # Tangential stress
tau_r_theta = (-1 - 1.0*chi)*E # Shear stress
tau_xx = .5*tau_r*(np.cos(2*theta)+1) -1.0*tau_r_theta*np.sin(2*theta) + .5*(1-np.cos(2*theta))*tau_theta
tau_xy = .5*np.sin(2*theta)*(tau_r - tau_theta) + 1.0*tau_r_theta*np.cos(2*theta)
tau_yy = .5*(1-np.cos(2*theta))*tau_r + 1.0*tau_r_theta*np.sin(2*theta) + .5*(np.cos(2*theta)+1)*tau_theta
tan_2B = (2.*tau_xy)/(1.0*tau_xx - 1.0*tau_yy)
beta1 = .5*np.arctan(tan_2B)
beta2 = .5*np.arctan(tan_2B) + np.pi/2.
return beta1, beta2
# Functions to plot beta as a vector field in the Cartesian plane
def stress_beta1_cartesian(X,Y,chi,F,a):
return stress_trajectory_cartesian(X,Y,chi,F,a)[0]
def stress_beta2_cartesian(X,Y,chi,F,a):
return stress_trajectory_cartesian(X,Y,chi,F,a)[1]
#Used to return the directions of the betas
def to_unit_vector_x(angle):
return np.cos(angle)
def to_unit_vector_y(angle):
return np.sin(angle)
下面的代码绘制了应力轨迹：
# Note that R_min is taken as the radius of the hole here
# Using R_min for a in these functions under the assumption that we don't want to analyze stresses across the hole
def plot_stresses_cartesian(F,chi,R_min):
Y_grid, X_grid = np.mgrid[-5:5:100j, -5:5:100j]
R_grid = np.sqrt(X_grid**2. + Y_grid**2.)
cart_betas1 = stress_beta1_cartesian(X_grid,Y_grid,chi,F,R_min)
beta_X1s = to_unit_vector_x(cart_betas1)
beta_Y1s = to_unit_vector_y(cart_betas1)
beta_X1s[R_grid<1] = np.nan
beta_Y1s[R_grid<1] = np.nan
cart_betas2 = stress_beta2_cartesian(X_grid,Y_grid,chi,F,R_min)
beta_X2s = to_unit_vector_x(cart_betas2)
beta_Y2s = to_unit_vector_y(cart_betas2)
beta_X2s[R_grid<1] = np.nan
beta_Y2s[R_grid<1] = np.nan
fig = plt.figure(figsize=(5,5))
#streamplot
ax=fig.add_subplot(111)
ax.set_title('Stress Trajectories')
plt.streamplot(X_grid, Y_grid, beta_X1s, beta_Y1s, minlength=0.9, arrowstyle='-', density=2.5, color='b')
plt.streamplot(X_grid, Y_grid, beta_X2s, beta_Y2s, minlength=0.9, arrowstyle='-', density=2.5, color='r')
plt.axis("image")
plt.xlabel(r'$\chi = $'+str(round(chi,1)) + ', ' + r'$F = $'+ str(round(F,1)))
plt.ylim(-5,5)
plt.xlim(-5,5)
plt.show()
plot_stresses_cartesian(0,1,1)
最佳答案
我认为您只需要为不想考虑的区域提供NaN值。我在下面生成了一个简单的示例。
import numpy as np
import matplotlib.pyplot as plt
Y, X = np.mgrid[-5:5:100j, -5:5:100j]
R = np.sqrt(X**2 + Y**2)
U = -1 - X**2 + Y
V = 1 + X - Y**2
U[R<1] = np.nan
V[R<1] = np.nan
plt.streamplot(X, Y, U, V, density=2.5, arrowstyle='-')
plt.axis("image")
plt.savefig("stream.png", dpi=300)
plt.show()
有剧情