Commit 8b6bf8c0 authored by François Trigaux's avatar François Trigaux
Browse files

🚧 added force_scaling, works for static but not for dynamic

parent f0e1f41d
......@@ -99,6 +99,7 @@ WRITE(*,*) error
CALL Preprocess(nkp,nelem,ndof_el,member,material,frame,coord,curvature,&
& dof_con,memb_info,error)
WRITE(*,*) error
IF(error/='') RETURN
ALLOCATE(x(nsize),STAT=allo_stat)
......@@ -133,7 +134,6 @@ IF(analysis_flag==3)THEN
ENDIF
! Initial angular velocity and linear velocities of the a frame
!-----------------------------------------------------------------------
!omega_a=omega_a0;v_root_a=v_root_a0
......
......@@ -45,6 +45,7 @@ REAL(DBL)::eCabhalfL(NDIM,NDIM) ! eCab*dL/2
REAL(DBL)::eCTCabhalfL(NDIM,NDIM) ! eCTCab*dL/2
REAL(DBL)::UiDot(NDIM),ThetaDot(NDIM),CTCabPdot(NDIM),CTCabHdot(NDIM)
REAL(DBL)::force_scaling=1.0d0,mass_scaling=1.0d0
!=============================================
......@@ -140,6 +141,13 @@ IF(init_flag/=0) THEN !needed for the initial step and time marching
! WM: End calculation of ElemEqn(NDOF_ND+4:NDOF_ND+6)
!-----------------------------------------------------------------------------
ENDIF
! Divide all f_u and f_theta terms by a force scaling
ElemEqn(1:6) = ElemEqn(1:6) / force_scaling
ElemEqn(ndof_el+1:ndof_el+6) = ElemEqn(ndof_el+1:ndof_el+6) / force_scaling
IF(ndof_el==18) ElemEqn(12:18) = ElemEqn(12:18) / mass_scaling
END FUNCTION ElemEqn
!***************************************************
......@@ -363,6 +371,17 @@ IF(ndof_el==18) THEN
ENDIF
! Add force_scaling for certain entries of the jacobian
! -----------------------------------------------------
ElemJac(1:6,1:6) = ElemJac(1:6,1:6) / force_scaling ! f_(u1,th1)/(u,th)
ElemJac(ndof_el+1:ndof_el+6,1:6) = ElemJac(ndof_el+1:ndof_el+6,1:6) / force_scaling ! f_(u2,th2)/(u,th)
ElemJac(7:12,7:12) = ElemJac(7:12,7:12) * force_scaling
ElemJac(ndof_el+7:ndof_el+12,7:12) = ElemJac(ndof_el+7:ndof_el+12,7:12) * force_scaling
IF(ndof_el==18) THEN
ElemJac(1:6,13:18) = ElemJac(1:6,13:18) / force_scaling * mass_scaling
ElemJac(ndof_el+1:ndof_el+6,13:18) = ElemJac(1:6,13:18) / force_scaling * mass_scaling
ENDIF
IF(init_flag==1) THEN !needed for the initial step calculating the initial conditions
ElemJac(:,1:6)=0.D0
ElemJac(1:3,1:3)=dL*.5d0*I3
......@@ -431,6 +450,12 @@ ElemM(16:18,4:6)=-MATMUL(TRANSPOSE(eCab),((4.0D0-DOT_PRODUCT(theta,theta)*0.25D0
! WM: End calculation of ElemM(16:18,4:6)
!-----------------------------------------------------------------------------
! force_scaling and mass_scaling
elemM(1:6,1:6) = elemM(1:6,1:6) / force_scaling
elemM(18+1:18+6,1:6) = elemM(18+1:18+6,1:6) / force_scaling
elemM(1:6,13:18) = elemM(1:6,1:6) / force_scaling * mass_scaling
elemM(18+1:18+6,13:18) = elemM(18+1:18+6,1:6) / force_scaling * mass_scaling
elemM(13:18,1:6) = elemM(13:18,1:6) / mass_scaling
END SUBROUTINE ElemMass
!***************************************************
......@@ -501,8 +526,8 @@ ELSE
ENDIF
ENDIF
Fi =x_elem(7:9)
Mi =x_elem(10:12)
Fi =x_elem(7:9) * force_scaling
Mi =x_elem(10:12) * force_scaling
IF (ndof_el==18) THEN
ePi =x_elem(13:15)
......@@ -511,8 +536,8 @@ ENDIF
!eFlex(4,4)=1/(1/eFlex(4,4)+(1/eFlex(5,5)+1/eFlex(6,6))*eFlex(1,1)*Fi(1))
!write(*,*)eflex(4,4)
e1gamma=MATMUL(eFlex(1:3,:),x_elem(7:12))+e1
kappa=MATMUL(eFlex(4:6,:),x_elem(7:12))
e1gamma=MATMUL(eFlex(1:3,:),x_elem(7:12) * force_scaling)+e1
kappa=MATMUL(eFlex(4:6,:),x_elem(7:12) * force_scaling)
IF(ndof_el==18)THEN
Vi=MATMUL(eMass(1:3,:),x_elem(13:18))
......
......@@ -36,6 +36,8 @@ INTEGER::ne
INTEGER,ALLOCATABLE::irn(:),jcn(:)
REAL(DBL),ALLOCATABLE::coef(:)
REAL(DBL) :: force_scaling=1.0d0,mass_scaling=1.0d0
!=============================================
CONTAINS
......@@ -443,7 +445,7 @@ IF(dof_con(term_pt)==0) THEN
IF(pre_dof>NSTRN) THEN ! force type conditions
free_dof=pre_dof-NSTRN
rhs(nrow_kp+free_dof-1)=rhs(nrow_kp+free_dof-1)+kp_cond(j) ! prescribed conditions are external force/moments at the boundary
rhs(nrow_kp+free_dof-1)=rhs(nrow_kp+free_dof-1)+kp_cond(j) / force_scaling ! prescribed conditions are external force/moments at the boundary
rhs(nrow_kp+pre_dof-1)=rhs(nrow_kp+pre_dof-1)-flag*x_pt(j)
ELSE IF(pre_dof<=NSTRN) THEN ! displacement type conditions
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment