问题背景：在工程结构设计中， 对管节点焊接结构进行疲劳评估， 需要通过参考点线性递推得到焊趾处热点应力。参考论文《管节点热点应力参数化分析方法》中的计算方法，本文对管节点相贯线的曲线建立数学方程， 并通过数值计算方法， 对距离相贯线任意距离的椭圆环线进行离散， 然后通过关键点建模， 对网格生成的大小以及位置进行精确控制。

理论知识：

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<=i*theta<0.5*pi or 1.5*pi<i*theta<2*pi):
                res=nsolve(I4f-d,0)
                coords_wy[i][0] = res
                coords_wy[i][2] = sqrt(r2 ** 2 - res ** 2)
                coords_wy[i][1] = (-u_0(i * theta) * (res - coords[i][0]) - \
w_0(i * theta) * (coords_wy[i][2] - coords[i][2])) / v_0(i * theta) + coords[i][1]
            else:
                if (0.5*pi<i*theta<1.5*pi):
                    res = nsolve(I4f+d, 0)
                    coords_wy[i][0] = res
                    coords_wy[i][2] = sqrt(r2 ** 2 - res ** 2)
                    coords_wy[i][1] = (-u_0(i * theta) * (res - coords[i][0]) - \
w_0(i * theta) * (coords_wy[i][2] - coords[i][2])) / v_0(i * theta) + coords[i][1]
                else:
                    if math.isclose(i*theta,0.5*pi):
                       coords_wy[i][0]=0
                       coords_wy[i][1]=coords[i][1]+d
                       coords_wy[i][2]=r2
                    else:
                        coords_wy[i][0] = 0
                        coords_wy[i][1] = coords[i][1]-d
                        coords_wy[i][2] = r2
            # 求解外延线坐标,将res代入法平面与主管外表面交线l0的方程，得出外延点坐标
            #coords_wy[i][0]=res
            #coords_wy[i][2] = sqrt(r2**2-res**2)
            #coords_wy[i][1]=(-u_0(i*theta)*(res-coords[i][0])-w_0(i*theta)*\(coords_wy[i][2]-coords[i][2]))/v_0(i*theta)+coords[i][1]
        return coords,coords_wy
    # 根据坐标点绘制三维散点图
    def plot_coordinates(self, coordinates):
        ax = plt.subplot(projection='3d')
        ax.set_title('3d_image_show')  # 设置本图名称
        ax.scatter3D(coordinates[:, 0], coordinates[:, 1], coordinates[:, 2], c='r')
        ax.set_xlabel('X')  # 设置x坐标轴
        ax.set_ylabel('Y')  # 设置y坐标轴
        ax.set_zlabel('Z')  # 设置z坐标轴
        plt.show()
joint1=IntersectingLine()
joint1.radius_1=200
joint1.radius_2=500
joint1.phi=(1/6)*pi
joint1.step_theta=2*pi/60
joint1.step_d=50
coords1,coords2 = joint1.calculate_coordinates(joint1.radius_1, joint1.radius_2, joint1.phi, joint1.step_theta, joint1.step_d)
coords_x=[str(x[0]) for x in coords1]
coords_y=[str(y[1]) for y in coords1]
coords_z=[str(z[2]) for z in coords1]
for i in range(0,len(coords_x)):
   apdl='K'+','+str(i+1)+','+coords_x[i]+','+coords_y[i]+','+coords_z[i]
   with open('data.txt','a',newline="") as f:
       f.write(apdl)
       f.write('\r\n')
coords2_x=[str(x[0]) for x in coords2]
coords2_y=[str(y[1]) for y in coords2]
coords2_z=[str(z[2]) for z in coords2]
for j in range(0,len(coords2)):
    apdl2 = 'K' + ',' + str(j+i+2) + ',' + coords2_x[j] + ',' + coords2_y[j] + ','\ ++coords2_z[j]
    with open('data.txt', 'a', newline="") as f:
        f.write(apdl2)
        f.write('\r\n')