stress.c

边界元程序

#include "cbox11.h"

void Stress( Xp,Yp,X1,Y1,X2,Y2,dx11,dy11,dx12,dy12,
             dx22,dy22,sx11,sy11,sx12,sy12,sx22,sy22)
float Xp,Yp,X1,Y1,X2,Y2,*dx11,*dy11,*dx12,*dy12;
float *dx22,*dy22,*sx11,*sy11,*sx12,*sy12,*sx22,*sy22; 
{
	float Xg[5],Yg[5],Z[5],W[5]; 	/*[dimension increased by 1]*/
	float Ax,Bx,Ay,By,nx,ny,slope,Perp,sgn,SR,FA,AL,Denom,rx,ry,Ra;
	int i;

	Z[1] = 0.86113631;
	Z[2] = -Z[1];
	Z[3] = 0.33998104;
	Z[4] = -Z[3];
	W[1] = 0.34785485;
	W[2] = W[1];
	W[3] = 0.65214515;
	W[4] = W[3];
	
	Ax = (X2-X1)/2.0;
	Bx = (X2+X1)/2.0;
	Ay = (Y2-Y1)/2.0;
	By = (Y2+Y1)/2.0;
	SR = sqrt(SQ(Ax)+SQ(Ay));
	nx = (Y2-Y1)/(2*SR);
	ny = (X1-X2)/(2*SR);
	
	if( Ax ) 
	{
		slope = Ay/Ax;
		Perp = fabs((slope*Xp-Yp+Y1-slope*X1)/sqrt(SQ(slope)+1));
	}
	else
		Perp = fabs(Xp-X1);
        
	/*[Determine the direction of the outward normal ]*/

	sgn = (X1-Xp)*(Y2-Yp)-(X2-Xp)*(Y1-Yp);

	if ( sgn < 0 ) 
		Perp = -Perp;
        
	(*dx11) = 0.0;
	(*dy11) = 0.0;
	(*dx12) = 0.0;
	(*dy12) = 0.0;
	(*dx22) = 0.0;
	(*dy22) = 0.0;
	(*sx11) = 0.0;
	(*sy11) = 0.0;
	(*sx12) = 0.0;
	(*sy12) = 0.0;
	(*sx22) = 0.0;
	(*sy22) = 0.0;

	/*[ Compute displacement and stress coefficients ]*/
	

	FA = 1.0-4.0*nu;
	AL = 1.0-2.0*nu;
	Denom = 4.0*pi*(1.0-nu); 	
  
         
	for(i=1;i<=4;i++)
	{
		Xg[i] = Ax*Z[i]+Bx;
		Yg[i] = Ay*Z[i]+By;
		Ra  = sqrt(SQ(Xp-Xg[i])+SQ(Yp-Yg[i]));
		rx = (Xg[i]-Xp)/Ra;
		ry = (Yg[i]-Yp)/Ra;
		(*dx11) += (AL*rx+2*cube(rx))*W[i]*SR/(Denom*Ra);
		(*dy11) += (2*SQ(rx)*ry-AL*ry)*W[i]*SR/(Denom*Ra);
		(*dx12) += (AL*ry+2*(SQ(rx))*ry)/(Denom*Ra)*W[i]*SR;
		(*dy12) += (AL*rx+2*rx*SQ(ry))/(Denom*Ra)*W[i]*SR;
		(*dx22) += (2*rx*SQ(ry)-AL*rx)/(Denom*Ra)*W[i]*SR;
		(*dy22) += (AL*ry+2*cube(ry))/(Denom*Ra)*W[i]*SR;
		(*sx11) += (2*Perp/Ra*(AL*rx+nu*2*rx-4*cube(rx))+
                       4*nu*nx*SQ(rx)+AL*(2*nx*SQ(rx)+2*nx)-
                       FA*nx)*2*mu/(Denom*SQ(Ra))*W[i]*SR;
		(*sy11) += (2*Perp/Ra*(AL*ry-4*SQ(rx)*ry)+
                       4*nu*nx*rx*ry+AL*2*ny*SQ(rx)-
                       FA*ny)*2*mu/(Denom*SQ(Ra))*W[i]*SR;
		(*sx12) += (2*Perp/Ra*(nu*ry-4*SQ(rx)*ry)+2*nu*
                       (nx*ry*rx+ny*SQ(rx))+AL*(2*nx*rx*
                       ry+ny))*2*mu/(Denom*SQ(Ra))*W[i]*SR;
		(*sy12) += (2*Perp/Ra*(nu*rx-4*rx*SQ(ry))+2*nu*
                       (nx*SQ(ry)+ny*rx*ry)+AL*(2*ny*rx*ry+
                       nx))*2*mu/(Denom*SQ(Ra))*W[i]*SR;
		(*sx22) += (2*Perp/Ra*(AL*rx-4*rx*SQ(ry))+4*nu*
                       ny*rx*ry+AL*2*nx*SQ(ry)-FA*nx)*2*mu/
                       (Denom*SQ(Ra))*W[i]*SR;
		(*sy22) += (2*Perp/Ra*(AL*ry+2*nu*ry-4*cube(ry))+
                       4*nu*ny*SQ(ry)+AL*(2*ny*SQ(ry)+2*ny)-
                       FA*ny)*2*mu/(Denom*SQ(Ra))*W[i]*SR;

     }   


 }