2016-01-20 38 views
1

我使用streamplot为了绘制在一个开放的圆圈应力轨迹。我不希望在圆的半径内分析应力轨迹,原因有两个:(1)应力不会通过空气传播,因为它们会穿过孔周围的介质;(2)数学不会允许它。我一直在搞一个面具的想法,但我一直没有能够得到它的工作。可能有更好的方法。有谁知道我如何绘制这些轨迹,而不需要在孔的半径内绘图?我有效地需要某种命令来告诉流图在到达洞的外半径时停止,但是也知道在哪里再次拾起。下面的第一个代码只是用来推导压力轨迹方向的数学运算。我把这个作为参考。在此之后,我绘制了轨迹。Streamplot掩盖一个圆形区域

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) 
+0

难道足以粘贴在streamplot顶部的白色(不透明)圈? – unutbu

+0

不幸的是没有。起初我是这么做的,所以我没有意识到搞乱了什么。当r变为零时,我在数学中的许多地方除以零(注意:c = a/r几次出现)。这是导致错误的绘图的孔,这是检查@unutbu – Ned

+0

什么是理想的(正确的值)'stress_trajectory_cartesian'应该返回作为'r'变为零的最重要的地方附近?或者我应该说'r unutbu

回答

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() 

随着剧情

enter image description here

+0

这似乎是一个肯定的想法!出于某种原因,当我适应这一点代码来适应我的问题时,我没有得到期望的效果。我不是只将nan值应用于半径<1范围内的区域,而是将-1赋予-1和1之间的间隔。我不明白,也不知道为什么它会这样做。我在@nicoguaro – Ned

+0

@Ned上面引入了我的更新代码,如果您更改(或删除)“minlength”,那么您可能会获得更接近我的图像的东西。 – nicoguaro

+0

谢谢! @nicoguaro – Ned

相关问题