基于C语言的材料力学模型计算实现(扭转变形篇)

本章节分为传动轮轴模型与无传动轮轴模型，两种模型的程序均可计算扭转角、切应力并检验轴是否满足条件，而传动轮轴模型会提供多一个功能：根据轮的转速与功率计算扭矩。

扭转变形程序算法流程图如下：

Part 1.传动轮轴模型

例题1.如图2-1所示的机床传动轴有三个齿轮，主动轮Ⅱ输入功率为1.52kw，从动轮Ⅰ、Ⅲ输出功率均为0.76kw，轴转速n=184r/min，轴的直径D=50mm，轴切变模量G=80GPa.

2-1 例题1图

①.[φ’]=1°/m，请检验该轴是否满足条件.                                      

②.许用切应力=2×10^6Pa，请检验轴是否满足条件.                                              

解题思想与算法：

A.头文件与宏定义、全局变量声明部分：          

#include<stdio.h>
#include<math.h>
#include<stdbool.h>//布尔型态值,true=1或非0整数,false=0 
#define N 5// 
#define Pi acos(-1.0)//圆周率 
double M[N];//外力偶矩 
double T[N];//叠加扭矩 
double P[N];//齿轮功率 
double l[N];//各段长度 
double Angle[N];//扭转角 
double F[N];//切应力

B.主函数与外部函数部分：

(1).按照提示输入对应数据.

运行界面：

代码片段：

int main(){
	bool Choice,Answer;
	int Number,type;
	double n,G,Angle_MAX,F_MAX;
	char Letter[N+1]={'O','A','B','C','D','E'};
	printf("是否有齿轮?是,请输入1;否,请输入0:");
	scanf("%d",&Choice);
if(Choice==true){
	printf("请输入齿轮数目(≤5):");
	scanf("%d",&Number);
	printf("请输入转速(r/min):");
	scanf("%lf",&n); 
	printf("是否输入功率?是,请输入1;否,请输入0:"); 
	scanf("%d",&Choice);

(2). 重新标记各位置字母，以左边第一个齿轮为原点O，往右数的齿轮依次标为A与B，由于要兼容多种计算，程序中的分段数默认为齿轮数，即除了OA、AB段之外，还要增加一段BC，BC的距离长短可任意输入(一般输入0)，不影响结果.(标记方式图如2-2所示)

2-2 标记方式图

程序会根据用户输入的功率与转速调用TorsionCouple(已在代码段标红) 函数得到O、A、B三点扭矩，并将各点外力偶矩存入数组M中，将各点扭矩累加，并分别存入扭矩数组T中；

运行界面：

代码片段：

if(Choice==true){
	printf("请从起始端至末端依次输入齿轮功率(以从右往左看,顺时针记功率为正)、各段距离/m\n");
	for(int i=0;i<Number;i++){
		printf("齿轮%d功率/kw:",i+1);
		scanf("%lf",&P[i]);
		printf("%c%c距离/m:",Letter[i],Letter[i+1]);
		scanf("%lf",&l[i]);
		M[i]=TorsionCouple(P[i],n);
	}
	printf("各段轴的外力偶矩:\n");
	for(int i=0;i<Number;i++){
	printf("%c%c段外力偶矩:%lf\n",Letter[i],Letter[i+1],M[i]);
	T[i]=M[i];
	if(i>0)T[i]=T[i]+T[i-1];//叠加各处外力偶矩
	}
	}

TorsionCouple函数：

double TorsionCouple(double p,double n){//外力偶矩 
	double m;
	m=(9549*p)/n;// 公式M=9549*P/n 
	return m;
}

(3).根据用户需要选择是否计算扭转角，如果需要计算扭转角，则程序执行如下操作：调用W_Ip函数，给W_Ip函数的形参choice传值为true，根据输入的轴的截面形状、参数，计算极惯性矩Ip值。接着，基于用户输入的轴切变模量与(2)中输入的各段轴长短(相邻两点距离)以及(2)中已求出的各处外力偶矩M(相邻两点外力偶矩求差值，得到∆M，即下式中的T)，结合计算得到的Ip值，套用公式φ’=T/(GIp)计算各段轴单位长度扭转角，如果用户选择了输出扭转角的功能，则通过公式∆φ= φ'l计算出各段轴扭转角。最后，用户输入许用扭转角大小，程序对已知的各段扭转角取绝对值进行比较，若全部单位长度扭转角均小于许用单位长度扭转角大小，则输出“刚度满足条件”；否则输出“刚度不满足条件”；

运行界面：(BC段在此题中是多余段，其计算结果无实际意义)

代码片段：

        printf("是否计算刚度?是,请输入1;否,请输入0:");
	scanf("%d",&Choice);
	if(Choice==true){
	printf("请输入轴切变模量G/GPa:");
	scanf("%lf",&G);
	printf("是否计算扭转角?是,请输入1;否,请输入0:");
	scanf("%d",&Choice);
    for(int i=0;i<Number;i++){
		printf("%c%c段截面形状:a.实心圆,输入0;b.空心圆,输入1;c.薄壁筒,输入2:",Letter[i],Letter[i+1]);
		scanf("%d",&type);  
	 if(i>0)AngleDerivative[i]=T[i]*180/(G*pow(10.0,9)*W_Ip(type,true)*Pi); 
else AngleDerivative[i]=M[i]*180/(G*pow(10.0,9)*W_Ip(type,true)*Pi);
		printf("%c%c段单位长度扭转角°/m:%lf\n",Letter[i],Letter[i+1],AngleDerivative[i]);
		if(Choice==true){
		Angle[i]=AngleDerivative[i]*l[i];
		printf("%c%c段扭转角(°):%lf\n",Letter[i],Letter[i+1],Angle[i]); 
		}
	}
	printf("请输入许用单位长度扭转角°/m:");
	scanf("%lf",&Angle_MAX);
	Answer=true;//初始化 
	for(int i=0;i<Number;i++){
		if(Angle_MAX<fabs(AngleDerivative[i]))Answer=false;
	}
	if(Answer==true)printf("刚度满足条件\n");
	else printf("刚度不满足条件\n");
	}

W_Ip函数：形参type传值为用户选择的截面形状代号，当type=0时，函数计算实心圆截面；当type=1时，函数计算空心圆截面；当type=2时，函数计算薄壁筒截面.当形参choice的布尔值为真时，函数计算Ip值；当choice的布尔值为假或者0时，函数计算Wp值。

double W_Ip(int type,bool choice){//扭转截面系数/m^3 
	double Wp,Ip,D,d,R,t;//D为大径或直径,d为小径,R为薄壁筒半径,t为厚度 
	if(!type){
	printf("请输入实心圆直径D/m:");
	scanf("%lf",&D);
	Ip=(Pi*pow(D,4))/32;
	Wp=(Pi*pow(D,3))/16;} 
	else if(type==1){
	printf("请输入空心圆大径D/m:");
	scanf("%lf",&D);
	printf("请输入空心圆小径d/m:");
	scanf("%lf",&d);
	Ip=(Pi*(pow(D,4)-pow(d,4)))/32; 
	Wp=(Pi*pow(D,3))*(1-pow((d/D),4))/16; 
	}
	else{
	printf("请输入薄壁筒直半径R/m:");
	scanf("%lf",&D);
	printf("请输入薄壁筒厚度t/m:");
	scanf("%lf",&t);
	Wp=2*Pi*pow(R,2)*t;
	}
	if(choice==true)return Ip;
	else return Wp; 
}

(4).根据用户需要选择是否计算切应力，如果需要计算切应力，则程序执行如下操作：仍调用(3)中的W_Ip函数，但给W_Ip函数的形参choice传值为false，根据输入的轴的截面形状、参数，计算扭转截面系数Wp值。接着，基于已在(2)中求出的各段扭矩T(各段M的总叠加值)，根据公式=T/Wp得到各段切应力。最后，用户输入许用切应力大小，程序对已知的各段切应力取绝对值进行比较，若全部切应力均小于许用切应力大小，则输出“切应力满足条件”；否则输出“切应力不满足条件”，程序结束。

运行界面：(BC段在此题中是多余段，其计算结果无实际意义)

代码片段：

        printf("是否计算切应力?是,请输入1;否,请输入0:");
	scanf("%d",&Choice);
	if(Choice==true){
	for(int i=0;i<Number;i++){
		printf("%c%c段截面形状:a.实心圆,输入0;b.空心圆,输入1;c.薄壁筒,输入2:",Letter[i],Letter[i+1]);
		scanf("%d",&type);
		F[i]=T[i]/W_Ip(type,false); 
	    printf("%c%c段切应力/Pa:%lf\n",Letter[i],Letter[i+1],F[i]);
	}
	printf("请输入许用切应力/Pa:");
	scanf("%lf",&F_MAX);
	Answer=true;
	for(int i=0;i<Number;i++){
		if(F_MAX<fabs(F[i]))Answer=false;
	}
	if(Answer==true)printf("切应力满足条件\n");
	else printf("切应力不满足条件\n");
	} 
	return 0;
}

Part 2.无传动轮轴模型

无传动轮轴模型的大致内容与传动轮轴模型相同，区别在于没有齿轮，不通过输入功率、转速的参数求解外力偶矩大小，而是直接输入各点处的外力偶矩，其余流程均一致。

运行界面：

代码片段：(else语句内为无传动轮轴模型的代码内容)

int main(){
	bool Choice,Answer;
	int Number,type;
	double n,G,Angle_MAX,F_MAX;
	char Letter[N+1]={'O','A','B','C','D','E'};
	printf("是否有齿轮?是,请输入1;否,请输入0:");
	scanf("%d",&Choice);
	if(Choice==true){
	//2.1中齿轮参数输入与函数调用、计算的内容
	}
	//非传动轮轴模型特有的代码片段：
	else{printf("请输入外力偶矩数目(≤5):");
	scanf("%d",&Number);
	printf("请从左至右输入各外力偶矩:\n");
	for(int i=0;i<Number;i++){
		printf("位置%d外力偶矩/Nm:",i+1);
		scanf("%lf",&M[i]); 
		printf("%c%c距离/m:",Letter[i],Letter[i+1]);
	    scanf("%lf",&l[i]);
	    T[i]=M[i];
		if(i>0)T[i]=T[i]+T[i-1];
	}
	}

Part 3.扭转变形程序的应用

变式1.某小型水电站水轮机容量为50kw，转速为300r/min，钢轴直径为75mm。若在正常运转下且只考虑扭矩作用，其许用切应力=20MPa。试校核轴的强度。

解题思路：此题属于传动轮模型，只计算1个齿轮，运行界面会默认要求输入对应段轴长，而轴强度计算无需轴长距离的参数，可输入任意数值(一般输入1)。其他参数根据界面提示输入(注意:直径75mm要按照界面要求换成m的单位，即0.075m)，即可得到该轴最大切应力与强度是否符合许用标准的结果。

结论：该轴最大切应力约为19.21MPa，满足强度条件.

变式2.如左图所示阶梯型圆杆，AE段为空心，外直径D=140mm，内直径d=100mm；BC段为实心，直径d=100mm。外力偶矩MA=18kN·m，MB=32kN·m，MC=14kN·m。已知：=80MPa，[φ']=1.2°/m，G=80GPa。试校核该轴的强度和刚度。

解题思路：此题属于无传动轮轴模型，由于程序默认的点标号从起始段开始数分别为O、A、B、C…因此，将A、E、B、C四点分别视为O、A、B、C(如右图所示)，MA即为Mo。由于各点均要求输入外力偶矩，而点E(对应界面的A)无外力偶矩，因此直接输入0即可。题目并未给出AE、EB、BC实际距离，可输入任意数值(本题输入1)，不会对结果产生影响；而对于多出来CD段(对应题目C点往右任意距离的虚拟轴段)，输入距离参数(本题输入0)与输出结果无实际意义，通过输入界面要求的参数，即可求解。

运行界面：

结论：该轴最大单位长度扭转角约为1.02°，满足刚度条件；该轴最大切应力约为71.30MPa，满足强度条件.

变式3.已知钻探机钻杆(如左图所示)的外径D=60mm，内径d=50mm，功率P=7.355kW，转速180r/min，钻探入土深l=40m，钻探材料的G=80GPa，许用切应力=40Mpa，[φ']=1.0°/m。假设土壤对钻杆的阻力是沿长度均匀分布的，试校核其强度并求出两端截面的相对扭转角。

解题思路：如图2-5所示，以刚进入土层的截面O为原点，钻杆底端为A点，将钻杆逆时针旋转90°横置，作扭矩图得：

由扭矩图与相对扭转角积分结果可知，本题的相对扭转角不是∆φ=，而是∆φ=。由于危险截面O上的外力偶矩Me与GIp是通过用户输入参数、计算机计算得到的，不方便进行修改。因此要想得到本题的相对扭转角，直接将OA之间的距离输入值改为，即输入20，即可得到两端截面的相对扭转角，不会影响其他计算结果。正常输入界面要求的其他参数，即可借助Me计算出钻杆内最大切应力，并进行强度校核。

结论：钻杆内最大切应力约为17.78MPa，满足强度条件；两端截面的相对扭转角约为8.48°.

       扭转变形部分的内容到此告一段落了，如有错漏，欢迎指正本渣(本人物理真的菜QAQ)。下一篇博客研究的是拉压变形，欢迎继续阅读~