PlaneRigidFrameCalcer

最新推荐文章于 2022-02-25 23:45:16 发布
dTianx
最新推荐文章于 2022-02-25 23:45:16 发布
阅读量1.9k
点赞数
本文链接：https://blog.csdn.net/dtianx/article/details/244351
版权
Visual Basic / Visual Basic.NET 专栏收录该内容
1 篇文章 0 订阅
订阅专栏
1Option Base 1
2Option Explicit
3
4Dim NodeCount As Integer
5Dim ElementCount As Integer
6Dim UnknownNodeDisplacementCount As Integer
7
8Dim NodeLoadCount As Integer
9Dim NonNodeLoadCount As Integer
10
11Dim ElementNodeNoArray() As Integer
12Dim NodeDisplacementArray() As Integer
13Dim NodeDisplacementNoArray() As Integer
14
15Dim EAArray() As Double
16Dim EIArray() As Double
17
18Dim XCoordArray() As Double
19Dim YCoordArray() As Double
20
21Dim NodeLoadArray() As Double
22Dim NonNodeLoadArray() As Double
23
24Dim ElementIndex As Integer
25Dim NodeSumLoadArray() As Double
26
27Dim ElementLength As Double
28Dim SinValue As Double
29Dim CosValue As Double
30Dim DisplacementNo As Integer
31
32Dim StiffnessMatrixGloble() As Double
33Dim ElementSupportForceArray(6) As Double
34Dim ElementCoordConvMatrix(6, 6) As Double
35Dim ElementLocaVecter(6) As Integer
36Dim ElementEquareNodeLoad(6) As Double
37Dim ElementStiffnessGlobleArray(6, 6) As Double
38Dim ElementStiffnessLocalArray(6, 6) As Double
39
40'all consts of a single element
41Function GetElementConsts(ElementNo As Integer, _
42                        ByRef ElementLength As Double, _
43                        ByRef SinValue As Double, _
44                        ByRef CosValue As Double _
45                        )
46
47    Dim NodeStartNo, NodeEndNo As Integer
48    Dim XLenDiff, YLenDiff As Double
49    
50    'first, we got to know which two node(No.) formed the element
51    NodeStartNo = ElementNodeNoArray(1, ElementNo)
52    NodeEndNo = ElementNodeNoArray(2, ElementNo)
53    
54    'got the X and Y coords, then we can get its length
55    XLenDiff = XCoordArray(NodeStartNo) - XCoordArray(NodeEndNo)
56    YLenDiff = YCoordArray(NodeStartNo) - YCoordArray(NodeEndNo)
57    
58    'All three to be ret
59    ElementLength = Sqr(XLenDiff ^ 2 + YLenDiff ^ 2)
60    SinValue = YLenDiff / ElementLength
61    CosValue = XLenDiff / ElementLength
62
63End Function
64
65'Form the Element location vetors
66Function FormElementLocVetor(ElementNo As Integer, _
67                            ByRef ElementLocVetor() As Integer _
68                            )
69    Dim NodeStartNo, NodeEndNo As Integer
70        
71    'first, we got to know which two node(No.) formed the element
72    NodeStartNo = ElementNodeNoArray(1, ElementNo)
73    NodeEndNo = ElementNodeNoArray(2, ElementNo)
74    
75    Dim index As Integer
76    For index = 1 To 3
77        'first 3 items must be the displacement of the start node, then the end node
78        ElementLocVetor(index) = NodeDisplacementNoArray(index, NodeStartNo)
79        ElementLocVetor(index + 3) = NodeDisplacementNoArray(index, NodeEndNo)
80    Next index
81                                
82End Function
83                            
84'Form Element Stiffness Matrix
85Function FormElementStiffnessMatrix(ElementNo As Integer, _
86                                    ByRef ElemStifMatrix() As Double, _
87                                    ElementLength As Double _
88                                    )
89    Dim G, G1, G2, G3 As Double
90    G = EA(ElementNo) / ElementLength
91    G1 = 2# * EI(ElementNo) / ElementLength
92    G2 = 3# * G1 / ElementLength
93    G3 = 2# * G2 / ElementLength
94    
95    ElemStifMatrix(1, 1) = G
96    ElemStifMatrix(1, 4) = -G
97    ElemStifMatrix(4, 4) = G
98    ElemStifMatrix(2, 2) = G3
99    ElemStifMatrix(5, 5) = G3
100    ElemStifMatrix(2, 3) = G2
101    ElemStifMatrix(2, 6) = G2
102    ElemStifMatrix(3, 5) = -G2
103    ElemStifMatrix(5, 6) = -G2
104    ElemStifMatrix(3, 3) = 2# * G1
105    ElemStifMatrix(6, 6) = 2# * G1
106    ElemStifMatrix(3, 6) = G1
107    
108    'do item exchanging
109    Dim index, subindex As Integer
110    For index = 1 To 5
111        For subindex = 1 + index To 6
112            ElemStifMatrix(subindex, index) = ElemStifMatrix(index, subindex)
113        Next subindex
114    Next index
115    
116End Function
117
118'form Element coords convert matrix
119Function FormElementCoordConvMatrix(SinValue As Double, _
120                                    CosValue As Double, _
121                                    ByRef ElementCoordConvMatrix() As Double _
122                                    )
123    ElementCoordConvMatrix(1, 1) = CosValue
124    ElementCoordConvMatrix(1, 2) = SinValue
125    
126    ElementCoordConvMatrix(2, 1) = -SinValue
127    ElementCoordConvMatrix(2, 2) = CosValue
128    
129    ElementCoordConvMatrix(3, 3) = 1#
130    
131    ElementCoordConvMatrix(4, 4) = CosValue
132    ElementCoordConvMatrix(4, 5) = SinValue
133    
134    ElementCoordConvMatrix(5, 4) = -SinValue
135    ElementCoordConvMatrix(5, 5) = CosValue
136    
137End Function
138
139'form Element Force on
140Function ElementSupportForce(NonNodeLoadNo As Integer, _
141                            ElementLength As Double, _
142                            ByRef ElementSupportForceArray() As Double _
143                            )
144    Dim LoadType As Integer
145    Dim Offset As Double
146    Dim LoadValue As Double
147    
148    LoadType = NonNodeLoadArray(2, NonNodeLoadNo)
149    Offset = NonNodeLoadArray(3, NonNodeLoadNo)
150    LoadValue = NonNodeLoadArray(4, NonNodeLoadNo)
151    
152    Dim C, G, B, S As Double
153    
154    C = Offset / ElementLength  'a/l
155    G = C ^ 2                   '(a/l)^2
156    B = ElementLength - Offset  'l - a
157    
158    Select Case LoadType
159        Case 10
160            S = LoadValue * Offset * 0.5                                '(q * a) / 2
161            ElementSupportForceArray(2) = -S * (2# - 2# * G + C * G)    ' -(qa/2) * (2 - 2 * (a/l)^2 + a/l * (a/l)^3 //right!
162            ElementSupportForceArray(5) = -S * G * (2# - C)             '-(qa/2) * (a/l)^2 * (2 - a/l)  //right
163            S = S * Offset / 6#                                         '(qa^2)/2 / 6 => (qa^2) / 12
164            ElementSupportForceArray(3) = -S * (6# - 8# * C + 3# * G)   ' - [ (qa^2) / 12 ] * ( 6 - 8a/l + 3a^2/l^2 )     //right
165            ElementSupportForceArray(6) = S * C * (4# - 3# * C)         '[ (qa^2) / 12 ] * a/l * ( 4 - 3a/l)    //good
166        'All above is right to P59 => Table 10 - 1
167        
168        Case 20
169            S = B / ElementLength
170            ElementSupportForceArray(2) = -LoadValue * S * S * (1# + 2# * C)    '-q * [(l-a)/l]^2 * ( 1 + 2 * (a/l)^2)
171            ElementSupportForceArray(5) = -LoadValue * G * (1# + 2# * S)        '-q * (a/l)^2 * [ 1 + 2 * (l-a)/l ]
172            ElementSupportForceArray(3) = -LoadValue * S * S * Offset           '-q * (l-a)/l * a
173            ElementSupportForceArray(6) = LoadValue * B * G                     'q * (l-a) * (a/l)^2
174        
175        Case 30
176            S = B / ElementLength
177            ElementSupportForceArray(2) = 6# * LoadValue * C * S / ElementLength
178            ElementSupportForceArray(5) = -ElementSupportForceArray(2)
179            ElementSupportForceArray(3) = LoadValue * S * (2# - 3# * S)
180            ElementSupportForceArray(6) = LoadValue * C * (2# - 3# * C)
181        
182        Case 40
183            S = LoadValue * Offset * 0.25
184            ElementSupportForceArray(2) = -S * (2# - 3# * G + 1.6 * G * C)
185            ElementSupportForceArray(5) = -S * G(3# - 1.6 * C)
186            S = S * Offset
187            ElementSupportForceArray(3) = -S * (2# - 3# * C + 1.2 * G) / 1.5
188            ElementSupportForceArray(6) = S * C * (1# - 0.8 * C)
189            
190        Case 50
191            ElementSupportForceArray(1) = -LoadValue * Offset * (1# - 0.5 * C)
192            ElementSupportForceArray(4) = -0.5 * LoadValue * C * Offset
193        
194        Case 60
195            ElementSupportForceArray(1) = -LoadValue * B / ElementLength
196            ElementSupportForceArray(4) = -LoadValue * C
197            
198        Case 70
199            S = B / ElementLength
200            ElementSupportForceArray(2) = LoadValue * G * (3# * S + C)
201            ElementSupportForceArray(5) = -ElementSupportForceArray(2)
202            S = S * B / ElementLength
203            ElementSupportForceArray(3) = -LoadValue * S * Offset
204            ElementSupportForceArray(6) = LoadValue * G * B
205        Case Default
206    End Select
207    
208End Function
209
210
211Private Sub Class_Initialize()
212    NodeCount = 0
213    ElementCount = 0
214    UnknownNodeDisplacementCount = 0
215
216    NodeLoadCount = 0
217    NonNodeLoadCount = 0
218    Dim X As Integer
219    
220    ReDim XCoordArray(255) As Double
221    ReDim YCoordArray(255) As Double
222    
223    ReDim ElementNodeNoArray(2, 255) As Integer
224    ReDim EAArray(255) As Double
225    ReDim EIArray(255) As Double
226    ReDim NodeLoadArray(2, 255) As Double
227    ReDim NodeSumLoadArray(255 * 3) As Double
228    
229    ReDim NodeLoadArray(2, 255) As Double
230    ReDim NonNodeLoadArray(4, 255) As Double
231    
232    ReDim NodeDisplacementArray(3, 255) As Integer
233    ReDim NodeDisplacementNoArray(3, 255) As Integer
234    
235End Sub
236
237Public Function AddNode(X As Double, Y As Double) As Integer
238    NodeCount = NodeCount + 1
239    If NodeCount > 255 Then
240        ReDim Preserve XCoordArray(NodeCount) As Double
241        ReDim Preserve YCoordArray(NodeCount) As Double
242        ReDim Preserve NodeDisplacementNoArray(3, NodeCount) As Integer
243        ReDim Preserve NodeSumLoadArray(NodeCount * 3) As Double
244    End If
245    XCoordArray(NodeCount) = X
246    YCoordArray(NodeCount) = Y
247    AddNode = NodeCount
248End Function
249
250
251
252Public Function AddElement(StartNodeNo As Integer, EndNodeNo As Integer) As Integer
253    ElementCount = ElementCount + 1
254    If ElementCount > 255 Then
255        ReDim Preserve ElementNodeNoArray(2, ElementCount) As Integer
256    End If
257    ElementNodeNoArray(1, ElementCount) = StartNodeNo
258    ElementNodeNoArray(2, ElementCount) = EndNodeNo
259    
260    AddElement = ElementCount
261End Function
262
263Public Function AddElementLoad(ElementNo As Integer, LoadType As Integer, Offset As Double, LoadValue As Double) As Integer
264    NonNodeLoadCount = NonNodeLoadCount + 1
265    If NonNodeLoadCount > 255 Then
266        ReDim Preserve NonNodeLoadArray(4, NonNodeLoadCount) As Double
267    End If
268    NonNodeLoadArray(1, NonNodeLoadCount) = ElementNo
269    NonNodeLoadArray(2, NonNodeLoadCount) = LoadType
270    NonNodeLoadArray(3, NonNodeLoadCount) = Offset
271    NonNodeLoadArray(4, NonNodeLoadCount) = LoadValue
272    NonNodeLoadCount = NonNodeLoadCount
273End Function
274
275Public Function AddNodeLoad(DisplacementNo As Integer, LoadValue As Double) As Integer
276    NodeLoadCount = NodeLoadCount + 1
277    If NodeLoadCount > 255 Then
278        ReDim Preserve NodeLoadArray(2, NodeLoadCount) As Double
279    End If
280    
281    NodeLoadArray(1, NodeLoadCount) = DisplacementNo
282    NodeLoadArray(2, NodeLoadCount) = LoadValue
283    AddNodeLoad = NodeLoadCount
284End Function
285
286Public Function SetElementEAandEI(ElementNo As Integer, EA As Double, EI As Double) As Integer
287    If ElementNo > 255 Then
288        ReDim Preserve EIArray(ElementCount) As Double
289        ReDim Preserve EAArray(ElementCount) As Double
290    End If
291    
292    EAArray(ElementCount) = EA
293    EIArray(ElementCount) = EI
294    
295    SetElementEAandEI = ElementNo
296    
297End Function
298
299Public Function GetNodesCount() As Integer
300    GetNodesCount = NodeCount
301End Function
302
303Public Function GetElementCount() As Integer
304    GetElementCount = ElementCount
305End Function
306
307
308Public Function GetNodeLoadCount() As Integer
309    GetNodeLoadCount = NodeLoadCount
310End Function
311
312Public Function GetNonNodeLoadCount() As Integer
313    GetNonNodeLoadCount = NonNodeLoadCount
314End Function
315Public Function GetResult()
316    Dim index As Integer
317    'process all noteloads
318    For index = 1 To NodeLoadCount
319        'Get the Displacement Id of the nodeload
320        DisplacementNo = NodeLoadArray(1, index)
321        NodeSumLoadArray(DisplacementNo) = NodeLoadArray(2, index)
322    Next index
323    
324    'then all nonnodeloads (element load)
325    For index = 1 To NonNodeLoadCount
326        'which element the load is in
327        ElementIndex = NonNodeLoadArray(1, index)
328        'first, work out all const needed, such as ElementLength, sin, cos
329        GetElementConsts ElementIndex, ElementLength, SinValue, CosValue
330        'Get the support forces, as follow
331        'StartNode
332        '   |_____Force on X
333        '   |_____Force on Y
334        '   |_____Moment
335        ElementSupportForce index, ElementLength, ElementSupportForceArray
336    
337        'this forces and matrixes are in local crood, turn it to globle
338        FormElementCoordConvMatrix SinValue, CosValue, ElementCoordConvMatrix
339    
340        'After that, the forces and moments VALUE can be globle
341        'Then, we should put the matrix into the globle one, but
342        'we do not know where to place it
343        'now what we need may be a Element Localtion Vetor
344        FormElementLocVetor ElementIndex, ElementLocaVecter
345        'P = - ElementCoordConvMatrix ^T * ElementSupportForceArray`
346        Dim ElementCoordConvMatrix_T() As Double
347        MatrixTranspose ElementCoordConvMatrix, ElementCoordConvMatrix_T
348        MatrixMultiply ElementCoordConvMatrix_T, ElementSupportForceArray, ElementEquareNodeLoad
349        MatrixChangeSign ElementEquareNodeLoad
350        'Here, we got equare effect node load
351        
352        Dim index1  As Integer
353        For index1 = 1 To 6
354            If ElementLocaVecter(index1) <> 0 Then NodeSumLoadArray(ElementLocaVecter(index1)) = NodeSumLoadArray(ElementLocaVecter(index1)) + ElementEquareNodeLoad(index1)
355        Next index1
356
357    Next index
358    
359    For ElementIndex = 1 To GetElementCount
360        GetElementConsts ElementIndex, ElementLength, SinValue, CosValue
361        FormElementCoordConvMatrix SinValue, CosValue, ElementCoordConvMatrix
362        FormElementStiffnessMatrix ElementIndex, ElementStiffnessLocalArray, ElementLength
363        FormElementLocVetor ElementIndex, ElementLocaVecter
364        Dim ElementCoordConvMatrix_T2(6) As Double
365        Me.MatrixTranspose ElementCoordConvMatrix, ElementCoordConvMatrix_T2
366        Dim TempMatrix(6) As Double
367        MatrixMultiply ElementStiffnessLocalArray, ElementCoordConvMatrix, TempMatrix
368        MatrixMultiply ElementCoordConvMatrix_T2, TempMatrix, ElementStiffnessGlobleArray
369        
370        Dim index2 As Integer
371        For index2 = 1 To 6
372            Dim CurrentLocIndex As Integer
373            CurrentLocIndex = ElementLocaVecter(index2)
374            If (CurrentLocIndex = 0) Then GoTo NextIndex2 'just skip, like cycle in fortran or continue in cpp
375            Dim index3 As Integer
376            For index3 = 1 To 6
377                Dim NonzeroVetorIndex As Integer
378                NonzeroVetorIndex = ElementLocaVecter(index3)
379                If NonzeroVetorIndex = 0 Or NonzeroVetorIndex < index2 Then GoTo NextIndex3
380                
381                Dim JJ As Integer
382                JJ = NonzeroVetorIndex - CurrentLocIndex + 1
383                
384                StiffnessMatrixGloble(CurrentLocIndex, JJ) = StiffnessMatrixGloble(CurrentLocIndex, JJ) + ElementStiffnessGlobleArray(index2, index3)
385NextIndex3:
386            Next index3
387NextIndex2:
388        Next index2
389        
390    Next ElementIndex
391
392    Dim K As Integer
393    For K = 1 To GetUnknownNodeDisplacementCount
394        Dim IM As Integer
395        IM = K + NW - 1
396        If GetUnknownNodeDisplacementCount < IM Then IM = GetUnknownNodeDisplacementCount
397        
398        Dim I As Integer
399        For I = K + 1 To IM
400            
401            Dim L As Integer
402            L = I - K + 1
403            
404            Dim C As Double
405            C = StiffnessMatrixGloble(K, L) / StiffnessMatrixGloble(K, 1)
406            Dim JM As Integer
407            JM = NW - L + 1
408            
409            Dim index As Integer
410            For index = 1 To JM
411                StiffnessMatrix(I, index) = StiffnessMatrix(I, index) - C * StiffnessMatrix(K, index + I - K)
412                NodeSumLoadArray(I) = NodeSumLoadArray(I) - C * NodeSumLoadArray(K)
413            Next index
414        Next I
415        
416        NodeSumLoadArray(GetUnknownNodeDisplacementCount) = NodeSumLoadArray(GetUnknownNodeDisplacementCount) / StiffnessMatrixGloble(GetUnknownNodeDisplacementCount, 1)
417        Dim K2 As Integer
418        For K2 = 1 To GetUnknownNodeDisplacementCount - 1
419            I = GetUnknownNodeDisplacementCount - K2
420            JM = K2 + 1
421            If (NW - JM) Then JM = NW
422            Dim TempSum As Double
423            Dim TempIndex As Integer
424            For TempIndex = 2 To JM
425                'NodeSumLoadArray(I) = NodeSumLoadArray(I) -
426                TempSum = StiffnessMatrixGloble(I, TempIndex) * NodeSumLoadArray(I + TempIndex - 1)
427            Next TempIndex
428            NodeSumLoadArray(I) = NodeSumLoadArray(I) - TempSum
429            NodeSumLoadArray(I) = NodeSumLoadArray(I) / StiffnessMatrixGloble(I, 1)
430        Next K2
431        
432    Next K
433End Function
434
435'it works good!
436Public Sub MatrixMultiply(m() As Double, n() As Double, ReturnValue() As Double)
437     Dim I As Long, j As Long, K As Long, row As Long, column As Long, max As Long
438     row = 6
439     column = 1
440     max = 6
441     'ReDim ReturnValue(1 To row, 1 To column)
442     For I = 1 To row
443          For j = 1 To column
444               For K = 1 To max
445                    ReturnValue(I) = ReturnValue(I) + m(I, K) * n(K)
446               Next
447          Next
448     Next
449End Sub
450
451
452Public Function MatrixTranspose(m() As Double, ByRef ReturnValue() As Double)
453     Dim I As Long, j As Long, row As Long, column As Long
454     row = 6
455     column = 6
456     ReDim ReturnValue(1 To column, 1 To row)
457     For I = 1 To row
458          For j = 1 To column
459               ReturnValue(j, I) = m(I, j)
460          Next
461     Next
462End Function
463
464Public Function MatrixChangeSign(ByRef m() As Double)
465    Dim row, column As Long, I As Long, j As Long
466    On Error GoTo ErrStp
467    row = UBound(m, 1)
468    column = UBound(m, 2)
469    
470    GoTo skip
471ErrStp:
472    column = 1
473skip:
474    
475    row = 6
476    For I = 1 To row
477        If column = 1 Then
478            m(I) = -m(I)
479        Else
480            For j = 1 To column
481                m(I, j) = -m(I, j)
482            Next j
483        End If
484    Next I
485End Function
486
487Public Function SetNodeDisplacementNo(NodeNo As Integer, uNo As Integer, vNo As Integer, thitaNo As Single) As Integer
488    NodeDisplacementNoArray(1, NodeNo) = uNo
489    NodeDisplacementNoArray(2, NodeNo) = vNo
490    NodeDisplacementNoArray(3, NodeNo) = thitaNo
491    SetNodeDisplacementNo = NodeNo
492End Function
493
494Public Function GetUnknownNodeDisplacementCount() As Integer
495    GetUnknownNodeDisplacementCount = 3 * NodeCount
496End Function