管节点相贯线的空间坐标计算（附python源代码） |
您所在的位置：网站首页 › 圆柱体坐标方程 › 管节点相贯线的空间坐标计算（附python源代码） |
管节点相贯线的空间坐标计算（附python源代码）
|
问题背景：在工程结构设计中， 对管节点焊接结构进行疲劳评估， 需要通过参考点线性递推得到焊趾处热点应力。参考论文《管节点热点应力参数化分析方法》中的计算方法，本文对管节点相贯线的曲线建立数学方程， 并通过数值计算方法， 对距离相贯线任意距离的椭圆环线进行离散， 然后通过关键点建模， 对网格生成的大小以及位置进行精确控制。 理论知识： 1、管节点相贯线方程（正交&斜交），本文不考虑偏置。
2、线贯线外延线公式推导
3、第一曲线积分公式
4、高斯-勒让德公式
5、python源代码（本文仅生成相贯线及指定外延距离处的离散点） import numpy as np import matplotlib.pyplot as plt from scipy.optimize import fsolve from sympy import * import math import sys from mpl_toolkits.mplot3d import Axes3D # 圆柱1半径：radius_1; （支管） # 圆柱2半径：radius_2; （主管） # 点距： theta；（按弧度） # 支管轴线与主管轴线的夹角phi（0-pi）(弧度） # 外延距离: d class IntersectingLine(): def __init__(self): self.radius_1 = 8.0 self.radius_2 = 30.0 self.phi = 0.78539 self.step_theta = 2*pi/60 self.step_d=10 # 计算空间坐标，保存在 ndarray中 def calculate_coordinates(self, r1, r2, phi, theta,d): shape = ((int)(pi * 2 / theta), 3) coords = np.zeros(shape, dtype=float) coords_wy=np.zeros(shape, dtype=float) #定义相贯线法向量,取定theta_0后，相贯线处某一点的总体坐标为（x0,y0,z0),此处的法向量为 #（u_0,v_0,w_0) def u_0(theta_0): return -r1*sin(theta_0) def v_0(theta_0): a1 = r1* sin(phi) * cos(theta_0) b1 = r1*r1*r1* cos(phi) * sin(theta_0) * cos(theta_0) c1 = sqrt(r2 ** 2 - (r1 * cos(theta_0))**2) d1 = r1 * cos(phi)**2 * cos(theta_0) return a1 + (b1 / c1 + d1) / sin(phi) def w_0(theta_0): e1 = r1**3 * sin(theta_0) * cos(theta_0) f1 = sqrt(r2 ** 2 - (r1 * cos(theta_0))**2) return e1 / f1 for i in range(shape[0]): coords[i][0] = r1 * cos(i * theta) a1=sqrt(r2*r2-(r1*cos(i*theta))**2)*cos(phi) b1=r1*(cos(phi))**2*sin(i*theta) c1=r1 * sin(i * theta)*sin(phi) coords[i][1] = c1+(a1+b1)/sin(phi) coords[i][2] = sqrt(r2 * r2 - (coords[i][0])**2) # 法平面与主管外表面的交线l0的方程（椭圆曲线） #高斯积分，求上限 def f(x): y_x=-u_0(i*theta)/v_0(i*theta)+w_0(i*theta)*x/(v_0(i*theta)*\ sqrt(r2-x)*sqrt(r2+x)) z_x=-x/(sqrt(r2-x)*sqrt(r2+x)) return sqrt(1+y_x**2+z_x**2) #积分上下限，a是上限，b是下限 a=r1*cos(i*theta) b=symbols('b') # 定常积分采用高斯-勒让德求积公式，已知外延距离d，反算积分上限b I4f=0.5*(b-a)*(f((b-a)*0.5*0.9061798459+(b+a)*0.5)*0.2369268851+\ 0.2369268851*f((b-a)*0.5*(-0.9061798459)+(b+a)*0.5)+\ 0.4786286705*f((b-a)*0.5*(0.5384693101)+(b+a)*0.5)+\ 0.4786286705* f((b - a) * 0.5 * (-0.5384693101) + (b + a) * 0.5)+\ 0.5688888889*f((b+a)*0.5)) if (0 |
【本文地址】
今日新闻 |
推荐新闻 |