Диссертация (Разработка метода расчета сопротивления качению и теплообразования в массивных шинах при стационарных режимах движения), страница 26

Please, change the net of the construction");exit(1);}7374757677//**********************************//Переопределение области контакта *//**********************************Contact_TruePenetration (ContactIndex,&NormalFlag);7879} while (NormalFlag == 1);8081printf("****************************");178828384printf("\n* END OF PROBLEM SOLUTION *\n");printf("****************************\n");printf("1.Full number of contact cycles = %d\n", cttstep);858687free(Rcon);free(Kcon);8889909192//*****************************//Вывод результатов на печать *//*****************************PrintResult (Ucon,Qcon,ContactIndex);9394959697free(ContactIndex);free(Ucon);free(Qcon);}Процедура «Contact_InitialPenetration»1234void Contact_InitialPenetration (int* ContactIndex, int* NormalFlag)/* Функция, определяющая принадлежность элементоввозможной области контакта5678910111213141.

ContactIndex.........Массив, определяющий принадлежностьКЭ области контакта*/{intNelem;//Номер элементаdouble x[3][8];//Глобальные координаты узлов КЭinti,j,k;//Переменные циклаdouble nv[3];//Вектор нормали к поверхности барабанаdouble gapn[4];//Вектор внедрения15161718NormalFlag = 0;for (Nelem = 1; Nelem<=Ne; Nelem++){1920212223if (ContactIndex[Nelem-1]==0)continue;Position (Nelem, x);179242526272829303132333435for(i=0;i<Nint;i++)for(j=0;j<Nint;j++){Drum_Norm_And_Gap (Nelem, x, -1., GaussP[i], GaussP[j], nv, gapn);if (gapn[2]>=0){ContactIndex[Nelem-1] = 0;goto next;}}next: continue;}363738394041k=0;for (j=0;j<Ne;j++)k+=ContactIndex[j];printf("Initial number of elements in contact: %d\n", k);CEnumering (ContactIndex);4243}Процедура «Contact_TruePenetration»123void Contact_TruePenetration (int* ContactIndex, int* NormalFlag)/* Функция, определяющая принадлежность элементадействительной области контакта456789101112131415161718191.ContactIndex.........Массив, определяющий принадлежность КЭк области контакта*/{intNelem;//Номер элементаdouble x[3][8];//Глобальные координаты узлов КЭdouble xNew[3][8];//Новые координаты узлов КЭ после деформацииdouble u[24];//Вектор перемещений узлов КЭinti,j,k,num;//Переменные циклаintr,s,q;//Переменные циклаdouble nv[3];//Вектор нормали к поверхности барабанаdouble gapn[4];//Вектор внедренияdouble S[6],Sgp[9][6]; /*Напряжения, вычисляемыев гауссовых точках */double EL[6],SL[6];//История напряжений и деформаций18020212223242526*NormalFlag = 0;r = NeR-1;for (s=0;s<NeZ;s++){zerro1((double*)StressLast,Nint*Nint*Nint*6);zerro1((double*)StrainLast,Nint*Nint*Nint*6);2728293031323334353637383940414243444546474849for (q=0;q<NeF;q++){Nelem = NeR*NeZ*q+NeR*s+r+1;Position (Nelem, x);Displacement (Nelem, u);XNglob(Nelem, x, u, xNew);num = 0;for (k=0;k<6;k++){EL[k] = StrainLast[num][k];SL[k] = StressLast[num][k];}StressCylCS (Nelem, x, u, 0., 0., 0., EL, SL, S);for (k=0;k<6;k++){StrainLast[num][k] = EL[k];StressLast[num][k] = SL[k];}for (k=0;k<6;k++){Sgp[num][k] = S[k];}5051525354num = 0;for(i=0;i<Nint;i++)for(j=0;j<Nint;j++){555657585960Drum_Norm_And_Gap (Nelem, xNew, -1., GaussP[i], GaussP[j], nv, gapn);if((gapn[2]>0.0)&&(ContactIndex[Nelem-1] == 0)){goto next;}18161626364656667686970717273747576777879808182838485868788if (Sgp[num][0]>0.0){if(ContactIndex[Nelem-1] == 0){printf("Element %d wants pull into contact\n",Nelem);goto next;}ContactIndex[Nelem-1] = 0;printf("Exclude element %d\n",Nelem);*NormalFlag = 1;goto next;}}if (ContactIndex[Nelem-1] ==0){ContactIndex[Nelem-1] = 1;*NormalFlag = 1;printf("Element %d is a new contact element\n",Nelem);}next: continue;}}k=0;for (j=0;j<Ne;j++)if (ContactIndex[j]>=1) k++;printf("Number of elements in contact: %d\n", k);CEnumering (ContactIndex); //Локальная нумерация контактных КЭ}Процедура «BaseStiffnes»123void BaseStiffnes (long Nelem, double x[3][8],\double KBelem[24][24], double Felem[24])/* Процедура формирования базовой МЖЭ4567891.Nelem..................Номер конечного элемента2.x[3][8]................Глобальные координаты узлов конечного элемента3.KBelem[24][24].........Базовая матрица жесткости конечного элемента4.Felem..................Вектор внутренних узловых сил*/1011{1821213int i,j,k;double dh[8][3],detJ0;141516double B[6][24];double C[6][6];/*Производные функций формыпо глобальным координатам, Якобиан*///Матрица производных функций форм//Вычисление матрицы касательных модулей упругости1718192021222324252627for (i=0; i<Nint; i++)for (j=0; j<Nint; j++)for (k=0; k<Nint;k++){Jacobi(Nelem,GaussP[k],GaussP[j],GaussP[i],x,dh, &detJ0);LeniarSDmatrix (dh, B);Tangent_Matrix (C);AtBAmul(arrayptr(B),6,24,arrayptr(C),\Wgt[k]*Wgt[j]*Wgt[i]*detJ0, arrayptr(KBelem), 1);}282930if (ViscosityFlag == 1)EVS_Force (Nelem, x, Felem);3132}Процедура «Element_StiffnessAndForce»12345void Element_StiffnessAndForce (long Nelem, double x[3][8],\double Kelem[24][24],\double Felem[24], int* ContactIndex)/*Процедура построения матрицы жесткости ивектора внутренних узловых сил конечного элемента67891011121314151.Nelem..................Номер конечного элемента2.x[3][8]................Глобальные координаты узлов конечного элемента3.Kelem[24][24]..........Матрица жесткости конечного элемента4.Felem..................Вектор внутренних узловых сил элемента5.ContactIndex...........Вектор, содержащий информациюо принадлежности КЭ области контакта*/{int i,j,k;161718//************************************************************//Собираем матрицу жесткости конечного элемента*18319//************************************************************202122zerro1((double*)Kelem,24*24);zerro1(Felem,24);2324252627//************************************// Базовая МЖЭ и вектор узловых сил *//************************************BaseStiffnes (Nelem, x, Kelem, Felem);28293031323334353637//****************************************************// МЖЭ и вектор узловых сил, обусловленные контактом.*// Они сразу добавляются к базовой части*//****************************************************if (ContactIndex[Nelem-1] >= 1){ECN_StiffnesAndForce (Nelem, x, Kelem, Felem, 1);ECT_StiffnesAndForce (Nelem, x, Kelem, Felem, 1);}3839}Процедура «Drum_Norm_And_Gap»123456789101112void Drum_Norm_And_Gap (long Nelem, double x[3][8], double ksi,\double etta, double dzeta, double nv[3],\double gapn[4])/* Определение вектора нормали к поверхности барабана.1.Nelem...............Номер элемента2.x...................Координаты узлов элемента в глобальнойсистеме координат3.ksi,etta, dzeta.....Локальные координаты точки4.nv..................Вектор нормали к барабану5.gapn................Вектор, содержащий сведения о взаимномпроникании контактирующих тел*/13141516{double xg[3];int i,j,k;//Вектор глобальных координат1718XGauss (x, ksi, etta, dzeta, xg);18419202122nv[0] = xg[0]/sqrt(pow(xg[0],2)+pow(Ydrum-xg[1],2));nv[1] = (xg[1]-Ydrum)/sqrt(pow(xg[0],2)+pow(Ydrum-xg[1],2));nv[2] = 0.0;232425262728gapn[0]gapn[1]gapn[2]gapn[3]}====xg[0]-Rdrum*nv[0];xg[1]-Ydrum-Rdrum*nv[1];gapn[0]*nv[0]+gapn[1]*nv[1];sqrt(pow(gapn[0],2)+pow(gapn[1],2));Процедура «ECN_StiffnesAndForce»1234void ECN_StiffnesAndForce (long Nelem, double x[3][8],\double KCelem[24][24], double Fcon[24], int Flag)/* Формирование МЖЭ, находящегося в контакте с барабаном.Определяется лишь часть МЖ, обусловленная нормальным контактом.567891011121.Nelem...............Номер элемента2.x...................Координаты узлов элемента в глобальной системе координат3.KCelem..............Матрица жесткости элемента, обусловленная контактом4.Fcon................Вектор внутренних узловых сил от контакта5.Flag................Флаг учета элементов матрицы жесткости и вектор узловых сил*/{1314double gapn[4];1516171819doubledoubledoubledoublenv[3];N[3][3];H[3][24];dS;inti,j,k,l;2021/*Вектор, содержащий сведения о взаимномпроникании контактирующих тел*///Вектор нормали к барабану//Матрица, составленная из компонент вектора нормали//Матрица функций форм/*Отношение площадей в глобальной и в локальной системахкоординат*/2223242526272829//*************************************if (Flag == 0){zerro1(arrayptr(KCelem),24*24);zerro1(Fcon,24);}185303132//**************************************************//Интегрирование по гауссовым точкам поверхности*//**************************************************3334353637for (k=0; k<Nint; k++)for (l=0; l<Nint; l++){Drum_Norm_And_Gap (Nelem, x, -1., GaussP[k], GaussP[l], nv, gapn);38394041N[0][0] = nv[0]*nv[0]; N[0][1] = nv[0]*nv[1]; N[0][2] = nv[0]*nv[2];N[1][0] = nv[1]*nv[0]; N[1][1] = nv[1]*nv[1]; N[1][2] = nv[1]*nv[2];N[2][0] = nv[2]*nv[0]; N[2][1] = nv[2]*nv[1]; N[2][2] = nv[2]*nv[2];42434445464748495051Aria(Nelem,-1., GaussP[k], GaussP[l], x, &dS);Hmatrix (-1., GaussP[k], GaussP[l], H);AtBAmul(arrayptr(H),3,24,arrayptr(N),\kn*Wgt[k]*Wgt[l]*dS, arrayptr(KCelem), 1);for (j=0; j<23; j++)for (i=0; i<2; i++)Fcon[j]+= -kn*Wgt[k]*Wgt[l]*dS*gapn[2]*H[i][j]*nv[i];}}Процедура «ECT_StiffnesAndForce»1234void ECT_StiffnesAndForce (long Nelem, double x[3][8], double KCTelem[24][24],\double FTcon[24], int Flag)/*Формирование матрицы жесткости и вектора внешней нагрузки,обусловленных наличием сцепления в области контакта567891011121314151.Nelem...............Номер элемента2.x...................Координаты узлов элемента в глобальнойсистеме координат3.KTCelem.............Матрица жесткости элемента, обусловленнаятангенциальным контактом4.FTcon...............Вектор внутренних узловых сил, обусловленных контакта5.Flag................Флаг формирования матрицы жесткостии вектора узловых сил*/{1617inti,j,k,l,m;//Параметр цикла18618double DH[3][24];1920double tz[3]={0,0,1},T[3][3];2122double FC[24];/*Матрица производных функций форм,вычисленных в окружном направлении*//*Касательный вектор в осевом направлениии соответствующая ему матрица*///Вспомогательный вектор2324double Sint;2526272829doubledoubledoubledoubleLint;F[24];gapT0;Kelem[24][24];/*Коэффициент преобразованияплощадей*///Коэффициент преобразования длин//Вспомогательный вектор//Зазор в тангенциальном направлении//Вспомогательная матрица303132if(ContactIndex[Nelem-1]==1) return;333435363738if (Flag == 0){zerro2(arrayptr(KCTelem),24,24);zerro1(FTcon,24);}39404142//************************************************//Интегрирование по гауссовым точкам поверхности *//************************************************43444546474849505152535455565758for(i=0;i<Nint;i++)for(j=0;j<Nint;j++){//**************************// Окружное направление*//**************************StickMatrix (Nelem, x, -1., GaussP[i], GaussP[j], FC, F, &gapT0);Aria(Nelem,-1., GaussP[i], GaussP[j], x, &Sint);for (k=0; k<24; k++)FTcon[k]+=-kt*gapT0*F[k]*Wgt[i]*Wgt[j]*Sint;for (k=0; k<24; k++)for (l=0; l<24; l++)KCTelem[k][l] +=kt*Wgt[i]*Wgt[j]*Sint*F[k]*F[l];//************************//Осевое направление *187596061626364656667//************************CircleDerivative (Nelem, x, -1., GaussP[i], GaussP[j], DH);for (l=0;l<3;l++)for (m=0;m<3;m++)T[l][m] = tz[l]*tz[m];AtBAmul(arrayptr(DH),3,24,arrayptr(T),\kt*pow(OmegaT,2)*Wgt[i]*Wgt[j]*Sint, arrayptr(KCTelem), 1);}}Процедура «StickMatrix»123456789101112131415161718void StickMatrix (long Nelem, double x[3][8], double ksi ,double etta,\double dzeta, double FC[24], double F[24], double* gapT0)/*Вычисление матриц, необходимых для построения матрицы жесткостии вектора узловых сил в задаче сцепления1.Nelem...............Номер элемента2.x...................Координаты узлов элемента в глобальнойсистеме координат3.ksi,etta, dzeta.....Локальные координаты точки4.FC[24],F[24]........Вспомогательные вектора5.gapT0...............Зазор в тангенциальном направлении*/{inti,j,k;//Переменные циклаdouble xg[3];//Координаты текущей гауссовой точки шиныdouble psi;//Угловая координатаdouble tb[3];/*Касательный вектор к поверхности барабанав окружном направлении*/double H[3][24];//Матрица функций форм192021222324252627//Векторы и матрицы, вычисляемые для первого элемента рядаlongNelemLast;double xLast[3][8];double xgLast[3];double psiLast;double tbLast[3];double uLast[24];//Вектор перемещений282930XGauss (x, ksi, etta, dzeta, xg);psi = xg[0]/Ro;1883132333435363738tb[0] = 1.0L; tb[1] = xg[0]/Rdrum; tb[2] = 0.0L;Hmatrix (ksi, etta, dzeta, H);NelemLast = Nelem-(ContactIndex[Nelem-1]-1)*NeR*NeZ;Position (NelemLast, xLast);XGauss (xLast, ksi, etta, dzeta, xgLast);psiLast = xgLast[0]/Ro;tbLast[0] = 1.0L; tbLast[1] = xgLast[0]/Rdrum;tbLast[2] = 0.0L;Displacement (NelemLast, uLast);3940414243444546474849zerro1(F,24);for (k=0;k<24;k++)for (i=0;i<3;i++)F[k]+=tb[i]*H[i][k];*gapT0 = (Ro-OmegaD/OmegaT*Rdrum)*(psi-psiLast);for(i=0;i<3;i++)for(k=0;k<24;k++)*gapT0 +=-H[i][k]*uLast[k]*tbLast[i];return;}Процедура «EVS_Force»123void EVS_Force (long Nelem, double x[3][8],double FV[24])/*Функция, формирующая матрицу жесткостии вектор нагрузки вязкоупругого материала456781.Nelem...............Номер передаваемого элемента2.x...................Глобальные координаты узлов элемента3.FV..................Вектор вязких узловых сил*/91011121314{int i,j,k,l,m, num=0;int nume[3];double dt;double dh[8][3], detJ0;151617181920doubledoubledoubledoubledoubleB[6][24];u[24];E[6];EL[6],SL[6];ELrot[6],SLrot[6];//Переменные цикла//Положение элемента в конструкции//Приращение времени/*Производные функций формпо глобальным координатам и Якобиан*///Матрица производных функций форм//Перемещения узлов элемента//Вектор деформаций//Предшествующие напряжения и деформации/*Предшествующие напряжения и деформации1892122double Sv[6],S[6];(повернутые тензоры)*///Вязкие и полные напряжения в текущем элементе23242526Displacement (Nelem, u);NumDetem (Nelem, nume);dt = DTime[nume[0]];2728293031for (i=0;i<Nint;i++)for (j=0;j<Nint;j++)for (k=0;k<Nint;k++,++num){32333435363738394041424344454647484950515253Jacobi(Nelem,GaussP[i], GaussP[j], GaussP[k], x, dh, &detJ0);LeniarSDmatrix (dh, B);Deformation (B, u, E);for (l=0;l<6;l++){EL[l] = StrainLast[num][l];SL[l] = StressLast[num][l];}ConvertRotSS (Nelem, alpha, EL, ELrot);ConvertRotSS (Nelem, alpha, SL, SLrot);Stress_ViscoElastic(Nelem, ELrot, SLrot, E, dt, Sv, S);for (l=0;l<6;l++){StrainLast[num][l] = E[l];StressLast[num][l] = S[l];}for (l=0; l<24;l++)for (m=0; m<6; m++)FV[l] += -B[m][l]*Sv[m]*Wgt[i]*Wgt[j]*Wgt[k]*detJ0;}}Процедура «KRungeKutt»1234void KRungeKutt(long Nelem, double E[6], double DE[6],\double dt, double Spoint[6], double K[6])/* Функция, вычисляющая вспомогательные векторыметода Рунге_Кутта в момент времени alfa*dt561.Nelem...............Номер элемента1907891011121314151617182.E...................Вектор полных деформаций3.DE..................Вектор производных полных деформаций по времени5.Spoint..............Вектор полных напряжений*/{inti;//Переменная циклаintnum[3];//Место элемента в конструкцииdouble DSe[6],DSv[6];/*Вектор скоростей упругихи составной части вязких напряжений*/double Eaprxn[6];/*Тензор деформаций,посчитанный в промежуточный момент времени*/double C;//Мгновенный модуль упругости19202122232425DVscStr(Nelem, E, Spoint, DSv);C = C2+C4;Stress_Hook (DE,C1,C,DSe);for (i=0;i<6;i++)K[i] = (DSe[i] + DSv[i])*dt;}Процедура «Stress_ViscoElastic»123void Stress_ViscoElastic(long Nelem, double EL[6], double SL[6],\double E[6], double dt, double Sv[6], double S[6])/*Функция, вычисляющая вязкоупругие напряжения в текущем элементе456789101.2.3.4.*/{Nelem........Номер элементаEL,SL........Предшествующие напряжения и деформацииdt...........Шаг интегрированияSv,S.........Вычисляемые вязкие и полные напряжения11121314151617int l;doubledoubledoubledoubledouble//Переменная циклаEcalc[6],DEcalc[6];//Промежуточные деформации и их скоростиSpoint[6];//Промежуточное напряжениеK1[6],K2[6],K3[6],K4[6]; //Промежуточные векторыSA[6];//Напряжение звена AC;//Суммарный модуль181920/*Используется процедура интегрирования Рунге-Кутта 4-го порядка*/19121222324252627282930313233343536373839404142StrainAprxn(Nelem, E, EL, 0.0, dt, Ecalc, DEcalc);KRungeKutt(Nelem, Ecalc, DEcalc, dt, SL, K1);StrainAprxn(Nelem, E, EL, 0.5, dt, Ecalc, DEcalc);for (l=0;l<6;l++)Spoint[l] = SL[l]+K1[l]/2;KRungeKutt(Nelem, Ecalc, DEcalc, dt, Spoint, K2);for (l=0;l<6;l++)Spoint[l] = SL[l]+K2[l]/2;KRungeKutt(Nelem, Ecalc, DEcalc, dt, Spoint, K3);StrainAprxn(Nelem, E, EL, 1.0, dt, Ecalc, DEcalc);for (l=0;l<6;l++)Spoint[l] = SL[l]+K3[l];KRungeKutt(Nelem, Ecalc, DEcalc, dt, Spoint, K4);for(l=0;l<6;l++){S[l] =SL[l] + 0.166666666667*(K1[l]+2*K2[l]+2*K3[l]+K4[l]);}C = C2+C4;Stress_Hook (E,C1,C,SA);for (l=0; l<6; l++)Sv[l] = S[l]-SA[l];}Процедура «DVscStr»1void DVscStr(long Nelem, double E[6], double S[6], double DSv[6])234/* Функция, вычисляющая производную вязкой составляющейтензора напряжений в элементе5678910111.Nelem...............Номер элемента2.E...................Вектор полных деформаций элемента3.S...................Вектор полных напряжений в элементе4.DSv.................Вектор производных вязких напряжений*/{121314double SA[6];double SB[6],SBD[6],tauB;151617double ED[6];double EeB[6];1921819double EpB[6], LpB[3];double Lchain,dgamma;2021int i,j,k;2223//*********************************2425262728293031323334353637383940414243444546474849505152535455Stress_Hook (E,C1,C2,SA);for (i=0;i<6;i++)SB[i] = S[i]-SA[i];for (i=0;i<6;i++)if (i<3)SBD[i] = SB[i]-0.3333333333333333*(SB[0]+SB[1]+SB[2]);elseSBD[i] = SB[i];tauB = 0.0;for (i=0;i<6;i++)tauB += pow(SBD[i],2);tauB = sqrt(0.5*tauB);Strain_Hook_Div (SB, C4, EeB);for (i=0;i<6;i++)if (i<3)ED[i] = E[i]-0.3333333333333333*(E[0]+E[1]+E[2]);elseED[i] = E[i];for (i=0;i<6;i++)EpB[i] = ED[i]-EeB[i];PrincipalStrain (Nelem, EpB, LpB);for (i=0;i<3;i++)LpB[i] = pow((1.0+LpB[i]),2);Lchain = sqrt((LpB[0]+LpB[1]+LpB[2])/3.);dgamma = C5*(OmegaT/2/pi)/pow((Lchain-1+C8),C7)*pow(tauB,(C6-1));for (i=0;i<6;i++)if (i<3)DSv[i] = -2*C4*dgamma/sqrt(2.)*SBD[i];elseDSv[i] = -C4*dgamma/sqrt(2.)*SBD[i];}Процедура «DVscStr»1void StrainAprxn(long Nelem, double E[6], double EL[6],\19323456double alfa, double dt, double Ecalc[6], double DE[6])/* Функция, вычисляющая тензор деформаций ипроизводную полного тензора деформаций по временив момент времени tcalc=alfa*dt.

При этом предполагается,что деформации изменяются по линейному закону от элемента к элементу789101112131415161718191.Nelem...............Номер элемента2.E...................Вектор полных деформаций элемента в момент времени t3.EL..................Вектор полных деформаций элемента в момент времени t-dt4.alfa................Момент времени между t-dt и t,в который происходит вычисление производной5.dt..................Приращение времени движения элемента6.Ecalc...............Вектор деформаций в момент времени tcalc7.DE..................Производные полных напряжений в элементев момент времени tcalc*/{int i;//Переменная цикла20212223242526for (i=0;i<6;i++){Ecalc[i] = EL[i]+alfa*(E[i]-EL[i]);DE[i] = (E[i]-EL[i])/dt;}}.