Demo entry 5853376

Hydrostatik

   

Submitted by anonymous on Jul 20, 2016 at 12:31
Language: Fortran. Code size: 9.1 kB.

subroutine hydrostat_zyl_sub(awl,nwl,azy,ezy,r,L,V,SP)

integer, parameter           :: rk=4 !für double precision
!Eingabedaten
!!!!!!!!!!!!!!!!!!
!Wasserfläche
real(kind=rk),dimension(3),intent(in) :: awl    !Aufpunktvektor
real(kind=rk),dimension(3),intent(in) :: nwl    !Normalenvektor positiv von der Wasseroberfläche weg
!Zylinder
real(kind=rk),dimension(3),intent(in) :: azy    !Aufpunktvektor
real(kind=rk),dimension(3),intent(in) :: ezy    !Erzeugender Vektor positiv ins Zylinder-Innere
real(kind=rk),intent(in)              :: r      !Radius Zylinder
real(kind=rk),intent(in)              :: L      !Länge Zylinder

!Berechnungsvariablen
real(kind=rk),dimension(3) :: f             !Zylinderfußpunkt
real(kind=rk),dimension(3) :: k             !Zylinderkopfpunkt
real(kind=rk),dimension(3) :: nzy           !Normalenvektor Zylinder (vom Fußpunkt ins Innere positiv)
real(kind=rk),dimension(3) :: n0zy          !Normaleneinheitsvektor Zylinder
real(kind=rk),dimension(3) :: n0wl          !Normaleneinheitsvektor Wasserfläche
real(kind=rk)              :: s             !Skalarprodukt <nwl,ezy> für Abfrage Orientierung des Zylinders
real(kind=rk)              :: d             !Abstand Zylinderfußpunkt zur Wasserfläche
real(kind=rk)              :: e             !Abstand Zylinderkopfpunkt zur Wasserfläche
real(kind=rk)              :: h1,h2         !Höhe für Zylinderhuf
real(kind=rk)              :: alpha         !Winkel zwischen der Wasserfläche und der Grundfläche des Zylinders
real(kind=rk)              :: beta          !Winkel für Tankproblem
real(kind=rk)              :: g,g1,g2       !Strecke für Zylinderhuf
real(kind=rk)              :: b,b1,b2       !Strecke für Zylinderhuf vgl. Bronstein S.160
real(kind=rk)              :: phi,phi1,phi2 !Winkel für Zylinderhuf
real(kind=rk)              :: V1,V2         !Teilvolumen Zylinderhuf
real(kind=rk),intent(out)  :: V             !Zylindervolumen

!Berechnungsvariablen zusätzlich für den Schwerpunkt
real(kind=rk),dimension(4) :: nwl2          !transformierter Normalenvektor Wasserfläche
real(kind=rk),dimension(4) :: n0wl2         !transformierter Normaleneinheitsvektor Wasserfläche
real(kind=rk),dimension(4,4):: Mt           !Transformationsmatrix um nzy in die z-Achse des Koordinatensystems zu verschieben
real(kind=rk),dimension(4,4):: Mtinv        !Inverse Transformationsmatrix
real(kind=rk)              :: xls,yls,zls   !Schwerpunktkoordinaten im lokalen System
real(kind=rk),dimension(3) :: xl,yl,zl      !Lokale Koordinatenvektoren (Orthogonalbasis)
real(kind=rk),dimension(3) ::xl0,yl0,zl0    !Lokale Koordinateneinheitsvektoren (Orthonormalbasis)
real(kind=rk),dimension(4)  :: SP4          !Schwerpunktvektor Dimension 4
real(kind=rk),dimension(3),intent(out) :: SP !Schwerpunkt

!später aus e4 holen
real(kind=rk),parameter    :: pi=3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679   

!Normaleneinheitsvektor WL
n0wl=nwl/ (norm2(nwl))

!Orientierung des Zylinders zur Wasserfläche bestimmen
s=dot_product(nwl,ezy)

skalarprodukt: if (s==0) then    !Zylinderachse und Wasseroberfläche parallel (Tankproblem)
  d=dot_product(azy,n0wl)-dot_product(awl,n0wl) !Abstand zur Wasserlinie 
 abstand: if (d<-r) then !Zylinder komplett getaucht
  V=pi*L*r**2.
  SP=azy+(L/2.)*(ezy/norm2(ezy))
  return
 elseif (abs(d)<r) then !Zylinder teilweise getaucht
 beta=2.*acos(d/r)
 V=L*(r**2/2.)*(beta-sin(beta))
  abstand2: if (d>=0) then !weniger als der halbe Zylinder getaucht
   SP=azy+(L/2.)*(ezy/norm2(ezy))-((2*sqrt(r**2-d**2))**3*L)/(12.*V)*n0wl
   return
  elseif (d<0) then !mehr als der halbe Zylinder getaucht
    SP=azy+(L/2.)*(ezy/norm2(ezy))- ((2*sqrt(r**2-d**2))**3)/(12.*pi*r**2))*n0wl  !(pi*r**2*L**2/2.)*n0zy-(V-pi*r**2*L)*(L/2.*n0zy+((2*sqrt(2*r*h-h**2))**3*L)/(12.*V)*n0wl)*(1./V)
   return
  endif abstand2
 elseif(d>r) then !Zylinder komplett ausgetaucht
 write(*,*) "Kein getauchtes Volumen!"
  return
 endif abstand
elseif (s>0) then !Zylinder "normal" orientiert
 f=azy
 nzy=ezy
 n0zy=nzy/(norm2(nzy))
 k=azy+L*n0zy
elseif(s<0) then !Zylinder ist entgegengesetzt orientiert
 f=azy+L*(ezy/(norm2(ezy)))
 k=azy
 nzy=(-1)*ezy
 n0zy=nzy/(norm2(nzy))
endif skalarprodukt

d=dot_product(f,n0wl)-dot_product(awl,n0wl) !Abstand Zylinderfußpunkt zur Wasserfläche
e=dot_product(k,n0wl)-dot_product(awl,n0wl) !Abstand Zylinderkopfpunkt zur Wasserfläche
alpha = acos(dot_product(n0wl,n0zy))

winkel: if (abs(alpha)>0) then  
!Orthogonalbasis lokales Koordinatensystem
zl=nzy
yl(1)=zl(2)*nwl(3)-zl(3)*nwl(2)
yl(2)=zl(3)*nwl(1)-zl(1)*nwl(3)
yl(3)=zl(1)*nwl(2)-zl(2)*nwl(1)
xl(1)=yl(2)*zl(3)-yl(3)*zl(2)
xl(2)=yl(3)*zl(1)-yl(1)*zl(3)
xl(3)=yl(1)*zl(2)-yl(2)*zl(1)
!Orthonormalbasis lokales Koordinatensystem
xl0=xl/norm2(xl)
yl0=yl/norm2(yl)
zl0=zl/norm2(zl)
!Transformationsmatrix
Mt(1,:)=[xl0, 0.]
Mt(2,:)=[yl0, 0.]
Mt(3,:)=[zl0, 0.]
Mt(4,:)=[0.,0.,0.,1.]
!Inverse der Transformationsmatrix
Mtinv(:,1)=[xl0,0.]
Mtinv(:,2)=[yl0,0.]
Mtinv(:,3)=[zl0,0.]
Mtinv(:,4)=[f,1.]

abstand_d_e: if (d>=r*sin(alpha) .and. e>=r*sin(alpha)) then !Zylinder komplett ausgetaucht
 write(*,*) "Kein getauchtes Volumen!"
 return
elseif (-r*sin(alpha)>=d .and. -r*sin(alpha)>=e) then !Zylinder ist komplett eingetaucht
 V=pi*L*r**2.
 SP=f+(L/2.)*n0zy
 return
elseif (d<=-r*sin(alpha) .and. e>=r*sin(alpha)) then !schief abgeschnittener Zylinder
 V=pi*(abs(d)/cos(alpha))*r**2.
 
 !Wasserfläche transformieren
 nwl2=matmul(Mt,[n0wl, 1.])
 n0wl2 = [nwl2(1:3)/norm2(nwl2(1:3)),1.] !durch die Transformation ist die Länge nicht mehr 1
 xls=-(r**2*n0wl2(1))/(4.*abs(d))
 yls=-(r**2*n0wl2(2))/(4.*abs(d))
 zls=abs(d)/(2.*n0wl2(3))+r**2/(8.*n0wl2(3)*abs(d))*(n0wl2(1)**2+n0wl2(2)**2)
 !Rücktransformation
 SP4=matmul(Mtinv,[xls,yls,zls, 1.])
 SP=SP4(1:3)
 return
elseif (d<=-r*sin(alpha) .and. abs(e)<r*sin(alpha)) then !kleine Ecke ausgetaucht
 g=e/sin(alpha)
 b=r+g
 h=b*tan(alpha)
 phi=acos(-g/r)
 !Volumen
 V1=((h*r**3.)/(b))*(sin(phi)-((sin(phi)**3.)/3.)-phi*cos(phi))
 V=pi*L*r**2.-V1
 !Lokale Schwerpunktkoordinaten
 xls=-((2.*h)/(b)*(r**4/8.*(pi/2.-asin((r-b)/b))-(r**2*(r-b)*sqrt(2.*r*b-b**2))/8.+((b-r)*(2.*r*b-b**2)**(3./2.))/12.))/(V)
 yls=0.
 zls=(L**2/2.*pi*r**2-(L-(h**2)/(V1*b**2)*((pi*r**2)/4.*(r**2/4.-2.*r*b+r**2+b**2)+(5.*(b-r)*(2.*r*b-b**2)**(3./2.))/12.& 
  & +(-(5.*r**3)/8.+(13.*r**2*b)/8.-(3.*r*b**2)/2.+b**3/2.)*sqrt(2.*r*b-b**2)+(r**3*b-(5.*r**4)/8.-(b**2*r**2)/2.)*asin((r-b)/r)))*V1)/(V)
 !Rücktransformation
 SP4=matmul(Mtinv,[xls,yls,zls, 1.])
 SP=SP4(1:3)
 return
elseif (abs(d)<r*sin(alpha) .and. e>=r*sin(alpha)) then !kleine Ecke eingetaucht
 g=d/sin(alpha)
 b=r-g
 h=b*tan(alpha)
 phi=acos(g/r)
 !Volumen
 V=((h*r**3.)/(b))*(sin(phi)-((sin(phi)**3.)/3.)-phi*cos(phi))
 !Lokale Schwerpunktkoordinaten
 xls=-(2.*h)/(V*b)*(r**4/8.*(pi/2.-asin((r-b)/b))-(r**2*(r-b)*sqrt(2.*r*b-b**2))/8.+((b-r)*(2.*r*b-b**2)**(3./2.))/12.)
 yls=0.
 zls=(h**2)/(V*b**2)*((pi*r**2)/4.*(r**2/4.-2.*r*b+r**2+b**2)+(5.*(b-r)*(2.*r*b-b**2)**(3./2.))/12.+(-(5.*r**3)/8.+(13.*r**2*b)/8. &
              & -(3.*r*b**2)/2.+b**3/2.)*sqrt(2.*r*b-b**2)+(r**3*b-(5.*r**4)/8.-(b**2*r**2)/2.)*asin((r-b)/r))
 !Rücktransformation
 SP4=matmul(Mtinv,[xls,yls,zls, 1.])
 SP=SP4(1:3)
 return
elseif (abs(d)<r*sin(alpha) .and. abs(e)<r*sin(alpha)) then !Zylinderhuf minus Zylinderhuf
 !großer Zylinderhuf
 g1=d/sin(alpha)
 b1=r-g1
 h1=b1*tan(alpha)
 phi1=acos(g1/r)
 !kleiner Zylinderhuf
 g2=e/sin(alpha)
 b2=r-g2
 h2=b2*tan(alpha)
 phi2=acos(g2/r)
 !Volumen
 V1=(h1*r**3.)/(b1))* (sin(phi1)-((sin(phi1)**3.)/3.)-phi1*cos(phi1))
 V2=((h2*r**3.)/(b2))*(sin(phi2)-((sin(phi2)**3.)/3.)-phi2*cos(phi2))
 V=V1-V2
 !lokale Schwerpunktkoordinaten
 xls=(-(2.*h1)/(b1)*(r**4/8.*(pi/2.-asin((r-b1)/b1))-(r**2*(r-b1)*sqrt(2.*r*b1-b1**2))/8.+((b1-r1)*(2.*r*b1-b1**2)**(3./2.))/12.) &
  & +(2.*h2)/(b2)*(r**4/8.*(pi/2.-asin((r-b2)/b2))-(r**2*(r-b2)*sqrt(2.*r*b2-b2**2))/8.+((b2-r1)*(2.*r*b2-b2**2)**(3./2.))/12.))/(V)
 yls=0.
 zls=((h1**2)/(b1**2)*((pi*r**2)/4.*(r**2/4.-2.*r*b1+r**2+b1**2)+(5.*(b1-r)*(2.*r*b1-b1**2)**(3./2.))/12.+(-(5.*r**3)/8.+(13.*r**2*b1)/8. &
   & -(3.*r*b1**2)/2.+b1**3/2.)*sqrt(2.*r*b1-b1**2)+(r**3*b1-(5.*r**4)/8.-(b1**2*r**2)/2.)*asin((r-b1)/r))- &
   & (L+(h2**2)/(V2*b2**2)*((pi*r**2)/4.*(r**2/4.-2.*r*b2+r**2+b2**2)+(5.*(b2-r)*(2.*r*b2-b2**2)**(3./2.))/12.+(-(5.*r**3)/8.+(13.*r**2*b2)/8. &
   & -(3.*r*b2**2)/2.+b2**3/2.)*sqrt(2.*r*b2-b2**2)+(r**3*b2-(5.*r**4)/8.-(b2**2*r**2)/2.)*asin((r-b2)/r)))*V2)/(V)
 !Rücktransformation
 SP4=matmul(Mtinv,[xls,yls,zls, 1.])
 SP=SP4(1:3)
 return
endif abstand_d_e

else
abstand_d: if (d>=0 .and. e>=0) then !Zylinder komplett ausgetaucht
 write(*,*) "Kein getauchtes Volumen!"
 return
elseif (d<0 .and. e<0) then !Zylinder ist komplett eingetaucht
 V=pi*L*r**2.
 SP=f+(L/2.)*n0zy
 return
elseif (d<0 .and. e>0) then !Zylinder teilweise getaucht
 V=pi*(abs(d)/cos(alpha))*r**2.
 SP=f+(abs(d)/2.)*n0zy
 return
endif abstand_d
endif winkel

end subroutine hydrostat_zyl_sub

This snippet took 0.02 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).