{  FlexPDE  REL. 5.0.4,     VERIF. St2LH_1.pde }

  { *****************************************************
	
	Structural analysis, theory of elasticity (St)
	Two-dimensional problem  (2)
	Linear problem (L)
	Harmonic analysis (H)

               Plane stress
	Compound Double-T Region (steel and copper)
	Bottom boundary is fixed

    ***************************************************** }

  title "Harmonic Vibrations of a Double-T Region "

  variables
    Ur   { Ur,Ui are the real and imaginary parts of the displacement along axis X }
    Ui
    Vr   { Vr,Vi are the real and imaginary parts of the displacement along axis Y }
    Vi

  definitions
    Um = sqrt(Ur^2+Ui^2)     { Amplitude of the displacement along axis X }
    Vm = sqrt(Vr^2+Vi^2)     { Amplitude of the displacement along axis Y }
    Scal=0.2*globalmax(magnitude(x,y))/globalmax(magnitude(Um,Vm))

    {geometrical parameters in meters}
    h=2e-2      {thickness}
    b=16e-2	{half-heigth}
    l=5e-2	{half-width of cleat}

    nu = 0.29	{ Poisson's ratio for steel }
    E  = 2e11 	{Youngs modulus in Pascals for steel}
    rho = 7.86e3   {density in kg/m^3 for steel }
    mu  = E/2/(1+nu)	{shear module}
    lambdaL =  2*mu*nu/(1-2*nu)   {Lame's parameter}
	{Lame's parameters for plane stress }
    lambdaz = 2*mu*lambdaL/(lambdaL+2*mu)

    force=1e3

    Freq1=98.17
    Freq2=735
    Q=1000           {Q-factor }

    Fsum=Freq1+Freq2 {damping coefficients}
    alpha=pi*2*Freq1*Freq2/(Fsum*Q)
    beta = 1/(pi*2*Fsum*Q)

    omega = 2*pi*freq1 {calculates results for frequency freq1, for another frequency set another value}

  {real and imaginary parts of the strains}
    Srxx = dx(Ur)	
    Sryy = dy(Vr)
    Srxy = (dy(Ur)+dx(Vr))/2
    Sixx = dx(Ui)	
    Siyy = dy(Vi)	
    Sixy = (dy(Ui)+dx(Vi))/2	
  {real and imaginary parts of the stresses}
    Trxx= lambdaz*(Srxx+Sryy)+2*mu*Srxx-omega*beta*(lambdaz*(Sixx+Siyy)+2*mu*Sixx)
    Tryy= lambdaz*(Srxx+Sryy)+2*mu*Sryy-omega*beta*(lambdaz*(Sixx+Siyy)+2*mu*Siyy)
    Trxy=2*mu*(Srxy-omega*beta*Sixy)		
    Tixx= lambdaz*(Sixx+Siyy)+2*mu*Sixx+omega*beta*(lambdaz*(Srxx+Sryy)+2*mu*Srxx)
    Tiyy= lambdaz*(Sixx+Siyy)+2*mu*Siyy+omega*beta*( lambdaz*(Srxx+Sryy)+2*mu*Sryy)
    Tixy=2*mu*(Sixy+omega*beta*Srxy)		

  equations             	{ дифференциальные уравнения гармонических колебаний
                          для плоской задачи теории упругости при учете затухания по Релею}
     Ur:  dx(Trxx) + dy(Trxy) + rho*omega^2*Ur - omega*alpha*rho*Ui = 0
     Ui:   dx(Tixx) + dy(Tixy) + rho*omega^2*Ui + omega*alpha*rho*Ur = 0
     Vr:   dx(Trxy) + dy(Tryy) + rho*omega^2*Vr - omega*alpha*rho*Vi = 0
     Vi:   dx(Tixy) + dy(Tiyy) + rho*omega^2*Vi + omega*alpha*rho*Vr = 0

  boundaries
       Region 1        {upper region 1 from steel}
	start (h/2,0)
	natural(Ur)=0 natural(Vr)=0   natural(Ui)=0 natural(Vi)=0 {free boundary}
	line to (h/2,b-h) to (l,b-h) to (l,b)
            point load(Ur)=force  point load(Vr)=force
            line to (-l,b)
	 point load(Ur)=-force  point load(Vr)=-force
	line to (-l,b-h)
            natural(Ur)=0 natural(Vr)=0
	line to (-h/2,b-h)
 	line to (-h/2,0) to close

     Region 2         {lower region 2 with another material properties}
	{material properties for copper}
	nu=0.33  E=1.2e11  rho=8.9e3
	start (h/2,0)
	natural(Ur)=0 natural(Vr)=0   natural(Ui)=0 natural(Vi)=0   {free boundary}
	line to (h/2,-b+h) to (l,-b+h) to (l,-b)
	value(Ur)=0  value(Vr)=0  value(Ui)=0  value(Vi)=0 {fixed displacements}
	line to (-l,-b)
	natural(Ur)=0 natural(Vr)=0   natural(Ui)=0 natural(Vi)=0
	line to (-l,-b+h) to (-h/2,-b+h) to (-h/2,0) to close


  monitors
    elevation(Vr,Vi) from(-l,b) to (l,b)
      report(omega)

  plots
    grid(x+Scal*Ur,y+Scal*Vr)   as "Real displacement"
      report(freq1)
    grid(x+Scal*Ui,y+Scal*Vi)   as "Imag displacement"
      report(freq1)

  end